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 (Hauksson and Shearer, 2006) 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 ω² 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 ts⁢t⁢a⁢t⁢i⁢o⁢n* (Haberland and Rietbrock, 2001):

The ts⁢t⁢a⁢t⁢i⁢o⁢n* 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*=ti⁢j*⁢ts⁢t⁢a⁢t⁢i⁢o⁢n* is often described by the station-dependent κ 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 (Hauksson and Shearer, 2006):

We apply the spectrum-fitting procedure from Nugraha et al. (2010) using a grid search technique to estimate the above-mentioned attenuation parameters in Eq. 6. In evaluatingt^*, 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 Ω₀ 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).

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, ti⁢j* only depends on the Q⁻¹ 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 (Hauksson and Shearer, 2006), 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_pwhich gives the minimum data misfit. As shown in Figure 4A, the Q_p value of 300 gives the minimum root-mean-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.

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.

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 (Hauksson and Shearer, 2006; 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. (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 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).

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 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).