Research ArticleGEOPHYSICS

Thermal squeezing of the seismogenic zone controlled rupture of the volcano-rooted Flores Thrust

See allHide authors and affiliations

Science Advances  29 Jan 2021:
Vol. 7, no. 5, eabe2348
DOI: 10.1126/sciadv.abe2348


Temperature plays a critical role in defining the seismogenic zone, the area of the crust where earthquakes most commonly occur; however, thermal controls on fault ruptures are rarely observed directly. We used a rapidly deployed seismic array to monitor an unusual earthquake cascade in 2018 at Lombok, Indonesia, during which two magnitude 6.9 earthquakes with surprisingly different rupture characteristics nucleated beneath an active arc volcano. The thermal imprint of the volcano on the fault elevated the base of the seismogenic zone beneath the volcanic edifice by 8 km, while also reducing its width. This thermal “squeezing” directly controlled the location, directivity, dynamics, and magnitude of the earthquake cascade. Earthquake segmentation due to thermal structure can occur where strong temperature gradients exist on a fault.


The maximum magnitude of an earthquake on a fault is largely determined by the area of the fault surface located above the transition from brittle to visco-plastic deformation, a region known as the seismogenic zone (1, 2). Thus, gently dipping megathrusts produce significantly larger earthquakes than vertical strike-slip faults. The brittle-to-ductile transition (BDT) that defines the base of the seismogenic zone is highly sensitive to temperature, and the thermal structure of the crust generally defines the maximum seismogenic depth (1). An additional restriction on earthquake size arises from geometric features of complex fault surfaces, which can limit both the lateral and up-dip propagation of ruptures (3, 4). Finally, a general dearth of seismicity in the very shallow continental crust indicates a different, and possibly temperature-related, limitation on shallow faulting (5, 6).

Tectonic faults that develop within the thermally heterogeneous zone of an active continental volcanic arc can offer a unique window into the role of temperature in controlling seismicity. In such settings, geometrically simple faults that traverse the whole crust can be perturbed by strong lateral thermal gradients from large magmatic plumbing systems of isolated arc volcanoes. However, earthquakes in this tectonic setting are rare and generally not well monitored. In 2018, an unusual and destructive cascade of earthquakes occurred over a 3-week period at the island of Lombok, Indonesia (7, 8). These earthquakes were produced by cascading rupture of the Flores Thrust, which extends beneath the active Sunda volcanic arc (9). At Lombok, the south-dipping Flores Thrust intersects the root zone of the great Gunung Rinjani volcano and, at this longitude, slips at ~10 mm/year (10). Gunung Rinjani is a spatially isolated arc volcano that hosted one of the largest and most explosive [Volcanic Explosivity Index (VEI) 7] eruptions in the Holocene (11).

The 2018 sequence began with a moment magnitude (Mw) 6.4 earthquake on 28 July, followed by a Mw 6.9 event to the west on 5 August (hereby called 5A), and a Mw 6.2 and another Mw 6.9 (called 19A) to the east on 19 August (Fig. 1). In response to the 28 July foreshock, we rapidly deployed a temporary seismic array to Lombok, which became operational just before the first Mw 6.9 earthquake and was deployed for 34 days (fig. S1). Here, we combine the array observations with data from regional stations, teleseismic waveforms, and space geodetic measurements of surface displacement to characterize the rupture sequence (12). We determine precise hypocenters of the large events and their aftershocks using double-difference relocation. We invert for the best double couple focal mechanism for all Mw >5 earthquakes using regional broadband records, and independently constrain the origin depth by modeling broadband teleseismic P waves containing the depth phases. We derive coseismic rupture models by joint inversion of teleseismic waveforms and Interferometric Synthetic Aperture Radar (InSAR) data.

Fig. 1 Overview of the 2018 Lombok earthquake sequence.

(A) Timeline of seismicity, with red symbols showing Mw >6 events, orange symbols showing Mw >5 events, and other symbols showing aftershocks. Unfilled symbols are events that occurred before the local array was active and were relocated using regional stations. Triangles indicate events located more than 4 km above or below the 3D fault surface. (B) Map of Lombok showing seismicity relocated using a rapidly deployed array, focal mechanisms of large events derived from waveform inversion of regional broadband seismic stations, and coseismic slip distributions from joint inversion of geodetic and teleseismic data. Isotherms on the 3D surface of the Flores Thrust are shown as dashed lines. (C) Cross section of seismicity and coseismic slip, projected onto a vertical plane (A and B). Isotherms on the 3D fault surface are dashed lines.


Together, the aftershocks and coseismic slip patches define a narrow, middle to lower crustal seismogenic zone that varies laterally in depth and width and generally underlies the northern coast of Lombok (Fig. 1). The maximum depth of aftershocks gradually increases along strike away from the volcanic center, from ~15-km depth directly north of Rinjani’s summit crater, to a maximum depth of ~23 km between Lombok and Sumbawa (Fig. 1B). Similarly, the minimum depth of seismicity varies from ~10 km deep north of the summit crater to ~18 km deep to the east. The coseismic rupture models show that the two Mw 6.9 ruptures propagated unilaterally but in opposite directions along a ~25° southward dipping fault (Fig. 1).

While the epicenters of the Mw 6+ earthquakes are all located within 20 km of each other and have similar focal mechanisms (Fig. 1), the teleseismic records of the two Mw 6.9 events exhibit surprisingly different waveform complexity (Fig. 2), reflecting an inherently different rupture style. Compared with the single asperity of the 19A event, the 5A event has more than five smaller asperities, corresponding to rougher moment rate function and a longer duration. The 5A rupture propagated ~50 km to the west within a narrow depth range (~7 km), with a high aspect ratio and multiple major asperities. In contrast, the 19A rupture propagated unilaterally ~25 km eastward, rupturing to the depth range of 8 to 20 km with a smoother and broader distribution of slip characterized by a single major asperity.

Fig. 2 Rupture model of the two Mw 6.9 earthquakes from finite fault inversion (5A and 19A events).

(A) Depth profile of coseismic rupture constructed from joint inversion of teleseismic and InSAR data. The red stars indicate the hypocenters, and the contour lines display the iso-rupture time. Blue vectors indicate rake angles. Multiple asperities (a1–a5) are labeled for the 5A rupture. (B) Observed (black) and synthetic (red) teleseismic displacement waveforms for stations at representative azimuths. Station distance (lower) and azimuth (upper) in degree are displayed at the beginning of each waveform pair. Peak displacement of the data in micrometer is shown at the upper right. (C) Moment rate functions.

The confinement of coseismic slip areas and aftershocks to a narrow belt within the middle to the lower crust, and the strong lateral directivity of the ruptures suggest that upward propagation of rupture was limited in both cases by an unknown factor. Similar behavior during previous large thrust earthquakes that did not reach the surface was due to the existence of geometric ramp-flat structures [e.g., (3)]. To test whether the spatially restricted rupture pattern at Lombok is due to a geometric feature of the Flores Thrust, we constructed a three-dimensional (3D) fault model that combines the relocated hypocenters with constraints on shallow fault geometry from offshore seismic reflection data. The best-fitting fault is an ~25° dipping ramp thrust extending from ~6-km depth to at least 25-km depth over the whole study area, consistent with the dip of focal mechanism nodal planes, and rolls over into to an ~5° dipping décollement underlying the basin sediments of the Bali Sea at depths of ~3 to 6 km (Fig. 3). The data are well fit by a semiplanar structure, indicating that all ruptures were hosted on and terminated upward along a relatively planar fault surface.

Fig. 3 3D perspective diagrams showing the fault geometry, inferred thermal structure of the crust, and the seismicity distribution.

Relocated seismicity is colored by hypocentral depth. The thermal structure of the crust of Lombok, and the Flores Thrust, is indicated using the 80°, 250°, and 450°C isotherms. (A) View from the southwest, showing the region of middle crustal seismicity. (B) View from the southeast, showing the lateral transition from midcrustal to lower-crustal seismicity.

The restricted seismogenic zone at Lombok is instead best explained by thermal perturbation of the Flores Thrust by the deep magmatic system underlying Gunung Rinjani. Thermal control is implied by the spatial correlation between the raised base of seismicity and the volcanic center (Fig. 3). Strong lateral variation in the depth of seismicity is also seen in volcanic regions of Japan, with up to 10 km in relief of the seismogenic zone over horizontal distances of ~20 to 50 km from volcanoes (6, 1316). While no direct constraints are available for the temperature structure of the crust at Lombok, a reasonable thermal model can still be constructed. The temperature of the BDT in other active continental volcanic arc settings is ~450°C (6, 17). Historical seismicity around Lombok indicates a maximum regional seismogenic depth of 25 km (18). Assuming BDT temperature of 450°C therefore requires a cold regional geothermal gradient of ~18°C/km, which is consistent with regional geothermal gradient estimates from the accretionary southeastern boundary of the Sunda plate that underlies Lombok (19). A quartz-rich crustal rheology with a BDT temperature of 300°C (20, 21) would require an extremely cold crust that we deem implausible, so we prefer a value of 450°C that is consistent with a feldspar-rich composition (22). We therefore constructed a thermal model that has a maximum geothermal gradient at Rinjani’s crater, is radially symmetric about the summit crater, and tapers smoothly to the regional geothermal gradient of 18°C/km. Values for the maximum geothermal gradient at the crater and the rate at which the gradient decreases to the regional value are selected to fit the base of 95% of seismicity (D95) to the expected BDT temperature of 450°C (fig. S13).

The seismicity distribution that is spread out in depth (Figs. 1C and 3) becomes clustered in temperature space (Fig. 4). The seismogenic zone is defined by temperature, and not by the layered geological structure of the crust that would presumably not be correlated with Gunung Rinjani’s edifice. There is a close relationship between the distribution of seismicity and coseismic slip, with large earthquakes nucleating near the base of the seismogenic zone. Notably, a proportion of coseismic slip of the 5A event extends below the base of seismicity (Figs. 1C and 4B), indicating that some areas can rupture with dynamic weakening during large events but may be too hot to host aftershocks.

Fig. 4 Variation of the seismogenic zone with depth, temperature, and distance from the volcano, for an assumed 450°C BDT (brittle-ductile transition).

(A) Seismicity extends from the middle (10 km) to the lower crust (25 km). (B) Temperatures of hypocenters and coseismic slip areas cluster between ~250° and 450°C.

The shallow limit of the seismogenic zone is closely approximated by the ~250°C isotherm, although the exact value is dependent on the assumed BDT temperature (Figs. 3 and 4B). Below ~250°C, phyllosilicate minerals present in the fault zone (23, 24) are thought to play a major role in limiting seismogenesis (5, 25, 26). The coincidence of the upper limit of the seismogenic zone with an isotherm of the fault surface supports mechanisms that would promote coseismic slip on the Flores Thrust only at temperatures above 250°C. However, the lack of shallower coseismic slip during this earthquake sequence does not preclude the possibility of tsunamigenic rupture of the shallower ramp fault. The shallow regions of many subduction megathrust faults are largely aseismic, and coseismic slip commonly terminates at substantial depth; however, shallow coseismic slip sometimes occurs in earthquakes that nucleate along the shallow fault (27, 28). Stressing of the shallow ramp fault by deep slip during the 2018 cascade might have enhanced the likelihood of a large event on the shallow Flores Thrust, which may only rupture when exposed to higher loading rates.


In the absence of Gunung Rinjani, the seismogenic zone of the Flores Thrust at Lombok would be wider and deeper within the crust, similar to the far-field (Fig. 5). Volcanic heat has raised the shallow limit of the seismogenic zone into crust that would not normally be seismogenic in this region, possibly by thermal decomposition of phyllosilicates. This shallower seismogenic zone hosted the complex 5A rupture, showing that this area of the fault is more heterogeneous in friction or geometry than the deeper and likely longer-lived seismogenic zone in the east. The lateral segmentation of the 2018 rupture cascade was defined by the limited extent of the first large foreshock, which occurred within the elevated part of the seismogenic zone that is characterized by small regions of intense slip (asperities), where lateral rupture propagation is difficult. The subsequent Mw 6.9 events completed the rupture of the whole seismogenic zone via bilateral propagation away from this area of large rupture nucleation. While regional earthquake records in California suggest that large ruptures may propagate toward areas of raised BDT (29), our observations at Lombok clearly show the opposite behavior. The differences may be due to the poorer resolution in earthquake locations and the lower relative heat flow in California versus the large lateral temperature gradient that exists at Lombok.

Fig. 5 Conceptual illustration of thermal “squeezing” of the seismogenic zone of the Flores Thrust due to volcanic heat from Gunung Rinjani.

Where the fault has been heated near Gunung Rinjani, multiple slip asperities are concentrated within a shallow, narrow seismogenic zone. Ruptures (stars) nucleate near the base of this narrowed zone and propagate laterally, mainly within the confines of the seismogenic zone but extending slightly upward and downward. Farther from the heat source, the “typical” seismogenic zone of the crust is wider and significantly deeper, and ruptures can be larger and smoother.

Thermal effects on the seismogenic zone are best revealed where strong lateral temperature gradients are imprinted onto a geometrically simple, active fault plane. This uncommon geological situation occurs at Lombok, where a relatively cold continental crust hosts both a crustal-scale thrust fault and an isolated, large arc volcano. Thermal “squeezing” of the seismogenic zone at Lombok therefore has exerted first-order controls on the distribution, segmentation, and rupture character of the earthquake cascade.


We conduct a comprehensive analysis of the earthquake sequence and the Flores Thrust at Lombok using seismic data from local, regional, and teleseismic networks, in addition to seismic reflection data, InSAR data, and thermal modeling. In the following sections, we describe each method used in our analysis of the earthquake sequence.

Temporary Lombok seismic array and regional network

We deployed a temporary array on Lombok following the first Mw 6.4 earthquake on 28 July 2018. The array was active from 5 August 2018 to 7 September 2018; hence, it recorded both Mw 6.9 earthquakes, a Mw 6.2 event, and their aftershocks sequences. The local array was composed of 16 three-component stations: 9 temporary broadband seismometers, 6 temporary short-period Z-land nodes, and 1 additional permanent broadband station. Short period nodes were deployed directly into the ground and covered with several centimeters of soil so as to ensure continuous GPS readings for clock drift correction. The temporary array was jointly deployed by the Earth Observatory of Singapore and Bandung Institute of Technology, while the permanent station on Lombok is operated by the Agency for Meteorology, Climatology, and Geophysics of Indonesia (BMKG). We also use regional data from the Indonesia nationwide network operated by BMKG.

Earthquake relocations

The temporary local array and nearby regional stations were used to pick seismic arrival times for earthquake relocations; in total, 22 stations were used (fig. S1). We use the regional earthquake catalog produced by BMKG to obtain event detections and origin times for events with Mw greater than 2. We first pick arrival times manually and obtain initial event locations by a 3D grid search method. During manual picking, arrival time picks are rated on quality to produce a high-quality travel time database. We subsequently refine picks by cross-correlation of P and S waveforms. Cross-correlation lag times are obtained for waveforms bandpass filtered between 0.5 and 4 Hz, windowed 0.2 s before and 3 s after the P and S picks, with cross-correlation coefficients greater than 0.8 retained.

Initial event locations are obtained using a grid search technique, which finds the best-fitting location such that the difference between observed and predicted travel times is minimized (e.g., (30)]. Events are then relocated using a double difference relative relocation method, which minimizes differential travel times between event pairs (31). Travel time picks and cross-correlation lag times are inverted simultaneously, with weighting of cross-correlation lags increasing with each iteration. Cross-correlations are weighted by the square of the cross-correlation coefficient. We use damping values of 80 and 90, which achieve a system condition number less than 50. We begin with 1160 events from 28 July to 7 September and obtained accurate relocations for 719 events. The 1D velocity model used is the best-fitting 1D P wave velocity model obtained from tomography at Lombok using the earthquake sequence, and we note that by using a double difference method, we reduce the dependence on the velocity model compared with absolute location techniques. Events that occurred before our temporary Lombok network was active are relocated using travel time picks from nearby regional stations only; we therefore calculate a second set of relocations using a lower station number requirement to determine locations for these events.

Location uncertainties from the inversion algorithm are less than 100 m, although we note that these are unreliable and do not account for the uneven station and event distributions (31). We therefore obtain a statistical uncertainty estimate by using bootstrap resampling, whereby a random 5% of events are removed and relocations are obtained; this resampling is repeated until all events are missing from one relocated dataset. Uncertainties are taken to be the largest distance from the ensemble of relocated points to the event location determined using all available data. The average uncertainty is 0.8 km east-west, 1.5 km north-south, and 1.5 km in depth. We also extensively benchmark earthquake depths calculated from event relocation with independent depth measurements obtained by waveform inversion of teleseismic waveforms containing the depth phase, as discussed in the next section.

Moment tensor waveform inversion and earthquake depth verification

We obtain focal mechanisms for all events with Mw greater than 5 by inverting regional waveforms. Of the 33 earthquakes with Mw greater than 5 between 28 July and 31 March 2019, we obtain robust focal mechanism solutions for 20 of them. We use waveforms recorded at over 60 regional stations; these stations are primarily from the Indonesia national network, supplemented with stations from the Global Seismic Network (GSN) to improve the azimuthal coverage. Waveform inversion is conducted by using the cut-and-paste (CAP) method, which cuts the seismograms into Pnl and surface wave segments and allows different time shifts between data and synthetics, reducing sensitivity to the velocity model (32). Before waveform inversion, we remove the instrument response and rotate horizontal components into the earthquake reference frame (radial and tangential components). The 1D velocity model used is extracted from a regional P wave model (33). Figure S2 shows the waveform inversion results for a Mw 5.9 event on 9 August 2018, with regional stations providing good azimuthal coverage to robustly constrain the focal mechanism.

To obtain accurate centroid depths, we then use the best focal mechanisms obtained previously and invert high-frequency teleseismic waveforms for the best-fitting depth. We use teleseismic waveforms to isolate the direct P wave and depth phases and therefore maximize sensitivity to centroid depth. We use only carefully selected high-quality broadband vertical seismograms. Figure S2 (C and D) shows that by using high-quality teleseismic data, a clear best-fit depth can be identified with minimal trade-off between depth and source duration. Figure S3 shows observed and modeled P waveforms for several Lombok earthquakes. For earthquakes with magnitudes less than Mw 6, we use high-frequency velocity data, while for earthquakes with magnitude larger than Mw 6, we use broadband displacement data, because at high frequencies, the waveforms are more complex for large earthquakes.

We compared centroid depths and focal mechanisms obtained by our waveform inversion, double difference relocation, and from the Global Centroid Moment Tensor (GCMT) catalog (fig. S4). Overall, earthquake depths between double difference relocation and waveform inversion are very consistent, with differences ranging between 0 and 3 km. There are larger differences between our results and the GCMT catalog, which is expected because our waveform modeling is conducted at higher frequencies than used in GCMT.

Seismic reflection data

We use marine multichannel seismic reflection profiles acquired by MultiClient Geophysical to characterize the shallow part of the Flores Thrust. The 2D seismic profiles were acquired in December 2014 using a 7-km-long streamer at 7.5-m water depth and an airgun source towed at 6.5-m water depth with a shot point interval of 25 m. The data were processed using a standard marine processing workflow, including de-noise, velocity analysis, and prestack time migration. Further-interpreted images of the seismic reflection data can be found in Yang et al. (34).

We interpret the data in the time domain and then depth convert specific horizons using the average velocity above the horizon, which is calculated from seismic migration velocities. We choose to depth convert specific horizons to reduce artifacts from the uncertain velocity model, which is obtained from seismic processing. We manually picked a décollement surface that separates highly deformed sedimentary layers on top from underlying reflectors with greater continuity. The décollement dips to the south toward Lombok and does not reach the surface to the north (fig. S5). The dip of the décollement increases beneath the higher topography of the Lombok shelf. We interpret this décollement as the shallow part of the Flores Thrust, which we use to construct a structural model of the fault.

Fault geometry model of Flores Thrust at Lombok

We construct a geometric model of the Flores Thrust using the relocated hypocenter catalog as an indication of fault location at depth, and a collection of points generated by spatial interpolation of the picked décollement from seismic reflection data (fig. S6). We first fit a 3D polynomial surface to all data points, resulting in a general estimate of the fault location. We exclude hypocenters located more than 4 km above or below this surface. The remaining data are binned into a series of evenly spaced (2.9 km apart), overlapping (extending 10 km east and west of the center line), north-south–oriented swaths. For each swath profile, we calculated the best fit to the data of a piecewise linear functionZ(y)=(y<b1):b2×(yb1)+b4(y>b1):b3×(yb1)+b4(1)

In this form, b1 represents the y coordinate of the breakpoint separating the ramp and flat faults, b4 represents the depth of the breakpoint, b2 is the slope of the steep ramp fault, and b3 is the slope of the shallow décollement. A nonlinear programming solver is used with initial estimates b1 = 9.103e + 06, b2 = 0.45, b3 = 0.05, and b4 = −6189. The slope of the resulting best-fit lines (b2) is consistent with the average dip of the fault plane ruptured during the large thrust earthquakes, as recorded by focal mechanisms. While the bilinear functional form has an unrealistically sharp fault bend, it is sufficient in this case because few data constrain the shallow part of the ramp fault.

InSAR data for the sequence

We processed three pairs of Sentinel-1 Synthetic Aperture Radar (SAR) images from the European Space Agency (ESA) acquired in descending track 32, and three pairs of ALOS-2 SAR images from the Japan Aerospace Exploration Agency (JAXA) acquired in ascending track 129 to map the surface displacement of the earthquake sequence (table S1). All of the Sentinel-1 images were acquired using the Terrain Observation by Progressive Scan (TOPS) mode. The ALOS-2 images were acquired in the StripMap mode.

From the Sentinel-1A/B TOPS images, we extracted all bursts covering Lombok island and processed each burst independently with a preseismic image as the reference. Interferometric image pairs were at first coregistered using a geometric approach with precise orbit ephemerides and the Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM). Then, an interferogram of each imaging burst was formed with its topographic component removed using the same SRTM DEM. All the interferograms were multilooked by a factor of 23 in range and 6 in azimuth, resulting in a ground pixel size of ~90 m. Burst-based interferograms were merged, with an azimuth phase ramp correction using the burst-overlap interferometry (35). The ALOS-2 images were coregistered using a similar geometric method, and interferograms were multilooked to about 90-m ground resolution as well. All interferograms were unwrapped using Statistical-Cost, the Network-Flow Algorithm for Phase Unwrapping (SNAPHU) (36).

We did not use all the possible InSAR images [such as in Salman et al., (7)], but select one ascending (ALOS-2) and one descending (Sentinel-1) InSAR image for each earthquake, and because these satellites have similar-looking angles to the ground, ascending/descending images from ALOS-2 is similar to that from Sentinel-1 [see more details in Salman et al., (7)]. For the slip model inversion, derived displacement measurements were downsampled using the quadtree structures (37). During the downsampling, we took the median value in each downsampling grid, further reducing the influence of outliers and possible localized unwrapping errors. For each downsampled data point, we calculated the look vector based on its geolocation and satellite orbital information.

Finite fault modeling

To derive finite source kinematic models, we use the method developed by (38), which allows the joint inversion of seismic waveforms and geodetic data. Since teleseismic and static data provide complementary constraints on the kinematic rupture process, it helps to suppress trade-off among model parameters. For each subfault, we solve for the slip amplitude and direction, rise time, and rupture velocity. For each parameter, we specify the bounds and a discretization interval. We solve for quadratic ramps in the InSAR data to correct for orbital errors that have not been removed through baseline reestimation and interseismic deformation.

We define the best-fit model as having the lowest objective function, given asEwf+WI×EI+WS×S+WM×M(2)where Ewf is the waveform misfit. EI is the geodetic misfit, S is a normalized, second derivative of slip between adjacent patches (smoothing), M is a normalized seismic moment, and WI, WS, and WM are the relative weighting applied to the geodic misfit, smoothing, and moment, respectively. The least squares misfits are calculated for the teleseismic and geodetic data. Here, we test different values of WI, and we found that setting the weight for the geodetic misfits twice larger than the waveform misfits did not significantly degrade the fits to the teleseismic or geodetic data between the individual and joint inversions given the normalizations schemes. The static green’s functions at free surface are calculated by using the same 1D velocity model from Crust1.0 (39). We use a simulated annealing algorithm (38) to find the best-fitting model parameters for the joint inversions for coseismic slip. This nonlinear, iterative inversion algorithm is designed to avoid local minima by searching broadly through parameter space in initial steps and then, in later iterations, to focus on regions that well fit the data.

Here, the subfault size is chosen as 2 km × 2 km, the rake angles have been constrained to be between 120° and 60°, and the rupture velocity is allowed to vary between 1.5 and 2.5 km/s; the rise time is assumed to be an arc of the cosine function, first quarter, and to vary between 0.4 and 9.0 s with 0.6-s steps. The slip amplitude can change from 0 to 3 m. The depth profiles of the slip models of the A5 and A19 events are presented in Fig. 2, along with the corresponding moment rate functions. The waveform and InSAR data fits can be found in figs. S7 to S10. Excellent fits between the data and synthetics indicate the robustness of the inversions. The digital format of the rupture model could be found in the online supplement material.

To further verify the resolution of the inversion, we also conducted checkerboard tests for the two inversion setups. We generated synthetic data for checkerboard-like slip models (fig. S11) and inverted these synthetic data with the same inversion parameters as we applied to the real data. The results (fig. S11) indicate that the best resolution of the model is located around the coastal line, where the largest slips of the two Mw 6.9 earthquakes are located, which also overlap with most of the aftershocks.

Thermal modeling

We adopt a simple thermal model of the crust parameterized as a spatially variable, linear geothermal gradient. The regional geothermal gradient is assumed to be a constant value. A radially symmetric geothermal gradient field, which is centered on the caldera of Gunung Rinjani, is added to this regional background. The temperature T(Z, D) at a given depth (Z) and distance from the caldera (D) is therefore expressed asT (Z,D)=Z (erf(Dd2)(ggmin)+g)+T0(3)where g is the geothermal gradient at Rinjani’s caldera, gmin is the regional geothermal gradient, T0 is the temperature at the surface of Earth (taken as 4°C because of the presence of deep water above much of the study area), and d is the horizontal distance over which the decrease in geothermal gradient from g to gmin occurs. This simplistic model form is supported by regional observations of active volcanic arcs, which show arching of both the seismogenic depth and the Curie temperature isotherm over wavelengths similar to that of volcanic topography (6). We assume that the seismogenic depth generally corresponds with the 450°C isotherm, consistent with high-temperature brittle-ductile transitions generally observed in active volcanic arcs (6). This assumption requires the far-field geothermal gradient gmin = 18°C/km, to explain the deepest observed seismicity.

The range of geothermal gradients adopted in this study is consistent with a regional heat flow map of Southeast Asia (fig. S12). For sites in the source data that have heat flow data but no geothermal gradient data (40), we estimated the geothermal gradient using a mean thermal conductivity inferred from a linear regression of all sites containing both data types, i.e., g = k × Q, where Q is heat flow and k is the mean thermal conductivity. There is a clear signature of low heat flow in southwest Sulawesi that extends south-westward to Lombok, following the inferred trend of the basement geology, namely, the eastern edge of the East Java-West Sulawesi Block (19). The crust in this belt likely consists of metamorphic rock from accreted oceanic crust and sediment and is therefore expected to have lower radiogenic heat production than the hotter, more granitic terranes of Sunda.

In the absence of local heat flow data, we determine values for the maximum geothermal gradient (g) and the distance over which the geothermal gradient decays to background levels (d) by conducting a grid search over the expected range of values. Motivated by the D95 description of seismogenic depth, we select the suite of models that have 95% of hypocenters located at temperatures less than the 450°C, given the expected isotherm at the base of the seismogenic zone (fig. S13). There is a direct trade-off between the parameters g and d (fig. S13). Using the southern limit of relocated seismicity as a guide, we select values of g = 48°C/km and d = 27.5 km. While other value combinations will change the specific shape of the 450°C isotherm, the changes are not sufficient to affect our interpretation (fig. S14). Setting the base of seismicity to the 450°C isotherm matches the top of seismicity to the 250°C isotherm for the entire range of g and d values, as shown by the suite of models in fig. S14. The depth of the 80°C isotherm also matches with the depth limit of shale diapirs imaged on the seismic reflection data (e.g., Fig. 1 and fig. S5), suggesting that a temperature limit for shale mobility also controls the structural style of the shallowest part of the fault.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank P. Tong for providing the 1D Lombok velocity model, R. Mallick for the useful discussion, and A. Loasby, H. Zeng, J. Yao, P. Maung Maung, Q. Shi, W. Fadil, and W. Chen for arrival time picking. We thank MultiClient Geophysical for providing seismic reflection data, BMKG for providing seismic data from regional stations, ESA for providing the Sentinel-1 SAR images, and JAXA for providing the ALOS-2 SAR images under the Research Project ER2A2N050. Teleseismic data were downloaded through IRIS DMC. Funding: This work comprises Earth Observatory of Singapore contribution no. 317. Support is from the National Research Foundation Singapore, the Singapore Ministry of Education under the Research Centres of Excellence initiative, and an NTU Presidential Fellowship awarded to K.L. (no. 04INS000845A620). K.B. is supported by the Earth Observatory of Singapore (no. 04MNS001950A620). T.W. is supported by the National Natural Science Foundation of China (no. 41974017) and the National Key Research and Development Program of China (no. 2019YFC1509204). A.D.N. is supported by “Hibah Penelitian Dasar Unggulan Perguruan Tinggi, Kemenristek/BRIN 2019-2021”and “Hibah Riset Institut Teknologi Bandung 2019-2020.” S.W.i acknowledges partial funding from the Indonesian Ministry of Research and Technology/National Agency for Research and Innovation, and Indonesian Ministry of Education and Culture under World Class University Program. Author contributions: K.L., M.M., S.We., and T.W. curated and analyzed the data. M.M., A.D.N., Z.Z., and S.Wi. deployed the temporary array. K.L., K.B., and S.We. developed the methodology and interpreted the results. K.L. and K.B. wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Data are available in the NTU Data Repository ( Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article