In this work, we perform a comprehensive investigation of the rupture process, aftershock sequence, and thermal environment for the 2024 Mw 7.4 Calama earthquake. We use a range of seismic data, including accelerograms and seismograms from a dense local network, along with teleseismic recordings (Fig. 1 and supplementary Fig. S1), to constrain the spatiotemporal evolution of the event. Additionally, we use coseismic displacement from continuous GNSS data (Fig. 1) to invert for the final static slip distribution. Our results reveal that the 2024 event ruptured multiple asperities along a steeply dipping fault plane, extending across an unexpectedly large depth range of over ~50 km. This implies that the rupture might have reached beyond the cold core of the slab. Comparison with geodynamic models indicates that the rupture indeed extended beyond the region where conditions are favorable for efficient serpentine dehydration. These observations suggest a possible rupture transition from dehydration embrittlement to shear thermal runaway for large intermediate-depth earthquakes.

Local strong motion data and teleseismic waveforms show a complex pattern of P and S waves indicative of rupture of different sub-sources for the 2024 Calama event (supplementary Fig. S2). This motivates a subevent inversion for the spatiotemporal rupture process of the M7.4 earthquake based on teleseismic P and SH waves from 54 global seismic stations and regional full waveforms from 47 strong motion stations (supplementary Fig. S1). Our method partitions complex rupture observations into a series of simpler point source subevents with their own timing, location, source duration, and moment tensors. A nonlinear inversion is implemented by Markov Chain Monte Carlo sampling, and we iteratively increase the number of subevents until the waveforms are fit sufficiently well (See the Methods section for more details). Locations and moment tensors of subevents are free parameters and provide information on fault geometries when those are not otherwise constrained.

We find that the 2024 Calama M 7.4 earthquake can be best modeled with five distinct subevents, where the best-fit subevent number is optimized based on the trade-off between waveform misfit and model complexity (Supplementary Fig. S3). In particular, a 5-subevent model provides noticeably better fits to complex waveform features compared to models with fewer subevents (Supplementary Fig. S3). The five subevents span over ~20 s, with an eastward progression over a horizontal distance of ~25 km (Fig. 2a). In contrast, the vertical extent of the rupture is much larger, with the subevents spanning a depth range of ~50 km (Fig. 2b). The rupture progression initiated with a M 6.6 subevent centered at 125 km depth (E1), close to our relocated hypocenter at 128 km. The rupture then triggered deeper subevents E2 (M 6.9) and E3 (M 7.0), which are close, slightly east of E1 but at a larger depth of 142 km. As the rupture advances, subevents E4 (M 7.0) and E5 (M 7.2) occurred even deeper, with centroids at 160 km and 174 km, respectively. Waveform fits (Fig. 2d, supplementary Figs. S4, S5) and model uncertainties (Fig. 2c, supplementary Fig. S6) indicate that the source parameters for each subevent are well resolved. The wide, asymmetric duration distributions for later subevent parameters (supplementary Fig. S6), as truncated by the limits of the search algorithm, indicate that their durations are not tightly constrained with long-period data. However, subevent timings and locations are well constrained. Given that the rupture may have initiated at depths shallower than the centroid of E1 and extended to depths larger than the centroid of E5, the whole depth range ruptured by the M 7.4 earthquake likely exceeded 50 km.

The compact horizontal yet substantial vertical extent of the subevents suggests that the rupture propagated along a steeply eastward-dipping fault plane. This geometry aligns well with the steep nodal plane of the subevent moment tensors (Fig. 2b). We also relocated the aftershocks for 2024 Calama event by manually identifying the P and S phases of seismograms recorded at local seismic stations for all aftershocks reported by the Chilean Seismic Network (CSN) using a 3D velocity model; results are shown in Figs. 1c, 2c, and supplementary data 1. Our detected and relocated aftershock seismicity indicates an eastward-dipping band spanning depths of 120-180 km, which closely matches the subevent hypocenters (Fig. 2b). The early, shallow subevents E1-E2 account for only 20% of the total moment release but are associated with high aftershock density. In contrast, the later, deeper subevents (E3-E5) contribute the bulk (80%) of the seismic moment, yet have lower aftershock density and more scattered aftershock distributions (see histogram in Fig. 2b). This suggests a more complete stress drop in the larger asperities ruptured during the later stages of the earthquake, given that the larger seismic moment over a compact rupture area implies a higher stress release. This interpretation is supported by the depletion of aftershocks for the later rupture stage, as a thermally driven comprehensive rupture would leave minimal residual stress to trigger aftershocks.

To further validate the rupture characteristics and depth extent for the 2024 Calama, we conducted directivity analysis on the broadband recordings. If there were no vertical rupture propagation, the apparent duration of both downgoing and upgoing body waves would be identical, i.e., no Doppler effect. However, for a unilaterally propagating rupture along dip, we expect a directivity effect: shorter and longer duration in the direction of rupture and opposite to it, respectively. To constrain source durations, we use two local earthquakes, the 2017/04/15 M 6.2 and 2022/12/10 M 5.6, to generate empirical Green's functions (EGFs) (Fig. 1a), and deconvolve the downgoing P and upgoing SH waves for the shared regional and global seismic stations for the apparent source time functions (STFs). Since the smaller events have similar depth and focal mechanisms, the derived EGFs allow removal of the shared path and site effects. Following deconvolution, we correct the stretching on the apparent source duration caused by the horizontal rupture toward the east as inferred from the subevent model. We then invert for the rupture dimension and velocity through regression, as those are linear functions of the vertical slowness and the corrected source durations for each station and phase.

The apparent STFs for downgoing P and upgoing SH waves have distinct durations and shapes (Fig. 3). Convolved with EGFs, they provide excellent fits to the observed waveforms of the Calama earthquake (Supplementary Figs. S7-S10). The downgoing STFs show coherent phases that align with the five subevents determined from the waveform inversion, with total duration of ~12-13 s (Fig. 3a). In contrast, the upgoing STFs have much longer durations of ~30 s, which suggests pronounced downward rupture directivity (Fig. 3b). The upgoing STFs also contain over 10 pulses, which reflect smaller scale asperities superimposed on the 5 major, lower frequency subevents ruptured coseismically (Fig. 3b). These details are preserved along the shorter travel distance and less-attenuating upgoing path, whereas the long teleseismic path for the downgoing waves acts as a low-pass filter. Applying a low-pass filter below 0.4 Hz to the upgoing STFs reveals comparable longer period subevent pattern (Supplementary Fig. S11). Meanwhile, the temporal gaps in the upgoing source time functions are necessitated by the waveform data in the deconvolution (Supplementary Fig. S12). Therefore, we interpret these high-frequency asperities as real, resolvable features of the rupture. Our regression analysis on the observed STF durations yields an optimal rupture length of 78 km and average velocity of 4.2 km/s, along a fault plane with 71° eastward dip. This angle coincides with the average dip of the 5 subevents' steeper fault planes. Moreover, the large downdip rupture dimension aligns with the depth range of the subevents and aftershocks (120-180 km) and explains the anomalously broad teleseismic sP and sS depth phases (supplementary Fig. S4). Our relocated aftershock distribution aligns well with the sub-vertical plane inferred from the subevents (Fig. 2). We also searched for precursor events and further aftershocks using template matching. No precursors were observed and the aftershock productivity is high compared to other similar events in the area (supplementary Figs. S13, S14).

The mainshock produced horizontal displacements of ~1 cm extending more than 250 km from the epicenter and a subsidence of ~20 cm above the rupture. Using these GNSS-derived displacements, we estimated the geodetic co-seismic slip distribution (Fig. 4). We tested both the sub-horizontal and sub-vertical planes based on the focal mechanism and our proposed hypocenter (Supplementary Fig. S15). Acknowledging the inherent challenge of resolving fault geometry and the detailed slip distribution for such a deep source, we found that both planes can fit the GNSS observations (supplementary Figs. S16-S18). However, the fit for the vertical plane model is slightly better and has relatively smoother inferred coseismic slip. The insignificant misfit differences suggest that geodetic data may not tightly constrain fault geometry. Nonetheless, the effective slip distribution for the sub-vertical plane is consistent with our subevent inversion, directivity analysis, and aftershock patterns. In conjunction, these provide much stronger seismological constraints on the rupture geometry.

The large depth extent of the 2024 Calama event renders it an important probe for understanding the mechanisms of intermediate-depth earthquakes in general. The event likely ruptured a substantial portion, if not reached the base of the subducting lithosphere. This is further supported by the distinct character of later waveform arrivals in the strong motion data, which necessitates deep depths of subevents E4-E5, as indicated by our forward modeling tests (Supplementary Fig. S19). The rupture could extend well beyond the inferred 600 °C to 650 °C isotherm at depths of 120-180 km, which bound the region capable of hosting dehydration embrittlement, a key mechanism invoked to explain the intermediate-depth seismicity in northern Chile. Compared with other intraplate earthquakes, 2024 Calama had a significant number of aftershocks, which suggests that at least part of its rupture occurred in a colder, brittle regime. This is more analogous to shallower events such as the 2007 Michilla M 6.7 and 2020 Calama M 6.9 earthquakes, rather than events in high-temperature zones (>600 °C) where aftershocks are typically suppressed.

We further assess the temperature environment of the region by computing steady-state, semi-kinematic thermal geodynamic models for the subduction zone along a 2D cross section. The slab geometry is constrained by multiple independent datasets, including regional seismicity, seismic tomography, and the latest high-resolution CSN earthquake catalog. Our models also incorporate the incoming plate thermal structure inferred from plate age, and an adiabatic temperature gradient (see Methods section). As for all such models, isotherms are strongly affected by plate age and convergence rate, the product of which is a scaled Peclet number, the "thermal parameter"; it determines how deep cold isotherms are subducted before they diffuse away. Our results indicate that although the M 7.4 earthquake likely nucleated near the cold core of the slab, its ~78 km downdip rupture extended far beyond the 650 °C isotherm that bounds the zone where serpentine can efficiently dehydrate to olivine (Fig. 5). This finding is robust to variations in thermal modeling setups and other methodological choices. Adopting alternative slab geometries (e.g., Slab2) or varying slab age still implies that 2024 Calama ruptured beyond the serpentinite dehydration isotherm (Supplementary Fig. S20). Even under the alternative hypothesis that the slab resides deeper, the thickness of its cold core is still well constrained. Meanwhile, our thermal model is consistent with other thermal models for the wider region; those also feature a cold core too thin to contain the entire Calama rupture (Supplementary Fig. S21). Note that other hydrous minerals such as chlorite, talc, and amphiboles can also dehydrate under varying pressure-temperature conditions. However, most amphiboles dehydrate at lower temperatures and thus release fluids at shallower depths, while chlorite and talc are present in substantially lower abundances than serpentine. As a result, antigorite serpentine remains the dominant contributor to fluid release and large-scale faulting. Therefore, once the rupture extended beyond into warmer regions, mechanisms other than dehydration embrittlement likely governed further rupture propagation.

We suggest that the rupture likely transitioned into thermal runaway, involving localized shear heating that weakens the rock and sustains rupture propagation. Thermal runaway alone is unlikely to be a general mechanism for earthquakes, due to the lack of a spontaneous self-localization. However, slow deformation over geological time scales allows significant strain to accumulate. Once a triggering mechanism like dehydration embrittlement can initiate localized failure, shear thermal runaway, even at high temperatures, can act as a positive feedback loop of shear heating that may substantially weaken the fault, generating large slip asperities, and as a consequence more complete moment release. Our observations agree with such a transition: the second rupture stage occurs at depths with sparse historical seismicity and low aftershock density (Fig. 1), yet releases 80% of the total seismic moment (E3-E5) (Fig. 2). This suggests that a secondary mechanism such as shear thermal runaway can play a critical role for fault weakening at larger depths, including for intermediate-depth earthquakes.