Seismic Attenuation Tomography From 2018 Lombok Earthquakes, Indonesia

Local earthquake data was used to determine a three-dimensional (3D) seismic attenuation structure around the aftershock source region of the 2018 Lombok earthquake in Indonesia. The aftershocks were recorded by 13 seismic stations from August 4 to September 9, 2018. The selected data consist of 6,281 P-wave t∗ values from 914 events, which had good t∗ quality in at least four stations. Our results show that the two aftershock clusters northwest and northeast of Lombok Island have different attenuation characteristics. A low P-wave quality factor (low-Qp), low P-wave velocity (Vp), and high ratio of P-wave velocity and S-wave velocity (Vp/Vs), which coincide with a shallower earthquake (<20 km) northwest of Lombok Island, might be associated with a brittle area of basal and imbricated faults influenced by high fluid content. At the same time, the high-Qp, low Vp, and low Vp/Vs, which coincide with a deeper earthquake (>20 km) northeast of Lombok Island, might be associated with an area that lacks fluid content. The difference in fluid content between the northwest and northeast regions might be the cause of the early generation of aftershocks in the northwest area. The significant earthquake that happened on August 5, 2018, took place in a region with moderate Qp, close to the contrast of high and low-Qp and high Vp, which suggests that the earthquake started in a strong material before triggering the shallower aftershocks occurring in an area affected by fluid content. We also identified an old intrusive body on the northeast flank of the Rinjani volcano, which was characterized by a high-Qp, high-velocity, and a high Bouguer anomaly.

Local earthquake data was used to determine a three-dimensional (3D) seismic attenuation structure around the aftershock source region of the 2018 Lombok earthquake in Indonesia. The aftershocks were recorded by 13 seismic stations from August 4 to September 9, 2018. The selected data consist of 6,281 P-wave t * values from 914 events, which had good t * quality in at least four stations. Our results show that the two aftershock clusters northwest and northeast of Lombok Island have different attenuation characteristics. A low P-wave quality factor (low-Qp), low P-wave velocity (Vp), and high ratio of P-wave velocity and S-wave velocity (Vp/Vs), which coincide with a shallower earthquake (<20 km) northwest of Lombok Island, might be associated with a brittle area of basal and imbricated faults influenced by high fluid content. At the same time, the high-Qp, low Vp, and low Vp/Vs, which coincide with a deeper earthquake (>20 km) northeast of Lombok Island, might be associated with an area that lacks fluid content. The difference in fluid content between the northwest and northeast regions might be the cause of the early generation of aftershocks in the northwest area. The significant earthquake that happened on August 5, 2018, took place in a region with moderate Qp, close to the contrast of high and low-Qp and high Vp, which suggests that the earthquake started in a strong material before triggering the shallower aftershocks occurring in an area affected by fluid content. We also identified an old intrusive body on the northeast flank of the Rinjani volcano, which was characterized by a high-Qp, high-velocity, and a high Bouguer anomaly.

INTRODUCTION
Lombok Island is located in West Nusa Tenggara and is part of the Sunda Arc, Indonesia, which is controlled by the subduction of the Indo-Australian Oceanic Plate beneath the Eurasian Plate with a relative motion of 70 ± 1.0 mm/year to the north (Koulali et al., 2016). The tectonic activity in this region forms the Flores Back-arc Thrust (FBT), which lies north of Flores, Sumbawa, and Bali (Hamilton, 1979;Silver et al., 1983). The term FBT, was first used by Silver et al. (1983) to map the forwardmost thrust to reach the surface above the decollement of the zone of the back-arc thrust north of Flores Island. According to the interpretation of 2D seismic reflection profiles (Silver et al., 1983;Yang et al., 2020), the complicated geological structure in the FBT is expressed by low-angle main basal faults and high-angle imbricate thrusts and its associated faults.
Two major geological structures significantly control the seismicity of Lombok Island (Figure 1): the megathrust in the south due to the subduction of the Indo-Australian Plate beneath the Eurasian Plate (Hamilton, 1979) and the FBT north of Lombok Island which dips to the south (Hamilton, 1979;Silver et al., 1983;Irsyam et al., 2017). The FBT occurs because of the accommodation of energy from the megathrust (Irsyam et al., 2017) which extends in an east-west direction with an estimated convergence rate of around 5.6-6.0 mm/year (Susilo et al., 2018). FBT extends from the north of central Flores Island to the northern part of Sumbawa Island (Hamilton, 1979). However, further research conducted by Silver et al. (1983) using seismic reflection data shows that this thrust extends further west to the Bali Basin.
A series of major earthquakes in 2018 occurred in the northern part of Lombok Island which caused severe damage to the Lombok region. On July 28, 2018, at 22:47 (UTC), the first earthquake (Mw 6.4) struck the island and on August 5, 2018, at 11:46 (UTC) a second, more massive earthquake (Mw 7.0) took place west of the first event with both epicenters located north of Rinjani volcano (Figure 1). These earthquakes were followed by light to moderate aftershocks for several days (≤Mw 5.5); they extended northwest and induced a moderate event (Mw 5.9) on August 9, 2018, at 05:25 (UTC). On August 19, 2018, other two strong earthquakes struck Lombok Island; however, the epicenters were located northeast of the island. The first earthquake occurring on August 19, 2018 (Mw 6.3), took place at 04:10 (UTC), and the second earthquake (Mw 6.9) took place at 14:56 (UTC). After these earthquakes, aftershock activities extended northeast and further south. In the previous study of Sasmi et al. (2020), a total of 3,259 aftershock events have been identified with magnitudes ranging from M w 1.7 to 6.7 and the depth between 5 and 25 km. The focal mechanism solutions from the Global Centroid Moment Tensor revealed that the Lombok earthquakes have a reverse fault mechanism with the strike generally oriented in an east-west direction parallel to the FBT (Sasmi et al., 2020).
The source area of the series of earthquakes occurring north of Lombok Island is an interesting study area as it produced a series of large earthquakes in a short amount of time (less than a month). Some studies have been conducted in this area; e.g., a study of the hypocenter and magnitude analysis of the 2018 Lombok aftershock, which suggests that spatial and temporal clustering of the seismicity may be linked to the segmentation of the Flores oceanic crust (FOC) (Sasmi et al., 2020). A bodywave travel time tomography study using a local network has also been done, which provide information about the characteristics of the mainshock and aftershock source regions (Afif et al., 2020). The recent study by Lythgoe et al. (2021) shows the influence of thermal squeezing in controlling the rupture in the Lombok seismogenic zone.
In this study, we applied seismic attenuation tomography to further investigate the characteristics of the source region. Seismic attenuation is the amplitude decrease of propagating waves due to the inelasticity and the inhomogeneities of the Earth, which are then parameterized using the dimensionless quality factor Q (Eberhart-Phillips, 2002;De Siena et al., 2010). This parameter primarily depends on temperature and fluid content (Karato, 2004;Wang et al., 2017); hence, seismic attenuation tomography is suitable for investigating very heterogeneous regions such as volcano, fault, and subduction zones (Eberhart-Phillips, 2002;Liu and Zhao, 2015;Wang et al., 2017). Several studies on attenuation tomography have been carried out to understand the physical state of fluids in a porous medium and reveal new details about the geometry of faults (Rietbrock, 2001;Hauksson and Shearer, 2006;Chiarabba et al., 2009;Muksin et al., 2013;Liu and Zhao, 2015;Komatsu et al., 2017;Wang et al., 2017). An attenuation tomography study of the southern California crust shows that the transition between high and low attenuation zones coincides with a major Late Quaternary fault . The attenuation study done in the Southwest Japan Arc by Liu and Zhao (2015) shows that a high attenuation anomaly could be identified in or around the active faults. Previous attenuation studies were also successfully carried out to further investigate in detail the source zone of the Mw 7.3 2016 Kumamoto earthquake, Japan (Komatsu et al., 2017;Wang et al., 2017). Wang et al. (2017) estimated the attenuation model of the source zone of this earthquake and found that the mainshock occurred in a low attenuation and high-velocity zone in the upper crust underlain by high attenuation and low-velocity anomalies.
In this study, we estimate the Qp structures beneath the source region of the 2018 Lombok earthquakes using the same local network as Sasmi et al. (2020) and Afif et al. (2020). We aim to reveal the generation of the destructive earthquakes, the causes of spatial clustering of the seismicity in the 2018 Lombok earthquake, and any other geological features that could be identified using attenuation results. The integration of Qp structures and the P-and S-wave velocity structure from Afif et al. (2020) could give a better understanding of the source region.

PREVIOUS GEOLOGICAL AND GEOPHYSICAL STUDIES
The interpretation of the seismic profile north of Lombok Island by Yang et al. (2020) reveals that there are three principal stratigraphic units (from deep to shallow levels): the basement,  (Irsyam et al., 2017;Afif et al., 2020). (C) Map of Lombok Island and the surrounding area which shows the distribution of the initial hypocenters from Afif et al. (2020). The hypocenters plotted here are the selected hypocenters from the spectral fitting process. The colored circle shows the depth of the earthquakes. The blue inverted triangles are the seismic stations, the big red triangle is the Rinjani volcano, the small red triangles are a group of volcanoes that are no longer active, and the yellow stars are the significant earthquakes.
the Upper Miocene limestone, and Pliocene to present-day young sediments. The basement depth is relatively deep in the western part, which implies that the accumulated sediment in the west is thicker than the accumulated sediment in the east. The active basal fault and imbricate fault in the western part are in the sediment layer, while the active fault intersects the basement in the east (Yang et al., 2020). Two strike-slip faults were also identified in the western and eastern parts of Lombok Island, which extend in a south-north direction: the Lombok Strait Strike-Slip Fault in the western part of Lombok Island which crosses the FBT, and the Sumbawa Strait Strike-Slip Fault in the eastern part of Lombok Island (Irsyam et al., 2017).
The currently active Rinjani volcano is located in the northern part of Lombok Island. The Rinjani caldera that can be seen today is the result of the paroxysmal eruption of the Samalas volcano. The AD 1257 Samalas volcanic eruption was one of the most massive, explosive Holocene eruptions and released a large amount of volcanic sulfur into the atmosphere globally (Lavigne et al., 2013). Based on the volcanic activity report from the Center for Volcanology and Geological Hazard Mitigation Indonesia the Rinjani volcano has not shown a significant increase in volcanic activity after this series of earthquakes and the last eruption occurred in 2016 (Global Volcanism Program, 2018).
Several geophysical studies have been conducted around the Lombok Islands. Based on the 1978 Flores earthquake and gravity data, McCaffrey and Nabelek (1984) found evidence of a 30 • dipping southward active thrusting of the FOC beneath Lombok, which also confirmed by seismic velocity tomography of Widiyantoro et al. (2011). From the gravity regional study done on Lombok Island, a high Bouguer gravity anomaly of around 140 mGal was found northeast of the Rinjani volcano, which can be interpreted as the expression of an igneous rock or magmatic conduit (Sukardi, 1979;Zubaidah et al., 2010). Afif et al. (2020) in their recent tomographic velocity study revealed that most of the significant earthquakes occurred at the edge of a high-velocity region where the stress concentration is expected to be the highest. The following aftershock is clustered in a low velocity and high-Vp/Vs zone. Afif et al. (2020) interpreted this zone as a highly fractured fault zone with a large amount of fluid. They suggest that the presence of fluid might accommodate the slip by lubricating and reducing the friction on the fault. After initiating some significant earthquakes, the relaxation of the stress in the subducting oceanic crust caused the generation of further events in areas where fluid infiltration had weakened the crust. They also interpreted a low velocity and high-Vp/Vs region located beneath the Rinjani volcano as caused by the presence of high temperatures and/or magma. A high-velocity and low-Vp/Vs on the eastern flank of Rinjani are interpreted as a cooled intrusive complex supported by a high Bouguer anomaly (Sukardi, 1979).

DATA AND METHOD
In this study, we use a temporary station network installed one week after the first large event of the 2018 Lombok earthquake sequence that occurred on July 29, 2018. Thirteen stations were deployed from August 4 to September 9, 2018, with the main purpose to monitor the aftershocks pattern and to obtain the subsurface image around the source zone. The coordinates of the stations are provided in Supplementary Table 1. We used aftershock data from the previous study done by Afif et al. (2020). We calculated the path-averaged attenuation factor (t * ) as a function of the travel time and dimensionless quality factor using the spectrum-fitting analysis before performing the tomography inversion.

Determining t *
The attenuation structure can be estimated based on the amplitude decay of seismic waves which can be quantified through the t * parameter. The observed amplitude spectrum A ij f from event i to station j can be described as (Scherbaum, 1990;Muksin et al., 2013): where f denotes the frequency, R j (f ) the site response, I j (f ) the instrumental response, and B ij (f ) the absorption spectrum. The far-field source spectrum S i (f ) can be expressed by the following equation (Muksin et al., 2013): where 0i is the long-period spectral level including geometrical spreading factor  and f ci is the source corner frequency. The high-frequency decay factor γ is assumed to be 2; hence, Eq. 3 is equivalent to a Brune type ω 2 source model (Brune, 1970). In terms of the whole path attenuation, the absorption spectrum B ij (f ) is given by Haberland and Rietbrock (2001) as: with t ij denoting travel time. The site response R j (f ) is described by a station correction operator t * station (Haberland and Rietbrock, 2001): The t * station operator is estimated simultaneously with the attenuation inversion as a travel time station correction in velocity inversion tomography (Nugraha et al., 2010). In addition, we set the transfer function of the seismometer to 1 since we only model the observed amplitude spectra in the passband of the seismometer and we are not using the absolute amplitudes (Haberland and Rietbrock, 2001). Using Eqs 1-4, the observed amplitude spectrum A ij f can be expressed as (Eberhart-Phillips, 2002;Muksin et al., 2013): The term t * = t ij * t * station is often described by the stationdependent κ parameter (Anderson and Hough, 1984). Based on numerical tests and observations by Rietbrock (2001) and Shearer et al. (2006), assuming γ is 2 and Q is frequency independent, this provides a good fit for the observed amplitude spectra : We apply the spectrum-fitting procedure from Nugraha et al.
(2010) using a grid search technique to estimate the abovementioned attenuation parameters in Eq. 6. In evaluating t * , the velocity recorded in the time domain is transformed into velocity amplitude spectrum in the frequency domain using Fast Fourier Transform (FFT). The P-wave spectrum is calculated from the vertical component of the seismogram in a window length of 2.56 s (0.56 s before and 2 s after the P-wave arrival) while the noise spectrum was taken 3.06 s before the arrival time of P-waves with the same window length (Figure 2) to calculate the signal-to-noise ratio (SNR). We first determine the average frequency corner f c for each event before calculating t * since there is a trade-off between f c and t * (Scherbaum, 1990). We obtained f c ranging from 0.2 to 13.3 Hz. The 0 and t * are determined for a frequency range of 1-25 Hz by fitting Eq. 6 to the P-wave observed spectrum using previously calculated f c . Two selection criteria were applied to ensure reliable estimation of t * ; each event needs to have at least four t * values and have SNR above 2-decibel (dB). We obtained 6,244 t * data from 914 events for P spectra (Figure 2). Since there are certain event criteria used in the data set, the number of earthquakes used for tomography attenuation differs from the velocity tomography of Afif et al. (2020). The minimum and maximum t * values used in the inversion are 0.001 and 0.1, which are similar to the values in the Southern California attenuation study done by Hauksson and Shearer (2006) and the Northern and Central California study by Lin (2014).

Attenuation Tomography
The relationship between the attenuation operator t * , quality factor Q, and velocity V can be expressed in the following equation: where ds is the distance along the ray path from hypocenter i to station j. Given the a priori velocity model, t * ij only depends on the Q −1 values along the ray path; hence, Eq. 1 is similar to travel time tomography with known origin time and travel times (Haberland and Rietbrock, 2001).
We apply the SIMULPS2000 code by Thurber and Eberhart-Phillips (1999) to obtain a 3D frequency-independent Q p structure from the t * values. The SIMULPS2000 code has been widely used in research on attenuation, such as in studies on the southern California crust , the western Japan subduction zone (Nugraha et al., 2010), the Tarutung geothermal fields (Muksin et al., 2013), and Kilauea volcano (Lin et al., 2015). For the model parameterization, we used a similar grid size as the 3D velocity model of Afif et al. (2020) with a horizontal spacing of 10 km and vertical spacing of 10 km from 0 to 60 km depth (Figure 3). We also use the hypocenter location and origin times determined simultaneously with the 3D velocity inversion of Afif et al. (2020) using Simulps12 (Evans et al., 1994).
Inversion usually depends on the initial models; hence, to obtain a robust Q p model, we searched for an optimal homogenous 1D Q p for the initial model. We ran a series of single-iteration inversions with Q p values varying from 100 to 700 and chose Q p which gives the minimum data misfit. As shown in Figure 4A, the Q p value of 300 gives the minimum rootmean-square (RMS) of data misfit and the RMS varies from 0.01999 to 0.02628 s.
In a tomographic inversion, damping parameters are needed to define how quickly the perturbation of the model can change (Lanza et al., 2020). The optimum values were determined from a trade-off curve between model variance and data variance by running a series of iteration inversions with a broad range of damping values (0.4-200), shown in Figure 4B, similar to the approach applied to the velocity model inversion (Lin, 2014). The inversion results show that the damping parameter of 200 has a very low Q p perturbation value (maximum perturbation of 0.1). At damping values 0.4 to 2, the results of the Q p anomaly pattern are similar, but with different perturbation values. Our results show that the smaller the damping, the stronger the perturbation of Qp so that each pattern shows an increasingly intense color. For Q p inversion here, we choose a damping parameter of 0.8 which produces an optimal compromise between minimizing the data variance versus increasing the model variance in the trade-off curve ( Figure 4B). We also provide the Q p inversion result using a damping parameter of 1 in the Supplementary Figure 1 which shows similar results to the damping parameter of 0.8.

Model Resolution
The quality of the Q model was evaluated using three different methods: Diagonal Resolution Elements (DRE) (Menke, 2018), Derivative Weight Sum (DWS) (Toomey and Foulger, 1989), and Ray Hit Count (RHC), which are calculated during the inversion. The DWS measures the ray coverage of the data by estimating the total ray length that samples a node. The DRE provides a better estimation of the model resolution for each node since it measures the relative contribution of the hypothetical "true" model value to the estimated value. The best solution corresponds to a DRE value of 1; whereas, zero represents no resolution. The quality of model resolution is shown by horizontal and vertical cross-sections which depict well-resolved parts of the study area. Figures 5, 6 show all three-resolution models. DRE for each profile has a range for values 0-0.8, while RHC has a range of values 0-2000, and DWS has a range of values between 0 and 900,000. We determined the well-resolved area quantitatively with the DRE value >0.1 (Lanza et al., 2020) bounded by the dashed line. All the resolution tests show that the area with good quality is located at 0 km down to a depth of 20 km where the largest well-resolved area appears at a depth of 10 km (Figure 5). However, at depth of 30 km, we observed that in the northwest area the DRE, DWS, and RHC still give resolution due to the deeper earthquake depth.

RESULTS AND DISCUSSION
The tomographic attenuation results are shown in Figure 7 for the horizontal map at different depths and Figure 8 for the vertical cross-sections (AA , BB , CC to DD marked in Figure 3A). In these figures, P-wave velocity (Vp) is plotted as percent velocity perturbation relative to the 1D initial velocity FIGURE 2 | Example of spectral amplitude fitting of one earthquake in ten stations (location is shown by inset map). For each plot, we show the normalized time-series waveform data of the vertical component (top panel) and the signal and noise spectra (bottom panel). The window used in the analysis is 2.56 s (0.56 s before the arrival time and 2 s after the arrival time). The signal used in the study is within the blue-shaded windows, and the noise is within the gray-shaded windows. The gray line is the noise spectrum, the black line is the signal spectrum, and the blue line is the model (calculated) spectrum.

FIGURE 3 | (A)
The vertical cross-section positions (AA , BB , CC , and DD lines) relative to grid node; the earthquakes plotted in the cross-sections are located within the areas indicated by the black rectangles (6 km in width); (B) the ray paths on 3D grid node; (C) the map view of ray paths; (D) the east-west cross-section of ray paths; (E) the northeast cross-section of ray paths. The inverted blue triangles indicate the stations; the gray lines indicate the ray paths; the open circles indicate hypocenters; the red plus marks indicate the grid nodes, the red triangle is the active volcano of Rinjani, the black lines are active faults, and the jagged black line is the FBT based on Irsyam et al. (2017). model, Vp/Vs are plotted as absolute values, and Qp is plotted as an absolute value. The blue and red colors represent the higher and lower values, respectively, for the Vp, Vp/Vs, and Qp. The results of the data analysis and the model resolution quality (Figures 5, 6) show a good resolution structure in some regions of Lombok Island and around the Rinjani volcano. From the tomographic results and the hypocenters location, we observe four interesting anomalies (R1-R4 as shown in Figures 7, 8) that we consider to be the most significant with acceptable resolution.
In order to interpret the velocity and attenuation tomography of the island of Lombok, a wide range of physical properties must be considered. In this region, we not only observe the subduction of the Indo-Australian plate in the south, but also the subduction of FOC beneath the Lombok Island which causes back-arc thrust faults in the north. Not far from the thrust faults zone (∼ 50 km), there is also the Rinjani volcano which is still active today. These geological settings cause the crustal conditions in this area to be unique and complex. Besides typical fractured and saturated rocks, there may be zones of partial melt with various melt percentages, accumulations of magma with various sizes, zones of hydrothermal alteration, zones with dense fracture, zones with high temperature, or other possibilities. Each of these crustal conditions has its own effect on seismic wave propagation. In general, a low-Qp together with high Vp/Vs anomaly have been interpreted as the existence of high temperature, partial melting, magma, or increase in fracture concentration in saturated rock (Sanders et al., 1995).
Many studies that examined the relationship between seismic wave properties and temperature at both high and low frequencies show similar results where the velocity and Q of seismic waves decrease with increasing temperature (Mizutani and Kanamori, 1964;Murase and McBirney, 1973;Kampfmann and Berckhemer, 1985;Christensen and Wepfer, 1989). Murase and McBirney (1973) measured rock properties on four igneous rocks (basalt, andesite, rhyolite, and tholeiite) over a range of temperatures including their melting interval at 100 kHz. They found that Vp decreased by about 50% and Qp decreased continuously from about 3000 to less than 10 from the solid temperature to melting temperature. Mizutani and Kanamori (1964) measured the physical properties of metal alloy and also found that Qp was about 90% lower in liquid compared to solid alloy.
The presence of fractures or pores in dry rock causes Vp and Vs to decrease where the reduction in Vp is greater than Vs so that Vp/Vs decreases (O'Connell and Budiansky, 1974). However, in saturated rock conditions, the existence of these fractures and pores will result in a decrease in Vp and Vs where the decrease in Vs is greater, hence, Vp/Vs becomes high (O'Connell and Budiansky, 1974;Moos and Zoback, 1983). Laboratory experiments by Mavko and Nur (1979) show that small amounts of water in fractured rocks can significantly increase attenuation. They also found that the low aspect ratio (i.e., flat) will cause relatively greater attenuation than more spheroidal porosity, therefore fractures have a relatively greater effect than pores.

Attenuation in the Source Area
The aftershock pattern of the 2018 Lombok earthquake could be divided into two clusters: the northwest cluster shallower than 20 km; the events of July 29, August 5, and August 9, 2018, took place in this area; and the northeast cluster which has a depth greater than 20 km and is the area where the August 19, 2018 events took place, as shown in Figure 1. The aftershocks are first concentrated in the northwest cluster (starting from August 4) before moving to the northeast cluster (starting from August 19), as also shown by previous study of Lythgoe et al. (2021).
In this study, we suggest that two factors could play an essential role in the 2018 Lombok earthquake sequence. The first factor that we suggest is the influence of fluid content which is similar to the 2016 Kumamoto, Japan earthquake (Komatsu et al., 2017;Wang et al., 2017) and the Tohoku, Japan earthquake (Liu et al., 2014). Based on experimental studies, Tenthorey et al. (2003) confirm that fluids could significantly affect the frictional strength and the stability of faults. From Figures 7, 8, our results show that the August 5, 2018 earthquake (Mw 7) is located in a region of moderate Qp (∼ 285), close to the contrast of high and low-Qp. This area also corresponds with the high-Vp anomaly. Most of the aftershocks that follow the major earthquake are at a shallower depth of 20 km, which is characterized by low-Qp and low Vp anomaly (Figure 8). Based on these properties, we suggest that the August 5, 2018 earthquake (Mw 7) was generated in a strong material layer of subducted FOC. However, the rupture area, identified by the aftershock distribution which has low-Qp and low Vp, might be affected by fluids, thus reducing fault friction, and could be the main cause of the aftershocks in this area (Tenthorey et al., 2003;Wang et al., 2017). Our attenuation structure findings support the previous seismic velocity study done by Afif et al. (2020) who state that large earthquakes take place in the boundaries of a strong material or at the barriers of the fluid content where stress has accumulated. The second factor that we suggest could play an essential role in this earthquake is the influence of thermal gradient. The recent study of Lythgoe et al. (2021) stated that there is an effect of thermal gradient from the deep magmatic system underlying Rinjani volcano and it controlled the rupture in the seismogenic zone.
We observe a significant low-Qp (∼ 270) anomaly at a depth between 0 and 10 km in the northwest cluster (indicated by R1 in Figures 7, 8). The low-Qp in R1 corresponds with the low-Vp anomaly and high-Vp/Vs (excluded the eastern portion of it). We suggest that at shallow depths (∼ 0 km), the low-Qp is correlated with sedimentary deposits in the Back-arc Basin. Sedimentary deposits are well known to be one of the causes of high attenuation Li et al., 2006;Bohm et al., 2013;Muksin et al., 2013). Bohm et al. (2013) found high attenuation in Kendeng Basin in Central Java. Muksin et al.  Afif et al. (2020), and Q p (right column) at depth of 0, 10, 20, and 30 km. The color bar shows the variations in V p , V p /V s , and Q p . The V p are plotted as percent perturbations relative to the 1D initial velocity model, and V p /V s are plotted as an absolute value. Blue and red colors represent higher and lower values. The black dashed line in the Q p tomogram indicates well-resolved areas based on DRE (>0.1). The yellow stars represent the hypocenters of the significant event; white dots indicate aftershock hypocenters; the red triangle is the Rinjani volcano. The captions listing R1, R2, R3, and R4 are the regions interpreted in this work.
(2013) found high attenuation in the Tarutung Basin, North Sumatra. The thick sedimentary deposit in R1 is formed by upper Miocene limestone overlain by the Pliocene to present-day young sediments with a depth range of 1-5.5 s two-way travel time from mean sea level (Yang et al., 2020). For the deeper depths (∼ 10 km), we suggest that the low-Qp (∼ 250) is associated FIGURE 8 | The vertical cross-sections AA , BB , CC to DD (see Figure 3A), which show the seismic velocity structure for V p (left column), V p /V s ratio (middle column), from Afif et al. (2020) and Q p (right column). The color bar shows the variations in V p , V p /V s , and Q p . V p plotted as percent perturbation relative to the 1D initial velocity model and V p /V s and Q p are plotted as absolute values. The blue and red colors represent higher and lower values. The black dashed line in the Q p tomogram indicates the well-resolved area based on DRE (>0.1). The white dots indicate the locations of hypocenters with a radius of 3 km from the centerline. The yellow stars represent the hypocenters of the significant events. The captions listing R1, R2, R3, and R4 are the regions interpreted in this work. with a brittle area; this is based on the evidence of basal and imbricated thrust fault from the 2D seismic profile north of Lombok Island found by Yang et al. (2020). This strong low-Qp anomaly is also associated with low-Vp anomaly and high-Vp/Vs as well as high seismicity. It corresponds to the anomaly caused by the presence of fractures in saturated rock based on laboratory experiments (O'Connell and Budiansky, 1974;Moos and Zoback, 1983). They found that in saturated rock conditions, the existence of fractures reduce Vp and Vs values where the decrease in Vs is greater, hence, Vp/Vs becomes high (O'Connell and Budiansky, 1974;Moos and Zoback, 1983). The fluid content in fractured rock provides a significant reduction in Q (Mavko and Nur, 1979). Once liberated, these fluids may drain through the active high angle imbricate thrust fault or migrate along the basal detachment fault (Ranero et al., 2008). Alternatively, the fluids could originate from the FOC dehydration beneath the Lombok Islands (Silver et al., 1983). This fluid might act as a driving factor for a smooth fault which then triggers seismicity.
The low-Vp, high-Vp/Vs, and low-Qp anomaly in R1 does not extend further to the east; based on the attenuation results (Figure 7), the anomaly has differing properties between the shallower earthquake cluster northwest and the deeper earthquake cluster to the northeast. In the northeast, the deeper earthquake cluster (indicated as R2 in Figures 7, 8) revealed high-Qp anomalies that correspond with low-Vp and low-Vp/Vs from ∼10 km onward. The low-Vp/Vs can be caused by fractures in the dry rock. The fractures reduce both Vp and Vs, but Vp is reduced more than Vs so that Vp/Vs has a low value (O'Connell and Budiansky, 1974). The high-Qp also supports the idea that fluids do not play a role in this zone. Thus, we suggest that the difference in Vp/Vs and the attenuation structure between R1 and R2 is mainly due to differences in the fluid content. Besides, we also argue that the difference in the fluid content between the northwest and northeast is causing the generation of aftershocks, primarily in the northwest. The presence of fluid can cause pressure to move first in the northwest area, before triggering earthquakes in the northeast due to the differences in pore pressure. Alternatively, the differences in attenuation characteristics and aftershocks distribution in R1 and R2 are probably due to the thermal gradient associated with the presence of the Rinjani volcano as stated in Lythgoe et al. (2021). The presence of a squeezing zone generated by a thermal push from the volcanic heat in the R1 might raise the seismogenic zone to the shallower crust so that the aftershocks in R1 (northwest area) are shallower than the aftershocks in R2 (northeast area).

Other Geological Features
The area north of the Rinjani volcano (indicated as R3 in Figures 7-9) at a depth of 0-15 km has a low-Vp, moderate to high-Vp/Vs, and low-Qp (∼ 230) at the depth ∼ 0-20 km. The Q value in this area is the lowest of all the areas. Seismicity in this area is relatively low and does not relate to clustered  Irsyam et al. (2017). The yellow star denotes the M 7.0 (August 5, 2018) which has a thrusting fault mechanism (Global CMT Catalogue). The white dots indicate the locations of hypocenters with a radius of 3 km from the centerline. The red lines indicate the faults from Supendi et al. (2020) and the black lines indicate the faults from Yang et al. (2020). The green area near the red lines (faults) indicates the high-temperature fluid content. The dark blue sediment layer represents the UpperMiocene limestone, the green sediment layer represents the Pliocene to present young sediment, and the light blue layer represents the ocean. Note that the difference in the number of hypocenters is due to the difference in the number of events before the spectral fitting procedure (A) and after spectral fitting procedure (B,C). earthquakes in the north of the island (R1). Based on the previous studies by Supendi et al. (2020) and Yang et al. (2020) this area is composed of intensely deformed rock containing dense fractures. These dense fractures cause both Vp and Vs to be low (O'Connell and Budiansky, 1974). The low-Qp and high-Vp/Vs suggest that this area also contains fluids. In saturated conditions, the decrease in Vs is greater than Vp, so that Vp/Vs becomes high (O'Connell and Budiansky, 1974;Moos and Zoback, 1983). The small amounts of fluids distributed in thin fractures can also significantly cause Q to decrease (Mavko and Nur, 1979). In addition, we suggest that the high temperature also plays a role in causing this area to have the lowest Q. The indication of the presence of the high temperature is supported by the existence of thermal manifestations and surface alteration in the northern part of the Rinjani volcano (Sundhoro et al., 2000;Hadi et al., 2007). Thermal manifestations that appear on the surface are hot springs with temperatures up to 70 • C (i.e., Sebau, Kukusan, Kalak, Segara Anakan I, Segara Anakan II, Segara Anakan III, and Segara Anakan IV) (Sundhoro et al., 2000;Hadi et al., 2007). Based on the velocity and attenuation characteristics and also an indication of high temperature, this region might contain fluids: either high temperature water-filled or melts (Nakajima and Hasegawa, 2006;Ramdhan et al., 2019). We suggest that the dense fractures at a depth of ∼ 0-20 km with high permeability may be the pathway for hot water to rise to the surface.
R4 is located in the eastern flank of the Rinjani volcano and is characterized by a high-Vp, low-Vp/Vs, and high-Qp anomaly at a depth of 0-10 km (Figures 7, 8). This region also correlates with a high Bouguer anomaly (Sukardi, 1979) and a high-Vs anomaly based on ambient noise tomography (Sarjan et al., unpublished). Indrastuti et al. (2019) found a similar velocity anomaly beneath the Sinabung volcano using data from before the 2013 eruption; they suggest that the anomaly is related to old, cold intrusive deposits from the previous eruption. Other previous studies that found high-Vp anomalies around volcanoes also interpreted these as older intrusive bodies (Laigle et al., 2000;Tanaka et al., 2002;Syracuse et al., 2011). In line with this study, we interpret the R4 anomaly as the seismic expression of an old intrusive body which relates to the ancient Samalas volcano on the northeast flank of the Rinjani volcano. According to the geological map of the Sembalun Bumbung area (Sundhoro et al., 2000;Hadi et al., 2007), the modern presence of the Samalas volcano is marked by a group of volcanoes that are no longer active in this region (i.e., the Gonduri, Pusuk, Nangi, Anakdare, and Batujang volcano) (Figure 1).

CONCLUSION
We successfully estimate the 3D Qp structures around the source region of the 2018 Lombok earthquake sequence. The interpretations of the Qp structures in this study are supported by the Vp and Vp/Vs structures from a previous study (Afif et al., 2020) and the 2D seismic profiles north of Lombok Island (Yang et al., 2020). The main findings of this study are summarized as follows: 1. The clustered northwest and northeast earthquakes have different velocity and attenuation properties. The low-Qp supported by low-Vp and high-Vp/Vs in the northwest area might be associated with the high fluid saturation in this area. Conversely, the high-Qp, low-Vp, and low-Vp/Vs in the northeast area might be associated with the lack of fluids. We also suggest that the existence of fluid content might be the cause of aftershocks that happened in the northwest area before moving northeast. 2. The significant Lombok earthquake on August 5, 2018 is located in a moderate-Qp and associated with a high-Vp zone. This feature indicates that the Lombok earthquake sequence first ruptured in a strong material of the subducted FOC and triggered shallower aftershocks in areas with fluid contents. 3. The low-Qp, low-Vp, and high-Vp/Vs anomaly beneath the northern area of the Rinjani volcano indicates hightemperature fluids. 4. The high-Qp, high-Vp, and low-Vp/Vs anomaly at the northeast flank of the Rinjani volcano is considered to be an old intrusive body related to the ancient Samalas volcano.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
ZZ, AN, ShW, DS, and MM conceived the survey and seismic study on AP, AN, AAu, BP, SR, HA, contributed to the writing of the manuscript. All authors contributed to the preparation of the manuscript. All authors have read and approved the final manuscript.

ACKNOWLEDGMENTS
This study was a collaboration between our colleagues at the Institut Teknologi Bandung (ITB), the Earth Observatory