Deep intra-slab rupture and mechanism transition of the 2024 M_w 7.4 Calama earthquake

Domain: earth_science
Published: 2025-08-30
Editorial outcome: minor_revision

Abstract
Whilesubductionzonehazardisdominatedbythemegathrust,intermediate-depth(70–300km)earthquakeswithintheslabcanlikewisehavecatastrophicimpacts.Theirphysicsremainsenigmatic,withsuggestedmechanismsincludingdehydrationembrittlementandthermalrunaway.Here,weinvestigatethe2024Chile,M 7.4intermediate-depthearthquakeandcomparetheruptureextentwithtemperatureconditionsfromthermo-mechanicalmodels.Werecordregionalgeodeticco-seismicdeformationandhigh-resolutionseismicityassociatedwiththistypeofevent.Ouranalysesrevealacomplexrupturespanninganexceptionaldepthrange,withdistinctasperitiespropagatingdeepintothesubductinglithosphere.Comparisonwiththermalmodelsshowsthatwhiletheruptureinitiatedwithinthecoldslabcore,itextendedwellbeyondthe~650°Cisothermthattypicallydelineatestheboundaryforefficientserpentinedehydration.Wesuggestthattherupturelikelyinitiatedwithdehydrationembrittlementwithinthecoldcorebutthenpropagatedintothewarmerregionsthroughshearthermalrunaway.Thisimpliesatransitionofmechanismsthatfacilitateslarge-scaleruptureandactivatestypicallyaseismic,high-temperatureslabregions.Ourfindingshighlighttheimportanceofconsideringinteractionsbetweenrupturemechanismsaswellasslabthermalandcompositionalsettingstobetterunderstandtheprocessesgoverningintermediate-depthearthquakes.

Introduction
Most earthquakes occur along subduction zones at convergent plate boundaries, where tectonic plates collide and one plate sinks beneath another. Among these events, intermediate-depth earthquakes are defined as those occurring between 70 and 300km depth range1,2. Although intermediate-depth earthquakes generally pose lower seismic hazard than shallower events, these events can still produce destructive ground accelerations3,4, potentially resulting in substantial fatalities and extensive economic losses5. An example is the M 7.8 1939 Chillán intermediate-depth earthquake6, which tragically led to a larger loss of life in Chile than any of the megathrust events on record. An improved understanding of the rupture behavior and maximum magnitude of intermediate-depth earthquakes is thus imperative for seismic hazard assessment.

Formostsubductionzones,intermediate-depthseismicityforms upperandlower(double)seismicplanesinsidetheslab7–9.Toexplain thestressesintheseplanes,forcesincludingslabpull10andtheflexure ofthesubductedlithosphere11havebeenproposed,whereunbending ofslabscanleadtodowndipcompressionalandtensionalstresswithin the upper and lower plane, respectively. However, given that the pressureandtemperatureatintermediate-depthsexceedthosetypical of the brittle-ductile transition12, the physical mechanism for shear faulting remains unclear. There are two major hypotheses: the embrittlement of rock caused by slab dehydration, where water releasedfromthesubductingplateincreasesfluidpressureandcauses thesurroundingrocktofracture1,13,andthermalrunaway,aninstability whereshearheatingsoftensrock,creatingafeedbackloopoffocused strain that culminates in rapid slip14,15. Observational efforts to test these hypotheses have typically focused on seismicity distributions9,16,17,whereasstudiesexaminingrupturepropertieshave revealed characteristics including relatively high stress drops, fast rupture velocities, and a systematic decrease in aftershock productivity with increasing depth, which may be influenced by slab dehydration18–20.However,theruptureextentoflargeintermediate-depthearthquakes,theirlocaltemperatureandmineralcompositions, and their interaction with the double seismic zone, remain poorly constrained.Thislimitsourcapabilitytotestthephysicalboundaries of the proposed source mechanisms including dehydration embrittlementandthermalrunaway.

OnJuly19,2024,anearthquakeofM 7.4struckNorthernChile near Calama, named here the 2024 Calama earthquake. This event occurredwithintheoceaniclithospherewheretheNazcaplatesubducts beneath South America. In the region, other similarly large intraplate,intermediate-deptheventshavebeenreported,theM 8.0 Calama 1950 and M 7.8 Tarapacá 2005 events (Fig. 1)6. The 1950 earthquakeoccurredalmostatthesamelongitudebutafewtenskmto theSouthfromCalama2024(Fig.1),albeitwithlocationuncertaintyof 50-100km21. Its rupture was constrained by analog teleseismic data only, indicating a normal faulting focal mechanism with poor constraints on rupture dynamics21. The 2005 Tarapacá event occurred further north with a subhorizontal rupture19,22,23. Cab- rera et al. 20 noted that deeper events located at elevated temperatures present relatively fewer aftershocks but some other studies24,25did notfind this pattern.

Atasimilarlatitudeasthe2024Calamaearthquake,severalwell- recordedeventsofmagnitudesbelow7,includingthe2007M 6.7and 2020M 6.9earthquakes,occurredwithinthe300°–400°Cisotherms in the oceanic crust and upper lithospheric mantle (Fig. 1)26,27. This patternsuggestsaroleofslabthermalstructureinconfiningruptures. However,the2024Calamaevent,withitssignificantlylargermagnitude,likelyinvolvedamorespatiallyextensiveruptureandfaultgeometry than previous earthquakes. Hence, it provides a unique opportunitytoinvestigateintermediate-depthearthquakes,providing animportantprobeofthethermodynamicregime.Forexample,the thermalstateoftheslab,whichaffectsmineralphasetransitions,plays a critical role in determining whether dehydration embrittlement is feasible,thusinfluencingruptureprocesses.

Inthiswork,weperformacomprehensiveinvestigationofthe rupture process, aftershock sequence, and thermal environment forthe2024M 7.4Calamaearthquake.Weusearangeofseismic data,includingaccelerogramsandseismogramsfromadenselocal network,alongwithteleseismicrecordings(Fig.1andsupplementaryFig.S1),toconstrainthespatiotemporalevolutionoftheevent. Additionally, we use coseismic displacement from continuous GNSSdata(Fig.1)toinvertforthefinalstaticslipdistribution.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 ~50km. This implies that the rupture might have reached beyond the cold core of the slab. Comparison with geodynamic models indicates that the rupture indeedextendedbeyondtheregionwhereconditionsarefavorable forefficientserpentinedehydration.Theseobservationssuggesta possible rupture transition from dehydration embrittlement to shearthermalrunawayforlargeintermediate-depthearthquakes.

Results
Subeventcharacterizationoftheruptureprocess
Localstrongmotiondataandteleseismicwaveformsshowacomplex patternofPandSwavesindicativeofruptureofdifferentsub-sources forthe2024Calamaevent(supplementaryFig.S2).Thismotivatesa subevent inversion for the spatiotemporal rupture process of the M 7.4earthquakebasedonteleseismicPandSHwavesfrom54global seismicstationsandregionalfullwaveformsfrom47strongmotion stations (supplementary Fig. S1). Our method partitions complex ruptureobservationsintoaseriesofsimplerpointsourcesubevents with their own timing, location, source duration, and moment tensors28–30. A nonlinearinversion isimplemented by Markov Chain Monte Carlo sampling, and we iteratively increase the number of subeventsuntilthewaveformsarefitsufficientlywell(SeetheMethods sectionformoredetails).Locationsandmomenttensorsofsubevents arefreeparametersandprovideinformationonfaultgeometrieswhen thosearenototherwiseconstrained.

We find that the 2024 Calama M 7.4 earthquake can be best modeled with five distinct subevents, where the best-fit subevent numberisoptimizedbasedonthetrade-offbetweenwaveformmisfit and model complexity (Supplementary Fig. S3). In particular, a 5-subeventmodelprovidesnoticeablybetterfitstocomplexwaveform featurescomparedtomodelswithfewersubevents(Supplementary Fig. S3). The five subevents span over ~20s, with an eastward progressionoverahorizontaldistanceof~25km(Fig.2a).Incontrast,the vertical extent of the rupture is much larger, with the subevents spanningadepthrangeof~50km(Fig.2b).Theruptureprogression initiatedwithaM 6.6subeventcenteredat125kmdepth(E1),closeto ourrelocatedhypocenterat128km.Therupturethentriggereddee- persubeventsE2(M 6.9)andE3(M 7.0),whichareclose,slightlyeast ofE1butatalargerdepthof142km.Astheruptureadvances,sube- ventsE4(M 7.0)andE5(M 7.2)occurredevendeeper,withcentroids at 160km and 174km, respectively. Waveform fits (Fig. 2d, supplementaryFigs.S4,S5)andmodeluncertainties(Fig.2c,supplementary Fig.S6)indicatethatthesourceparametersforeachsubeventarewell resolved.Thewide,asymmetricdurationdistributionsforlatersube- ventparameters(supplementaryFig.S6),astruncatedbythelimitsof the search algorithm, indicate that their durations are not tightly constrained with long-period data. However, subevent timings and locationsarewellconstrained.Giventhattherupturemayhaveiniti- ated at depths shallower than the centroid of E1 and extended to depthslargerthanthecentroidofE5,thewholedepthrangeruptured bytheM 7.4earthquakelikelyexceeded50km.

The compact horizontal yet substantial vertical extent of the subevents suggests that the rupture propagated along a steeply eastward-dippingfaultplane.Thisgeometryalignswellwiththesteep nodalplaneofthesubeventmomenttensors(Fig.2b).Wealsorelo- catedtheaftershocksfor2024Calamaeventbymanuallyidentifying thePandSphasesofseismogramsrecordedatlocalseismicstations for all aftershocks reported by the Chilean Seismic Network (CSN) usinga3Dvelocitymodel31;resultsareshowninFigs.1c,2c,andsup- plementarydata1.Ourdetectedandrelocatedaftershockseismicity indicatesaneastward-dippingbandspanningdepthsof120-180km, whichcloselymatchesthesubeventhypocenters(Fig.2b).Theearly, shallow subevents E1-E2 account for only 20% of the total moment releasebutareassociatedwithhighaftershockdensity.Incontrast,the later,deepersubevents(E3-E5)contributethebulk(80%)oftheseis- micmoment,yethaveloweraftershockdensityandmorescattered aftershock distributions (see histogram in Fig. 2b). This suggests a morecompletestressdropinthelargerasperitiesrupturedduringthe laterstagesoftheearthquake,giventhatthelargerseismicmoment 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 wouldleaveminimalresidualstresstotriggeraftershocks.

Validationofthefaultgeometryandruptureextent
Tofurthervalidatetherupturecharacteristicsanddepthextentforthe 2024 Calama, we conducted directivity analysis on the broadband recordings.Iftherewerenoverticalrupturepropagation,theapparent duration of both downgoing and upgoing body waves would be identical, i.e., no Doppler effect. However, for a unilaterally propa- gatingrupturealongdip,weexpectadirectivityeffect:shorterand longer duration in the direction of rupture and opposite to it, respectively.Toconstrainsourcedurations,weusetwolocalearth- quakes, the 2017/04/15 M 6.2 and 2022/12/10 M 5.6, to generate empirical Green’s functions (EGFs)32,33 (Fig. 1a), and deconvolve the downgoingPandupgoingSHwavesforthesharedregionalandglobal seismic stations for the apparent source time functions (STFs)34,35. Sincethesmallereventshavesimilardepthandfocalmechanisms,the derived EGFs allow removal of the shared path and site effects.

Followingdeconvolution,wecorrectthestretchingontheapparent sourcedurationcausedbythehorizontalrupturetowardtheeastas inferred from the subevent model. We then invert for the rupture dimensionandvelocitythroughregression,asthosearelinearfunc- tionsoftheverticalslownessandthecorrectedsourcedurationsfor eachstationandphase36–38.

TheapparentSTFsfordowngoingPandupgoingSHwaveshave distinct durations and shapes (Fig. 3). Convolved with EGFs, they provide excellent fits to the observed waveforms of the Calama earthquake(SupplementaryFigs.S7–S10).ThedowngoingSTFsshow coherentphasesthatalignwiththefivesubeventsdeterminedfrom the waveform inversion, with total duration of ~12-13s (Fig. 3a). In contrast,theupgoingSTFshavemuchlongerdurationsof~30s,which suggests pronounced downward rupture directivity (Fig. 3b). The upgoingSTFsalsocontainover10pulses,whichreflectsmallerscale asperitiessuperimposedonthe5major,lowerfrequencysubevents rupturedcoseismically(Fig.3b).Thesedetailsarepreservedalongthe shortertraveldistanceandless-attenuatingupgoingpath,whereasthe longteleseismicpathforthedowngoingwavesactsasalow-passfilter. Applyingalow-passfilterbelow0.4HztotheupgoingSTFsreveals thecomparablelongerperiodsubeventpattern(SupplementaryFig.S11). Meanwhile,thetemporalgapsintheupgoingsourcetimefunctions arenecessitatedbythewaveformdatainthedeconvolution(Supple- mentaryFig.S12).Therefore,weinterpretthesehigh-frequencyaspe- ritiesasreal,resolvablefeaturesoftherupture.Ourregressionanalysis on the observed STF durations yields an optimal rupture length of 78kmandaveragevelocityof4.2km/s,alongafaultplanewith71° eastward dip. This angle coincides with the average dip of the 5 subevents’steeperfaultplanes.Moreover,thelargedowndiprupture dimension aligns with the depth range of the subevents and after- shocks(120–180km)andexplainstheanomalouslybroadteleseismic sPandsSdepthphases(supplementaryFig.S4).Ourrelocatedafter- shockdistributionalignswellwiththesub-verticalplaneinferredfrom the subevents (Fig. 2). We also searched for precursor events and further aftershocks using template matching. No precursors were observedandtheaftershockproductivityishighcomparedtoother similareventsinthearea20(supplementaryFigs.S13,S14).

The mainshock produced horizontal displacements of ~1cm extendingmorethan250kmfromtheepicenterandasubsidenceof ~20cmabovetherupture.UsingtheseGNSS-deriveddisplacements, we estimated the geodetic co-seismic slip distribution (Fig. 4). We testedboththesub-horizontalandsub-verticalplanesbasedonthe focal mechanism and our proposed hypocenter (Supplementary Fig. S15). Acknowledging the inherent challenge of resolving fault geometryandthedetailedslipdistributionforsuchadeepsource,we foundthatbothplanescanfittheGNSSobservations(supplementary Figs.S16–S18).However,thefitfortheverticalplanemodelisslightly betterandhasrelativelysmootherinferredcoseismicslip.Theinsig- nificantmisfitdifferencessuggestthatgeodeticdatamaynottightly constrainfaultgeometry.Nonetheless,theeffectiveslipdistribution for the sub-verticalplane is consistent with our subevent inversion, directivity analysis, and aftershock patterns. In conjunction, these provide much stronger seismological constraints on the rupture geometry.

Ruptureacrossthedehydrationisothermsuggestsamechanism transition
The large depth extent of the 2024 Calama event renders it an importantprobeforunderstandingthemechanismsofintermediate- depthearthquakesingeneral.Theeventlikelyrupturedasubstantial portion,ifnotreachedthebaseofthesubductinglithosphere.Thisis furthersupportedbythedistinctcharacteroflaterwaveformarrivals inthestrongmotiondata,whichnecessitatesdeepdepthsofsube- ventsE4-E5,asindicatedbyourforwardmodelingtests(Supplemen- tary Fig. S19). The rupture could extend well beyond the inferred 600°Cto650°Cisothermatdepthsof120-180km,whichboundthe region capable of hosting dehydration embrittlement13,39, a key mechanism invoked to explain the intermediate-depth seismicity in northernChile17.Comparedwithotherintraplateearthquakes20,2024 Calamahadasignificantnumberofaftershocks,whichsuggeststhatat leastpartof its ruptureoccurred ina colder, brittleregime. Thisis moreanalogoustoshallowereventssuchasthe2007MichillaM 6.7 and 2020CalamaM 6.9 earthquakes26,rather than events in high- temperature zones (>600°C) where aftershocks are typically suppressed.

Wefurtherassessthetemperatureenvironmentoftheregionby computingsteady-state,semi-kinematicthermalgeodynamicmodels forthesubductionzonealonga2Dcrosssection40.Theslabgeometry is constrained by multiple independent datasets, including regional seismicity, seismic tomography, and the latest high-resolution CSN earthquake catalog31,41. Our models also incorporate the incoming plate thermal structure inferred from plate age, and an adiabatic temperaturegradient(seeMethodssection).Asforallsuchmodels, isothermsarestronglyaffectedbyplateageandconvergencerate,the productofwhichisascaledPecletnumber,the“thermalparameter”;it determineshowdeepcoldisothermsaresubductedbeforetheydif- fuseaway.OurresultsindicatethatalthoughtheM 7.4earthquake likely nucleated near the cold core of the slab, its ~78km downdip rupture extended far beyond the 650°C isotherm that bounds the zonewhereserpentinecanefficientlydehydratetoolivine(Fig.5).This findingisrobusttovariationsinthermalmodelingsetupsandother methodological choices. Adopting alternative slab geometries (e.g., Slab241) or varying slab age still implies that 2024 Calama ruptured beyond the serpentinite dehydration isotherm (Supplementary Fig.S20).Evenunderthealternativehypothesisthattheslabresides deeper,thethicknessofitscoldcoreisstillwellconstrained.Mean- while,ourthermalmodelisconsistentwithotherthermalmodelsfor thewiderregion42,43;thosealsofeatureacoldcoretoothintocontain theentireCalamarupture(SupplementaryFig.S21).Notethatother hydrous minerals such as chlorite, talc, and amphiboles can also dehydrateundervaryingpressure-temperatureconditions44–46.How- ever, most amphiboles dehydrate at lower temperatures and thus releasefluidsatshallowerdepths,whilechloriteandtalcarepresentin substantiallylowerabundancesthanserpentine13,45.Asaresult,anti- gorite serpentine remains the dominantcontributor tofluid release andlarge-scalefaulting.Therefore,oncetheruptureextendedbeyond intowarmerregions,mechanismsotherthandehydrationembrittle- mentlikelygovernedfurtherrupturepropagation.

Wesuggestthattherupturelikelytransitionedintothermalrun- away, involving localized shear heating that weakens the rock and sustainsrupturepropagation.Thermalrunawayaloneisunlikelytobe ageneralmechanismforearthquakes,duetothelackofaspontaneous self-localization14,15.However,slowdeformationovergeologicaltime scales allows significant strain to accumulate. Once a triggering mechanismlikedehydrationembrittlementcaninitiatelocalizedfail- ure,shearthermalrunaway,evenathightemperatures,canactasa positivefeedbackloopofshearheatingthatmaysubstantiallyweaken thefault,generatinglargeslipasperities,andasaconsequencemore complete moment release47–49. Our observations agree with such a transition:thesecondrupturestageoccursatdepthswithsparsehis- toricalseismicityandlowaftershockdensity(Fig.1),yetreleases80% ofthetotalseismicmoment(E3-E5)(Fig.2).Thissuggeststhatasec- ondarymechanismsuchasshearthermalrunawaycanplayacritical roleforfaultweakeningatlargerdepths,includingforintermediate- depthearthquakes.

Methods
Earthquakedetectionusingtemplatematching
We used template matching70 to detect unreported earthquakes aroundthe2024Calamamainshock,bymeansoftheFastMatched Filteralgorithm71.Thiswasdonebyanalyzingcontinuousbroadband velocitywaveformsofsixstationsoftheCSNnetworksituatednearthe epicenterofthemainshockandatdistancesoflessthan3°.Following previous regional studies20,26, we used three component records downsampledto40Hzandbandpassfilteredfrom1Hzto20Hz.For theapplicationoftemplatematching,aninitialcatalogisrequired,for whichweusetheseismicityreportedbyCSN(www.sismologia.cl)in space-timewindows,aswillbeexplainedbelow.Foreacheventofthe initialcatalog,P-waveandS-wavearrivalsforallselectedstationswere calculatedusingalocal1-Dvelocitymodel72.Toestimatethetemplates aswellasthemoveoutsusedintheanalysis,weusetheverticalcom- ponenttoscantheP-waveandthehorizontalcomponentstoscanthe S-wave71. In this sense, we cut the continuous waveforms from 2s beforethetheoreticalarrivaltimesuntil3safter.Thus,wemadeone template waveform for each component of every station selected, whichrepresents5softherecordedevents.

Foreveryevent,wescannedthecontinuousdatabymovingthe templatesonesample,andpreservingtheirmoveouts.Ateveryposi- tion of the sliding windows, we calculated an averaged cross- correlation coefficient (CC) over all stations and components. In ordertoconsideronlywaveformswithanacceptablesignal-to-noise ratio(SNR),weexcludedstationswithanSNRlowerthan2.0inthe vertical component; others are weighted such that the sum of all weightsoverthestationsisunity.Lastly,weremovedthosetemplates withfewerthan4stationswithacceptableSNRontheverticalcom- ponent. This resulted in a time series representing the similarity between the templates and the continuous data. We used a daily detectionthresholdof15timesthemedianabsolutedeviation(MAD) of the daily time-series of averaged cross-correlationcoefficients to define the potential detection of an earthquake with significantly similarwaveformstothetemplates,assumingtheyoccuratthesame hypocentrallocationastheeventusedforthedetection.

Given the possibility of obtaining consecutive detections sepa- rated by a time of less than the template waveform’s duration, we removedtheseconsecutivedetections,leavingtheonewithhighest CCandremovingthesurroundingdetections.Themagnitudeswere estimated by computing the median amplitude ratio between the waveform of the templates and the detections over the stations, assumingthatatenfoldincreaseinamplitudecorrespondstooneunit increaseinmagnitude73.

Forouranalysis,weconsideredtwocaseswithdifferentspatial windows to define the initial catalogs before running template matching. The first is for comparative purposes and considers the sameradiusofrupture(21km)adoptedbyHerreraetal.26forthe2020 Mw 6.8 Calama intraslab earthquake, and the second considers a greater radius,of 47km,calculated using expressions proposed for anotherMw7.4event74,basedonthesourceradiusestimatedbyWells &CopperSmith75.Inthiscontext,wemaketwodifferentcatalogsas follows:
ForDS1,theinitialcatalogcontainseveryearthquakereportedby CSNthatoriginatedbetweenonemonthbeforeandonemonthafter the mainshock, and located at a distance less than 21km from the reportedhypocenterofthemainshock.Withthiscriterion,33initial eventswereconsidered,ofwhich32wereusedastemplates,and150 werefinallyreportedafterrunningtemplatematching.ForDS2,the initialcatalogiscomposedofeveryearthquakereportedbyCSNand located in a distance of less than 47km from the location of the mainshock,inatemporalwindowbetweenonemonthbeforeandone monthafterthemainshock.Inthiscase,84initialearthquakeswere considered, 82 were adopted as templates, and a total of 314 final eventswerereportedafterapplyingtemplatematching(supplemen- taryFigs.S13andS14).

Subeventinversionforruptureprocess
Weapplyamultiple-subeventinversionmethod29,30toconstrainthe sourceparametersofthesubeventsfor2024Calamainsteadofkine- maticfinite-faultinversion,asthesubeventparameterizationavoids potentialbiases from assuming a fixed fault geometry and allowsa morecomprehensiveexplorationoftheparameterspacewithfewer freeparameters.Subeventinversionwasalsoselectedoverback-pro- jection,asourprimarygoalistoresolvedown-diprupturedirectivity,a characteristic to which back-projection is largely insensitive. The number ofsubevents is determined by iteratively adding subevents untiltheimprovementindatafitbecomesinsignificant.Eachsubevent ischaracterizedby5nonlinearparameters:horizontallocation,depth, centroid time, and centroid duration, except for the firstsubevent, whosehorizontallocationisfixedattheCSNhypocentertoprevent simultaneousshiftingofallsubeventlocationsandseismogramtim- ings. Additionally, each subevent has 5 linear parameters corre- spondingtoitsdeviatoricmomenttensorcomponents.

Theinversionemploysatwo-stagealgorithmtoefficientlyexplore parameterspace.Intheouterstage,wesearchfornonlinearlocation and timing parameters using a Metropolis-Hastings Markov Chain Monte Carlo(MCMC)sampler.DuringeachMCMCrandomwalk,one nonlinear parameter is perturbed while others remain fixed, which allowsahighacceptancerateandsamplingefficiency76.Foreachsetof nonlinearparameters,welinearlyinverttheseismicdatatoestimate subevent moment tensors77. This separation of nonlinear and linear parameter estimation reduces computational cost and enhances efficiency.

Wegenerate192MCMCchainsandretainthebest-performing48 chainsto mitigate the effects ofbadinitialsamples and avoidlocal minima.ABayesianframeworkisusedintheinversiontoincorporate data errors and model priors, which allows estimation of model uncertainties. Uniform prior distributions are assigned to nonlinear parameters, with horizontal location priors informed by the spatial densityofearlyaftershocks.Althoughinstrumentalnoiseisminimal, weintroducea10%empiricaldataerrorrelativetomisfittoaccountfor inaccuracies in wave propagation assumptions. Model uncertainties arequantifiedasthe95%confidenceintervalsoftheposteriorprob- abilitydensityfunctions,definedbythe2.5and97.5percentilevalues oftheMarkovChainsamples,withoutassumingaGaussianfittothem.

We use P-wave displacement and velocity waveforms from 54 teleseismicstations,SH-wavedisplacementwaveformsfrom53tele- seismic stations, and three-component velocity waveforms from 47 regionalstronggroundmotionstations.Datawereobtainedfromthe GlobalSeismicNetwork(GSN),theInternationalFederationofDigital SeismographNetwork(FDSN),andtheChileanNationalSeismological Center(CSN).Weightsareadjustedtoachievecomparablecontribu- tionstothemisfit,withaweightingratioof3:1:0.04forteleseismicP andSHwaves,andregionaldata,respectively.Instrumentresponses andlineartrendswereremoved,andteleseismicdatawereband-pass filtered between 0.02–0.2Hz, while regional data were filtered between0.02–0.15Hz.Weallowtimeshiftsofupto2sforteleseismic P waves, 5s for teleseismic SH waves, and 2s for regional data to account for path effects and picking errors. Green’s functions were calculatedusingalocal1Dvelocitymodel72forthecrustandtheIAS- PEI91modelfordeeperEarthstructure.TeleseismicGreen’sfunctions were computed using a hybrid method combining the propagator matrixapproachandraytheory78,79,whileregionalGreen’sfunctions were computed using the frequency-wavenumber integration method80.TheoptimalsubeventmodelfortheCalamaeventisshown inTableS1.

Geodeticinversionoftheslipdistribution
We used GNSS observations from CSN (https://doi.org/10.7914/SN/ C1),togetherwithstationsfromtheInternationalGNSSService(IGS) network81andthoseinstalledbytheANILLO-PRECURSOR82andIPOC83 projects. All GNSS data were processed using the differentiation strategy84, using the Bernese 5.4 software85. The daily results were accumulated to generate time series that are compatible with the ITRF2020 reference frame86. The resulting coseismic displacements aresummarizedinTableS2.

Basedontheearthquakefocalmechanismandhypocenterloca- tion,wemodeledtwofaultplanes:(1)asub-horizontalplane(dip=21°, strike=172°)and(2)asub-verticalplane(dip=69°,strike=341°).We compute Green’s functions for both planes using triangular elastic dislocation modeling87. We then perform a generalized least square inversion to estimate the coseismic slip distribution, where prior informationisdefinedbytheEqualPosteriorInformationCondition (EPIC)88.EPICappliesaspatiallyvariablesmoothingconstraintbased onadiscreteLaplacian,whichallowsregionsofthefaultplanewith denseGNSSdatacoveragetoretainmoredetailedslipfeatureswhile enforcing strongersmoothinginpoorlyconstrained faultareas.We applied positivity constraints to the dip-slip components and con- strainedthemaximumvaluesofthesliponthestrikedirection.Both fault plane models achieved good fits to the GNSS data, while the vertical plane produced slightly better RMSE residuals (1.48mm vs. 1.51mm).Additionally,theverticalfaultplanealignsmorecloselywith thespatialdistributionofobservedaftershocks.

Sourcetimefunctiondeconvolutionforrupturedirectivity andextent
We designed a time-domain deconvolution approachto extract the sourcetimefunctionfromobservedseismicwaveforms.Theproblem minimizes the mismatch between the observed data u(t), and the convolution of the source function s(t) with a predefined empirical Green’sfunctionG(t),

min s 1/2 ||G(t)*s(t) - u(t)||2 + λ*L12(s) (1)

inwhichL12(s)isamixedL-L normregularizationtermtoprevent overfittingofthedata.Themixednormtermisdefinedas

L12(s) = sum g=1 to M ||sg||2 (2)

wheres isthesubvectorofs(t)correspondingtogroupsizeg.Here, wechoosegas20samples(2swindowforsamplingrateof0.1s),and we choose λ as 0.05 to balance sparsity and smoothness of the resolvedsourcetimefunctions.

The optimization problem is solved using the CVXPY library89, with the non-negativity constraints to reflect the physical non- negativityofthesourcetimefunctions.

Totesttherobustnessoftheobservedgapsinthesourcetime functions,weincorporateasecond-orderdifferenceregularizationto enforce smoothing in the source time function by penalizing rapid changesincurvature,

min s 1/2 ||G(t)*s(t) - u(t)||2 + λ*L12(s) + β*||D s||2 (3)

wherethefinite-differenceoperatorisrepresentedas

||D s||2 = sum i=1 to N-2 (s i+2 - 2s i+1 + s i)2 (4)

Theresultswithβas0,102(moderateforcedsmoothing),and104 (strong forced smoothing) are shown in supplementary Fig. S12. Retrievingacceptablefitstothedatarequiresgapsinthesourcetime functions.

Afterweobtainthesourcetimefunctions,weestimatethesource durationsforeachtrace.Thesedurationsaredeterminedbyidentify- ingthemain,coherentenergyreleasethatisclearlyvisibleabovenoise and maintains a consistent shape across stations. We confirm its robustnessbycomparingourresults(Fig.3)tothosederivedfroma standardquantitativecriterion(STFamplitude≥25%ofthepeak,used by the IRIS DMC source time function products90) (Supplementary Fig.S22).Whilethistestresultsinqualitativelysimilarrupturepara- meters,therigidthresholdintroducesunnaturaltruncations,particu- larlyforteleseismicP-waves.Therefore,weretainourinitialselections as they provide a more physically meaningful measure of source duration.Thesefinaldurationsarethenusedtocorrectthehorizontal rupture directivity constrained from the subevent locations (25km towardstheeast).Wefurtherexaminetherelationshipbetweensource durationsandverticalslownessbyfittingalinearregressionmodelto observational data to estimate the vertical rupture directivity. The slopeinterceptoftheregressionlineinFig.3cprovidesanestimateof the inverse of the vertical rupturevelocity, and the x-axis intercept indicates the source duration, which allows the characterization of rupturedimensionalongthedip.Toquantifyuncertaintiesoftheslope andtheintercept,weconsidertwotimesthestandarddeviationerrors oftheirrespectiveestimates;thoseprovideupperandlowerbounds fortheestimatedrupturevelocityanddimension,respectively.

Thermalsimulationsfortheslabtemperatures
WeconductthermalsimulationsforsubductingslabsatChile,Japan, and the Aleutian subduction zones. We compute two-dimensional thermal models initialized from a plate cooling model of reference plate thickness of 95km91 based on ages of oceanic plates92, which represents the thermal state of the oceanic lithosphere atincipient subduction. For the slab surface geometry, we use a slab interface modelestimatedusingtheCSNseismicitycatalogforthetwonorthern Chiletransects(Fig.5a,b),andtheslabsurfacefromSlab241forthe othertwotransects(Fig.5c,d).ForthenorthernChiletransects,we alsotesttheslabsurfacefromSlab2;thisdoesnotaffectourconclu- sions (Supplementary Fig. S21). We project convergence velocity92 ontothetransects(TableS3)becausetheplatemotionvelocityisnot aligned with the direction perpendicular to the trench. With these initial conditions for each subducting slab, we use a finite element modeltoestimatethethermalstructure40.Thesemodelsincorporatea kinematicslabwithprescribedgeometryaswellasadynamicmantle wedge with a composite rheology including both diffusion and dis- locationcreep.Theslabdecouplesfromtheoverridingplateat50-km depth. The mesh resolution is 2km, with refinement to 1km at the mantlewedgeandslabsurface.Themantlepotentialtemperatureis 1300°Candadiabatictemperaturegradientis0.3°C/km91,93,94.Tests assumingthattheslabdecouplesfromtheoverridingplateat80-km depth show no substantial difference from the thermal simulation resultsusedinthisstudy.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors present different lines of evidence that the 2024 Calama intermediate-depth earthquake ruptured vertically almost 50 km, suggesting that the rupture went beyond the 600-650C isotherm and thus suggesting dehydration embrittlement cannot explain the rupture extent. I think the observations do support this interpretation, but I think some of the more speculative interpretation are not supported and can be part of the discussion but downplayed a little.\nI have some minor comments on various parts of the paper and a major comment on the interpretation of rupture of asperities.\n1) Lines 59-61\nThe authors state that much of the work has been focused on earthquake locations (which is true) but does not include other aspects people have looked at like stress drops, directivity and rupture velocity or aftershock productivity (which is later discussed). I would add those topics and references here.\n2) Lines 73-75\nThe smaller events around Calama are not shown in the Figure. By the way, the notation of Supp. figures is not clear, maybe because the paper was prepared for another journal? In any case, figrue numbers are not correct or point towards a different figure.\n3) Lines 91-93\nThe authors say the thermal models show that part of the rupture extended beyond the region where serpentinite dehydration is possible. But there is recent work that talks about other minerals that may suffer dehydration as well in different conditions (see work by Ferrand for example). This is just to be a bit more complete in the discussion, not to say everything has to be dehydration.\n4) Line 98. One example of Figure S2 that does not correspond.\n5) Lines 120-121.\nAlthough it seems plausible and robust, I would still like to see how much each of the 5 sources are required by the data. I know the methods section probably will talk a little about this, how you decide 5 events and not 4, but is there a reasonable criteria to choose 4 instead of 5? How is this done? Can you show the reader the difference? AIn the methods section this is somehow explained but no comparison is shown. I would like to see it, I think it would be nice for the readers to see it. It is not a deal breaker.\n6) LINES 124-126\nThe authors go directly to the aftershocks without any introduction as to how they are detected, relocated, etc. I short sentence explaining the basics would help. I see later in the text this is expanded and the methods section does it too, but it seems disconnected at this point in the text.\n7) Line 131-133.\nCan you expand a little as to why the observations discussed before support a complete stress drop of larger asperities? Because the moment is larger? Similar to the thermal suppresion of aftershocks.\n8) Lines 152-153.\nNot sure if I am reading correctly. But the figures showing the waveforms and comparison with the inversion are for the source inversion with 5 subevents, not the EGF you used for looking at directivity.\n9) Comparison of upgoing and downgoing waves (Fig 3).\nThis figure is trying to show that the STF duration is longer for upgoing waves than for downgoing waves, so a Doppler effect. But the authors suggest that the P waves clearly show the 5 subevents (ok, one could argue that is true), but for the SH waves I don't see that. And to me, that is confusing. Why would we have 5 subevents for downgoing waves and not the 5 subevents for the upgoing waves? The argument logic is not clear to me.\n10) Lines 165 and Figure.\nCan you show the expected pP, sP arrival times?\n11) Lines 263-265.\nAlthough the inference is interesting, it is highly speculative. The STF, to me, does not show a rupture by jumps between larger and smaller asperities. I am sorry but I do not agree that is shown in Figure 3. And I think this is where the authors interpret the data much beyond what it is showing. Maybe this is my bias, but I don't think this is supported.\n12) Lines 270-271.\nSimilar to the previous point, I think, you would need more evidence to support this. It might be interesting, but the data does not support that.", "reviewer_2": "The physics of intermediate depth earthquakes within the slab remains engimatic. In this work the authors investigate the rupture of the 2024 Chile Mw7.4 intermediate depth Calama earthquake, and compare the rupture extent with temperature conditions from thermo-mechanical models, and find that the Calama earthquake likely initated by dehydration embrittlement within the cold slab core, then propagated deeper depth, and transitioned from dehydration embrittlement to thermal rauway. In overall, this work is quite interesting and of substantial scientific importance. In order to further improve the impaction of this work, a minor to moderate revision is needed before being accepted. The major scientific issues needed to be further clarified are listed as following:\n1. In this work, the spatiotempral rupture process of the 2024 Calama earthquake is investigated with a subevent inversion instead of a Finite Fault Inversion (FFI). Please give a explanation for chosing subevent inversion instead of FFI, since most rupture process model are built with FFI method.\n2. As for determining rupture speed of the rupture, the authors conducted directivity analysis instead of Back Projection method, please give a explanation for this.\n3. From Fig. 3a, b, there are STF waveforms that last longer after the STFs labeled in red, which are followed by amplitudes about the same as the ones labeled in red, how were these red-labeled waveforms selected? This is very important and determines the duration of the rupture.\n4. In the directivity analysis, two local earthquakes, the 2017/04/15 Mw6.2 and 2022/12/10 Mw 5.6 are used as the EGFs, however, in the figures, the authors do not show the locations of these two earthquakes, please add them to the figures, such as in Figure 1. or Figure 2.\n5. The most important new finding in this paper is that the 2024 Calama earthquake ruptures over a small scale in horizontal direction while over a wide range in vertical direction, passing across the 650-degree isotherm. Thus the authors speculate that subevents at shallow depths (E1 to E3) may be due to the embrittlement of rock caused by slab dehydration, while those at deeper depths are due to a thermal runaway, which is of great significance. In overall, there are three main evidences to support this hypothesis: first, the STF functions exhibit an obvious time gap between the subevents; second, the distribution pattern of the aftershocks is more consistent with the rupture location of the sub-events, and the third, the earthquake was also found to have crossed the 650-degree isotherm through temperature simulations. However, these evidences are relatively circumstantial evidences. Can we compare the characteristics of seismic waveforms and their spectrums of the earthquakes occurring above the 650-degree isotherm due to dehydration embrittlement, as well as those of the earthquakes occurring at greater depths due to thermal runaway, with the theoretical seismic waveforms and spectral characteristics inspired by rupture of the shallow subevents of the present earthquake (E1~E3), as well as of the deeper subevents (E4~E5) and thus give possible direct evidences?\n6. In Figure 5. I can not find thin red lines (represent the 1000oC isotherms) in this figure. In addition, the background gray circles are too small and blur to be observed, please redraw this figure.\n7. From Extended Data Figure 4, the marginal probability density distributions of subevent rupture durations are wide and show a clear asymmetry, and it looks very strange to see these probability density distributions cut off abruptly in EP3 through EP5, so please check the results carefully to see if they make sense.\n8. Line 177~178, the authors mentioned that they tested both the sub-horizontal and sub-vertical planes based on the focal mechanism, I would like to see the geometries of the sub-horizontal and sub-vertical planes in three dimensional space as well as the hypocenter location of the mainshock.\n9. As mentioned in Line 256~265, the rupture speed of the 2024 Calama earthquake exceeds ~93% of the local shear wave velocity, which is speed falls in the forbidden zone between Rayleigh and shear wave speeds, so that the authors proposed that the 2024 Calama earthquake should rupture at both regular subshear and episodic supershear speeds. If so, maybe we can observe some kind of Doppler effects through seismic waveforms. If possible, the authors should try to find some evidences to verify the existence of the supershear rupture of this earthquake, such as using the method proposed by Vallée & Dunham (2012).\n10. Figure 4. The figure caption of Fig. 4a is wrong, for Fig. 4a, the blue and green arrows are observed and predicted horizontal component displaces instead of vertical displacements, please fix this error.\nSuggested reference:\nVallée, M. & Dunham, E.M. (2012). Observation of far-field Mach waves generated by the 2001 Kokoxili supershear earthquake. Geophys. Res. Lett., 39,. L05311.", "reviewer_3": "Jia et al. investigate the rupture characteristics of the large, intermediate-depth Calama earthquake, finding that it had a steeply dipping fault plane extending ~75 km down from the hypocenter. By themselves, the combined analyses are rigorous and convincing, making this event one of the best constrained intermediate depth earthquakes of this magnitude I have read about. But what makes the earthquake potentially compelling for publication in Nature Communications is its context within the thermal conditions of the slab. As the authors explain well, the rupture propagated from the cold core of the slab, where earthquakes are known to occur, out into the hotter regions of the slab where brittle, elastic-rebound processes are unexpected. In my opinion, the manuscript is convincing and in good shape, other than some minor comments, provided below.\nWhat I find myself wondering about is the thermal modeling and relationship between the Calama earthquake and temperature. [Full disclosure, I am not an expert in intermediate depth earthquakes or in the thermal state of the Nazca slab.] Fig. 5 of the manuscript shows that this section of the Nazca slab has substantially more earthquakes along the deeper cold-to-hot transition in the subducting slab than the other three settings shown. I have not looked closely at other subduction zones to see if this uniqueness extends globally, but I think this is something to probe further about the setting of the Calama earthquake since it is essential to the argument about rupture process. Could it be that the slab in the Calama region is actually slightly deeper than currently proposed, so these events (and the 2024 rupture) are still mostly within the cold slab core? Or do the authors know of alternative interpretations for the number of deeper intermediate earthquakes along the base of the slab here? Addressing this would go a long way towards solidifying the interpretive aspects of this study.\nMinor Comments:\nAn aspect of this problem the authors allude to but never address directly is the topic of strain rates. How it relates to the brittle-ductile transition is implicit in the rupture mechanisms for intermediate-depth earthquakes, but I would like to see it laid out more explicitly (cf., ll. 55-56). One thing in particular that strikes me is that at typical strain rates and temperatures of ~1000ºC, even with a mechanism like thermal runaway to localize strain it is hard to imagine having enough elastic strain energy stored for an earthquake. I understand that this is one of the key points in the manuscript. I think it would be helpful to lay it out for the reader.\nI am aware of several other thermal models for subducting slabs. Qualitatively, the results here look similar, but obviously the details matter for interpreting this event. It would be good to have a comparison (even a brief one) and some references to those other models (e.g., somewhere in ll. 198-205)\nDetailed Comments:\nll. 26-34: I think it would be better to disentangle observations and analyses from interpretations about rupture mechanisms here, focusing first on the detailed aspects of the event, then the comparison to thermal models, and finally the connection to dehydration embrittlement or shear thermal instability\nl. 42: What is the point the authors are trying to make in saying that “the number of earthquakes drops markedly from shallow to deep”? Could this statement be removed?\nll. 57-58: The references for these intermediate earthquake processes are helpful, but it might also be useful to give a brief (one sentence) description of each process too\nl. 60: What exactly do the authors mean by “environmental conditions”? Temperature? Other properties?\nll. 59-62: How would knowing things like rupture extent contribute to constraining the process? Making this explicit will help guide the reader through the rest of the manuscript.\nll. 67-68: For an M 7.8, a few tens of km is not that far away. What is the location quality of the 1950 event?\nl. 68: I think there is no need for “However” at the beginning of this sentence.\nl. 69: I am unsure about what “constrained only by analog teleseismic data” implies. That it only has a focal mechanism or moment tensor?\nl. 69: I think it would be good to put a period after “indicating a normal faulting earthquake” and then talk about the 2005 earthquake separately. I initially thought the points were connected, but they are not. Also, just a point to note: the concept of “normal faulting” is a little different for intermediate depth earthquakes, since the slab is dipping and surrounded by rocks (i.e., there is not a requirement for a principal stress axis to be perpendicular to a free surface).\nll. 85-88: I would divide this sentence into two, one describing the information obtained from seismic data and the other describing the information obtained from co-seismic displacements\nl. 90: Maybe I missed it, but I do not think the “core of the slab” has been defined yet\nll. 96ff: I find the description of these results really clear and easy to follow\nll. 135ff: I think some references might be useful here (e.g., for EGF approaches, the stretching of the source duration caused by horizontal rupture, or the linearity of the rupture dimension/velocity)\nll. 153-154: I confused “coherent phases” with seismic phases at first. Is there a better word here?\nl. 157: Some numbers would be helpful here: How many pulses? What size asperities?\nll. 174ff: I would make the geodetic analysis its own paragraph. Also, just a comment that at this depth I suspect that it is challenging to resolve a slip distribution, much like it is challenging to resolve fault geometry.\nl. 189: What is a “complete half-section”?\nll. 193-197: I am unsure of the argument being made about aftershocks and temperature here. Sorry my comment is not more helpful.\nll. 236-237: If this were the case, we might expect to see a temporal relationship between the megathrust earthquake cycle and intermediate depth earthquakes. I know this has been argued for the region further north (where the 2005 Tarapaca occurred), but I am unaware of a similar argument here. Also, in very large megathrust earthquakes, these stresses are recovered, so it might not be a long-term effect. I guess my point is that I am unsure about whether high coupling is a “key factor” or not.\nll. 244ff: It seems that there is a little jumping back and forth between intermediate and very deep earthquakes. The authors should be very explicit throughout this paragraph.\nll. 349-350: What is the specific criterion for no longer adding subevents? Is it based on misfit alone? Or is there something like the Akaike Information Criterion applied to select the best number of subevents?\nl. 371: What is the “width” of the probability density function defined by? Is a Gaussian fit to it?\nll. 407-409: As I said above, I suspect that the fault slip determined by static co-seismic displacements is poorly constrained\nFig. 1: In panel 1a, “Magnitude” is spelled wrong.\nFig. 5: The 1000ºC isotherm is very hard to see. Could the authors make it slightly bolder?\nFig. 6: Just a note: the scale of the overriding plate makes it seem as if everything is occurring at fairly shallow depths"}, "author_response": "Response: We sincerely appreciate the constructive feedback from the reviewers and editorial team, the response to which has, in our opinion, significantly improved the quality of our manuscript and strengthened our conclusions. In our revision, we have addressed all comments by conducting new analyses to validate our conclusions, including waveform analyses to confirm the deep rupture extent, validation of subevent numbers, tests on different criterions on deconvolved source time functions, and a comparison of our thermal model against other established models. We have also expanded our discussion on key topics such as the physical mechanisms permitting high-temperature rupture and the rationale for our choice of inversion methods. Furthermore, we have refined the text and figures throughout the manuscript for improved clarity and precision. Below, we provide detailed responses to the reviewers’ comments. Reviewer comments are shown in black, and our responses in blue. Line numbers refer to the revised manuscript with tracked changes.\nReviewer: 1\nThe authors present different lines of evidence that the 2024 Calama intermediate-depth earthquake ruptured vertically almost 50 km, suggesting that the rupture went beyond the 600-650C isotherm and thus suggesting dehydration embrittlement cannot explain the rupture extent. I think the observations do support this interpretation, but I think some of the more speculative interpretation are not supported and can be part of the discussion but downplayed a little.\nI have some minor comments on various parts of the paper and a major comment on the interpretation of rupture of asperities.\nResponse: Thank you for the positive evaluation. We appreciate the constructive comments and our responses are detailed below item by item. Your comments are shown in black, and our responses in blue. Line numbers refer to the revised manuscript with tracked changes.\n1) Lines 59-61\nThe authors state that much of the work has been focused on earthquake locations (which is true) but does not include other aspects people have looked at like stress drops, directivity and rupture velocity or aftershock productivity (which is later discussed). I would add those topics and references here.\nResponse: We have added a brief introduction on these related topics and references in L69-74.\n2) Lines 73-75\nThe smaller events around Calama are not shown in the Figure. By the way, the notation of Supp. figures is not clear, maybe because the paper was prepared for another journal? In any case, figrue numbers are not correct or point towards a different figure.\nResponse: Thank you for pointing this out. The 2007 magnitude 6.7 event was mistakenly labeled as 7.7 in Figure 1, and we have corrected this error. Regarding the figure notation, the confusion is likely due to the formatting requirements specific to Nature Communications, which distinguishes between Extended Data figures and Supplementary figures, both different from the main text figures. We have now carefully reviewed all figure references and ensured they correctly correspond to the appropriate figures.\n3) Lines 91-93\nThe authors say the thermal models show that part of the rupture extended beyond the region where serpentinite dehydration is possible. But there is recent work that talks about other minerals that may suffer dehydration as well in different conditions (see work by Ferrand for example). This is just to be a bit more complete in the discussion, not to say everything has to be dehydration.\nResponse: Thank you for your suggestion. We agree that other hydrous minerals such as chlorite, talc, and amphiboles can also dehydrate under varying pressure-temperature conditions. However, most amphiboles dehydrate at lower temperature and thus release fluids at shallower depths, while chlorite and talc are present in much lower abundances. As a result, antigorite serpentine remains the dominant contributor to fluid release and large-scale faulting.\nWe added additional references (e.g., Schmidt et al., 1998; Ferrand et al., 2017; Ferrand et al., 2021) and broadened our discussions in L270-275 to include other potential dehydration reactions.\n4) Line 98. One example of Figure S2 that does not correspond.\nResponse: Thank you for noting this. The previous Figure S1 showed waveforms after applying the absolute value, which may have caused confusion. To improve clarity and better illustrate the presence of the five distinct pulses, we have updated Figure S1 to display the stacked displacement waveforms without taking the absolute value.\nRevised Figure S1. Complex pulses in the broadband strong motion waveforms in displacement.\nThe waveforms in the left panel are sorted by azimuth, while those in the right panel are organized based on the directivity parameter. The red line in the right panel represents the stacked strong motion waveforms.\n5) Lines 120-121.\nAlthough it seems plausible and robust, I would still like to see how much each of the 5 sources are required by the data. I know the methods section probably will talk a little about this, how you decide 5 events and not 4, but is there a reasonable criteria to choose 4 instead of 5? How is this done? Can you show the reader the difference? In the methods section this is somehow explained but no comparison is shown. I would like to see it, I think it would be nice for the readers to see it. It is not a deal breaker.\nResponse: To address your point, we have added a comparison between models with different numbers of subevents. The new Figure S2 presents the relationship between model misfit and the number of subevents. While the misfit generally decreases with more subevents (as expected), the rate of reduction slows significantly after five. The 5-subevent model is therefore identified as the preferred solution, as it provides a good fit to the data without adding unnecessary complexity.\nFigure S2 also compares waveform fits for the 4- and 5-subevent models. The 5-subevent model provides noticeably better fits to the complex waveform features, particularly on the East and North components at stations such as C1.AF01, C1.A22F, and C1.A18F. These improvements support our selection of the 5-subevent model as the most appropriate representation of the data.\nWe have included this comparison as a new supplementary Figure S2.\nNew Figure S2. Justification for the choice of the 5-subevent model. (a) Variation of the waveform misfit as a function of the number of subevents. (b-c) Waveform fits for the 4-subevent and 5-subevent models, respectively. Observed data are shown in black and synthetic seismograms are in red.\n6) LINES 124-126\nThe authors go directly to the aftershocks without any introduction as to how they are detected, relocated, etc. I short sentence explaining the basics would help. I see later in the text this is expanded and the methods section does it too, but it seems disconnected at this point in the text.\nResponse: We agree that introducing the aftershock detection and relocation methods earlier in the text improves clarity and helps guide the reader. In the revised manuscript L153-158, we have moved the following sentence to precede the comparison between subevents and aftershocks:\n“We also relocated the aftershocks of the 2024 Calama event by manually identifying the P and S phases at local seismic stations for all events reported by the Chilean Seismic Network (CSN), using a 3D velocity model (Potin et al., 2024); results are shown in Figs. 1c, 2c, and Table S1.”\n7) Line 131-133.\nCan you expand a little as to why the observations discussed before support a complete stress drop of larger asperities? Because the moment is larger? Similar to the thermal suppresion of aftershocks.\nResponse: We have expanded the discussion to clarify the link between the observations and the inference of a complete stress drop in the larger asperities.\nSubevents with larger seismic moments but occurring over similar spatial extents suggest higher stress drops, as more energy is released over a comparable area. This is supported by the relative depletion of aftershocks in these regions, which implies that most of the accumulated stress was released during the main rupture.\nFrom a physical perspective, the rapid and localized nature of thermal runaway processes facilitate large stress release, leaving minimal residual stress to drive aftershock activity. Additionally, the high temperature and pressure conditions at these depths, along with ductile strain localization, can inhibit the formation of new fractures, further suppressing aftershock generation.\nWe have incorporated this expanded explanation into the revised manuscript L165-170 to improve clarity.\n8) Lines 152-153.\nNot sure if I am reading correctly. But the figures showing the waveforms and comparison with the inversion are for the source inversion with 5 subevents, not the EGF you used for looking at directivity.\nResponse: We believe the confusion may be related to Extended Data Figures 2-3, which indeed show waveform fits from the finite-source inversion using five subevents. The EGF-based directivity results are presented separately in Supplementary Figures 2-5. We have clarified this distinction in the text to avoid further confusion.\n9) Comparison of upgoing and downgoing waves (Fig 3).\nThis figure is trying to show that the STF duration is longer for upgoing waves than for downgoing waves, so a Doppler effect. But the authors suggest that the P waves clearly show the 5 subevents (ok, one could argue that is true), but for the SH waves I don't see that. And to me, that is confusing. Why would we have 5 subevents for downgoing waves and not the 5 subevents for the upgoing waves? The argument logic is not clear to me.\nResponse: We agree that the five subevents are more clearly delineated in the downgoing P wave ASTFs than in the raw upgoing SH wave ASTFs. This apparent discrepancy does not imply a different number of subevents for up- and downgoing waves, but rather reflects the significant difference in attenuation along their respective travel paths. The teleseismic P waves travel long distances through the Earth's mantle, which strongly attenuates high-frequency energy and acts as a natural low-pass filter. Consequently, the resulting P-wave ASTFs are inherently smooth and reveal the macroscopic, low-frequency features of the rupture process, in this case, the five primary subevents. In contrast, the upgoing SH-waves travel much shorter regional paths and experience much less attenuation. Hence, their STFs preserve more high-frequency content which reflects smaller-scale complexities and heterogeneities within the main subevents. This rich high-frequency signal tends to obscure the underlying, larger-scale rupture patterns.\nTo demonstrate that the fundamental 5-subevent pattern is present in both datasets, we have applied a low-pass filter at 0.4 Hz to the SH-wave moment rate functions. This filtering process mimics the attenuation experienced by the teleseismic P-waves, allowing for a more direct comparison of the underlying rupture pattern. As shown in the new supplementary Figure S7, low-pass filtering the SH-wave ASTFs reveals the same five macroscopic subevents seen in the P-wave ASTFs. This demonstrates that the underlying rupture pattern is consistent across both datasets.\nWe have added this clarification to L195-204 to strengthen our measurement of the STF durations.\nNew Figure S7. Comparison of upgoing SH apparent source time functions at different frequency bands. (a) Original regional SH-wave ASTFs. (b) Regional SH-wave ASTFs after applying a 0.4 Hz low-pass filter. After filtering, the five dominant subevents become more clearly visible and better agree with the pattern observed in the P-wave ASTFs.\n10) Lines 165 and Figure.\nCan you show the expected pP, sP arrival times?\nResponse: We have added the timing of pP, sP, and sS arrivals to the waveforms in the Extended Data Figure 2, and revised the text in L208 to clarify that that sP and sS durations are longer than the main P and S wave durations.\n11) Lines 263-265.\nAlthough the inference is interesting, it is highly speculative. The STF, to me, does not show a rupture by jumps between larger and smaller asperities. I am sorry but I do not agree that is shown in Figure 3. And I think this is where the authors interpret the data much beyond what it is showing. Maybe this is my bias, but I don't think this is supported.\nResponse: We agree that our original wording may have overstated the interpretation. To clarify, our intent was not to suggest discrete \"jumps\" between larger and smaller asperities, but rather that the rupture appears to propagate through a sequence of asperities with varying sizes and stress drops. This interpretation is based on the multiple pulses observed in the source time function (STF) and the spatial separation of major subevents.\nWe have revised the text in L348-350 to present this idea more cautiously.\n12) Lines 270-271.\nSimilar to the previous point, I think, you would need more evidence to support this. It might be interesting, but the data does not support that.\nResponse: Similar to 11), we better clarified the evidence and presented the interpretation more cautiously in L355-358.\nReviewer: 2\nThe physics of intermediate depth earthquakes within the slab remains engimatic. In this work the authors investigate the rupture of the 2024 Chile Mw7.4 intermediate depth Calama earthquake, and compare the rupture extent with temperature conditions from thermo-mechanical models, and find that the Calama earthquake likely initated by dehydration embrittlement within the cold slab core, then propagated deeper depth, and transitioned from dehydration embrittlement to thermal rauway. In overall, this work is quite interesting and of substantial scientific importance. In order to further improve the impaction of this work, a minor to moderate revision is needed before being accepted. The major scientific issues needed to be further clarified are listed as following:\nResponse: Thank you for the positive evaluation. We appreciate the constructive comments and our responses are detailed below item by item. Your comments are shown in black, and our responses in blue. Line numbers refer to the revised manuscript with tracked changes.\n1. In this work, the spatiotempral rupture process of the 2024 Calama earthquake is investigated with a subevent inversion instead of a Finite Fault Inversion (FFI). Please give a explanation for chosing subevent inversion instead of FFI, since most rupture process model are built with FFI method.\nResponse: We agree that Finite Fault Inversion (FFI) is a powerful and widely used technique for imaging detailed rupture processes. However, for our investigation of the 2024 Calama earthquake, we choose subevent modeling approach due to its advantages in cases with notable unknowns, which we believe leads to a more robust and reliable primary model of the source. Our rationale is based on the following points:\nFirst, fewer assumptions on fault geometry. FFI requires defining a specific fault plane with fixed strike and dip, which can bias results if the true geometry is uncertain, as is often the case early in an event like the Calama earthquake. The subevent model avoids this by not requiring a predefined fault, which allows us to image rupture without imposing a possibly incorrect fault geometry.\nSecond, objective imaging of the rupture evolution. The subevent approach lets the data drive the reconstruction of rupture progression without enforcing assumptions about rupture velocity or direction. On the other hand, FFI is more constrained by model assumptions on rupture spatiotemporal evolution.\nFinally, improved inversion robustness. FFI involves hundreds to thousands of unknown parameters (e.g., slip, rake, rise time on subfaults) in the inversion, which can make the inversion unstable or non-unique, especially with limited seismic data coverage. On the other hand, the subevent model has much fewer nonlinear parameters (25 for the 5-subevent inversion), making it more robust and less prone to overfitting.\nThus, for the 2024 Calama earthquake, we prioritized reliability and interpretability by using the subevent method to reduce assumptions and generate a robust first-order view of rupture evolution.\nWe briefly addressed this point in L438-441.\n2. As for determining rupture speed of the rupture, the authors conducted directivity analysis instead of Back Projection method, please give a explanation for this.\nResponse: Our primary objective here was to constrain the down-dip rupture characteristics, including directivity and rupture speed. Resolving vertical directivity fundamentally relies on observing the Doppler effect from a wide geometric range of stations, and specifically by contrasting the apparent source durations of downgoing waves (recorded at teleseismic distances) and upgoing waves (recorded at local or regional distances). For this reason, we chose the subevent modeling which provides the depth-wise flexibility and resolution of the rupture process, as well as the directivity analysis based on the deconvolution and comparison of downgoing and upgoing apparent source time functions (aSTFs). These methods allowed us to quantify the Doppler-like shortening and lengthening caused by rupture propagation. The systematic variation in STF duration as a function of station azimuth and takeoff angle was further modeled to provide a robust estimate of the rupture speed in the down-dip direction.\nThe Back Projection (BP) method, which typically uses only high-frequency P-waves from teleseismic arrays, is an excellent tool for tracking the horizontal migration of high-frequency energy release. However, its reliance on a narrow cone of downgoing teleseismic waves makes it geometrically insensitive to rupture propagation along the dip of the fault. It cannot effectively distinguish between a shallow and a deep source of radiation within its time window, and hence it cannot resolve vertical directivity or the down-dip rupture speed.\nWe briefly addressed this point in L441-443.\n3. From Fig. 3a, b, there are STF waveforms that last longer after the STFs labeled in red, which are followed by amplitudes about the same as the ones labeled in red, how were these red-labeled waveforms selected? This is very important and determines the duration of the rupture.\nResponse: We aim to isolate the main coherent energy release from the rupture from later coda signals. To achieve this, our determination of durations focus on identifying STF segments that are not only clearly distinct from pre-event noise but also show clear visibility and generally consistent shape across a wide range of stations.\nTo test the robustness of our result, we repeat the analysis using a fixed, quantitative criterion: we systematically define the STF duration window with amplitude ≥ 25% of the peak, standard used by the IRIS DMC for their public source-time function products. This allowed comparison with an established, representative standard.\nThe overall source durations and rupture parameters derived from these automatically truncated STFs were qualitatively similar to those from our preferred selections. However, we found that the rigid 25% threshold introduced significant unnatural truncations and artifacts, particularly for the teleseismic P-wave STFs where complex wave propagation can distort signal shape. This confirms our results’ robustness while highlighting that our selection offers a more physically meaningful measure for the source durations.\nWe've added a new supplementary Figure S13 and clarified this point in L531-539 for transparency.\nNew Figure S13. Robustness test of the directivity analysis using a fixed duration criterion. This model is derived using a fixed 25% peak amplitude threshold to define source durations. The resulting apparent source time functions (left and middle panels) exhibit significant unnatural truncations not present in our main model (Fig. 3). However, the plot of vertical slowness versus duration (right panel) yields a rupture velocity and length consistent with our preferred results. This test confirms that our directivity analysis is robust while highlighting that our manual duration selection provides a more physically meaningful result by avoiding artifacts from an automated threshold.\n4. In the directivity analysis, two local earthquakes, the 2017/04/15 Mw6.2 and 2022/12/10 Mw 5.6 are used as the EGFs, however, in the figures, the authors do not show the locations of these two earthquakes, please add them to the figures, such as in Figure 1. or Figure 2.\nResponse: We have added these two EGF events in Figure 1.\n5. The most important new finding in this paper is that the 2024 Calama earthquake ruptures over a small scale in horizontal direction while over a wide range in vertical direction, passing across the 650-degree isotherm. Thus the authors speculate that subevents at shallow depths (E1 to E3) may be due to the embrittlement of rock caused by slab dehydration, while those at deeper depths are due to a thermal runaway, which is of great significance. In overall, there are three main evidences to support this hypothesis: first, the STF functions exhibit an obvious time gap between the subevents; second, the distribution pattern of the aftershocks is more consistent with the rupture location of the sub-events, and the third, the earthquake was also found to have crossed the 650-degree isotherm through temperature simulations. However, these evidences are relatively circumstantial evidences. Can we compare the characteristics of seismic waveforms and their spectrums of the earthquakes occurring above the 650-degree isotherm due to dehydration embrittlement, as well as those of the earthquakes occurring at greater depths due to thermal runaway, with the theoretical seismic waveforms and spectral characteristics inspired by rupture of the shallow subevents of the present earthquake (E1~E3), as well as of the deeper subevents (E4~E5) and thus give possible direct evidences?\nResponse: The reviewer suggests comparing the characteristics of the shallow and deep subevents to those of other earthquakes at similar depths and locations. While we agree that such comparisons can be informative in principle, they are unlikely feasible in this case. Moderate-to-large events at these specific depths and locations are rare, which limits the ability to capture differences in thermal structure with sufficient resolution. Moreover, even if such events were available, any subtle source characteristics would be strongly affected by complex path and site effects, especially at the frequencies used in our analysis, making reliable “fingerprinting” of waveform features extremely difficult.\nThat said, we can still present direct seismological evidence that the rupture propagated from shallow to deep depths. We achieve this by demonstrating that the distinct timing and character of the waveform packets can only be explained by a significant increase in the rupture depth. To illustrate this, we have isolated the waveform contributions from the shallow (E1-E3) and deep (E4-E5) subevents and performed a forward-modeling test, as shown in the new supplementary figure S10. The results show that the observed timing and shape of the waveform pulses can only be explained by a significant increase in rupture depth. A shallow-only rupture model is fundamentally inconsistent with the later pulses in the data.\nTherefore, our waveform modeling supports our conclusion that the rupture extended across the 650°C isotherm and transitioned into a domain where a different mechanism, such as thermal runaway, shall be required.\nWe have added this analysis as a new supplementary Figure S10 and added discussion in L238-240 to clarify that our results are supported by a robust, waveform-based modeling analysis..\nNew Figure. S10. Contributions from different subevent groups (synthetics in red) to the strong-motion waveforms (observations in black). (a) Synthetic waveforms from shallow subevents (E1-E3 between ~120-140 km) match the early pulses in the strong motion waveform data. (b) Deeper subevents (E4-E5 between ~150-170 km) explain the later, coherent energy arrivals in the data. (c) Placing E4–E5 at a shallower depth (140 km) yields synthetic waveforms that can’t fit the timing and shape of the strong motion pulses, ruling out this scenario and supporting deep depths for the later rupture.\n6. In Figure 5. I can not find thin red lines (represent the 1000oC isotherms) in this figure. In addition, the background gray circles are too small and blur to be observed, please redraw this figure.\nResponse: We revised Figures 5 and S11 to enhance the visibility of the 1000°C isotherms.\n7. From Extended Data Figure 4, the marginal probability density distributions of subevent rupture durations are wide and show a clear asymmetry, and it looks very strange to see these probability density distributions cut off abruptly in EP3 through EP5, so please check the results carefully to see if they make sense.\nResponse: The abrupt cut-off at 10 seconds in the duration distributions for subevents E3-E5 comes from the parameter search boundary. We set the maximum allowable subevent duration to 10 seconds in our inversion priors. The reason is that the relatively longer period waveforms used in our inversion (periods ≥ 5 s) have limited sensitivity to the exact duration of a subevent once it exceeds ~5 seconds. In other words, a subevent with a 7-second or a 10-second duration produces a similar waveform fit. Therefore, the result should be interpreted as these subevents have a relatively long duration, but the precise value is not tightly constrained.\nThe wider, asymmetric nature of these distributions come from several fundamental aspects of the inverse problem itself. The relationship between our model parameters and the resulting seismic waveforms is highly non-linear due to complex wave propagation. This nonlinearity naturally warps the probability space, meaning that even with Gaussian data errors, the posterior distributions for the model parameters become non-Gaussian. This effect is compounded by the limited sensitivity of our long-period data to certain parameters; for durations longer than ~5 seconds, varying the duration has a minimal impact on the waveform fit, creating trade-offs with other parameters and leading to the wide, skewed posteriors shown.\nIn summary, the observed distributions are a realistic representation of the model uncertainty. They correctly show that while the later subevents had long durations, the precise value is not perfectly resolvable with the available data. This is a well-understood feature of waveform inversion, and these asymmetries do not affect the robustness of our primary conclusions regarding the subevents' locations, moments, and the overall rupture spatiotemporal evolution.\nWe have added discussions on this in L144-147.\n8. Line 177~178, the authors mentioned that they tested both the sub-horizontal and sub-vertical planes based on the focal mechanism, I would like to see the geometries of the sub-horizontal and sub-vertical planes in three dimensional space as well as the hypocenter location of the mainshock.\nResponse: We have added the fault geometry plot in 3D in a new supplementary Figure S9.\nNew Figure S9. Three-dimensional views of the slip distribution on the (a) sub-vertical and (b) sub-horizontal fault planes. Color shading indicates the slip amplitude, with warmer colors representing larger slip. The green star indicates the hypocenter location.\n9. As mentioned in Line 256~265, the rupture speed of the 2024 Calama earthquake exceeds ~93% of the local shear wave velocity, which is speed falls in the forbidden zone between Rayleigh and shear wave speeds, so that the authors proposed that the 2024 Calama earthquake should rupture at both regular subshear and episodic supershear speeds. If so, maybe we can observe some kind of Doppler effects through seismic waveforms. If possible, the authors should try to find some evidences to verify the existence of the supershear rupture of this earthquake, such as using the method proposed by Vallée & Dunham (2012).\nResponse: We fully agree that a rupture process of this nature should produce observable Doppler effects, and we argue that this is indeed what we observe and present in our manuscript. Our primary evidence is the strong directivity effect illustrated in Figure 3, in which seismic stations located in the direction of rupture propagation recorded waveforms with systematically shorter durations and higher amplitudes compared to stations at other azimuths. This compression of the signal in the forward direction is a classic signature of a source moving at a velocity approaching, and in this case exhibiting an average rupture speed ~93% of the local shear wave velocity.\nWe appreciate the suggestion to apply the method proposed by Vallée & Dunham (2012), which constrains supershear rupture by identifying the Mach cone in the horizontally propagating surface-wave field. However, its application is most suitable for large, strike-slip earthquakes with predominantly horizontal rupture propagation, where a coherent Mach front of Rayleigh waves is efficiently generated along the surface. For the 2024 Calama earthquake, the rupture propagated primarily down-dip on a steeply dipping fault plane. In this vertical rupture scenario, the horizontal projection of rupture velocity is significantly less than the surface-wave speed, thus will unlikely generate a classic Rayleigh wave Mach cone. Consequently, we believe that for this specific fault geometry, analyzing the durations of downgoing and upgoing body waves, as we have done in our study, provides the most robust and direct constraints on the rupture velocity. We have added discussion in L338-342 to clarify this methodological choice.\n10. Figure 4. The figure caption of Fig. 4a is wrong, for Fig. 4a, the blue and green arrows are observed and predicted horizontal component displaces instead of vertical displacements, please fix this error.\nResponse: We have fixed this error in the caption of Figure 4.\nSuggested reference:\nVallée, M. & Dunham, E.M. (2012). Observation of far-field Mach waves generated by the 2001 Kokoxili supershear earthquake. Geophys. Res. Lett., 39,. L05311.\nReviewer: 3\nJia et al. investigate the rupture characteristics of the large, intermediate-depth Calama earthquake, finding that it had a steeply dipping fault plane extending ~75 km down from the hypocenter. By themselves, the combined analyses are rigorous and convincing, making this event one of the best constrained intermediate depth earthquakes of this magnitude I have read about. But what makes the earthquake potentially compelling for publication in Nature Communications is its context within the thermal conditions of the slab. As the authors explain well, the rupture propagated from the cold core of the slab, where earthquakes are known to occur, out into the hotter regions of the slab where brittle, elastic-rebound processes are unexpected. In my opinion, the manuscript is convincing and in good shape, other than some minor comments, provided below.\nResponse: Thank you for the positive evaluation. We appreciate the constructive comments and our responses are detailed below item by item. Your comments are shown in black, and our responses in blue. Line numbers refer to the revised manuscript with tracked changes.\nWhat I find myself wondering about is the thermal modeling and relationship between the Calama earthquake and temperature. [Full disclosure, I am not an expert in intermediate depth earthquakes or in the thermal state of the Nazca slab.] Fig. 5 of the manuscript shows that this section of the Nazca slab has substantially more earthquakes along the deeper cold-to-hot transition in the subducting slab than the other three settings shown. I have not looked closely at other subduction zones to see if this uniqueness extends globally, but I think this is something to probe further about the setting of the Calama earthquake since it is essential to the argument about rupture process. Could it be that the slab in the Calama region is actually slightly deeper than currently proposed, so these events (and the 2024 rupture) are still mostly within the cold slab core? Or do the authors know of alternative interpretations for the number of deeper intermediate earthquakes along the base of the slab here? Addressing this would go a long way towards solidifying the interpretive aspects of this study.\nResponse: We appreciate the reviewer’s concern regarding potential uncertainties in the slab geometry. The Nazca slab geometry used in our thermal modeling is not arbitrary–it is constrained by multiple independent datasets, including regional seismic tomography and the distribution of hundreds of thousands of relocated earthquakes that collectively delineate the shape of the slab interface. In fact, the seismicity pattern itself provides one of the most reliable indicators of slab location.\nTo make it more clear, we have updated Figure 5 (and the supplementary Figure S11) with the latest high-resolution relocated catalog from the Centro Sismológico Nacional (CSN) of Chile. This improved regional seismicity catalog aligns well with the slab structure, which strengthens the simulated thermal profiles. Notably, both the Calama event hypocenter and the historical intermediate-depth seismicity remains generally confined within the serpentinite phase transition zone near 650 °C. This confirms that the event’s location relative to this thermal boundary is well constrained.\nEven under the alternative hypothesis that the slab lies slightly deeper, the thickness of the cold core is better constrained than its absolute depth. As shown in revised Supplementary Figure S11, even when adopting an alternative geometry (e.g., from the Slab2 model) and varying the slab age (which changes the effective thickness of the cold core), the along-dip rupture of the Calama earthquake should still extend beyond the serpentinite dehydration isotherm. This confirms that our main conclusion is robust to reasonable variations in slab geometry.\nWhile the continuous rupture of the 2024 Calama event provides clear evidence of propagation beyond this dehydration boundary, interpreting the smaller, scattered background earthquakes that also appear 'out-of-bounds' is more complex. The apparent 'smearing' of this seismicity could partly reflect inherent hypocentral uncertainties. Alternatively, the subducting slab is not homogeneous; inherited faults, variations in hydration, and localized thermal anomalies may create broader and more complex zones of potential failure where this smaller-scale seismicity can occur.\nWe have added discussions on this in L263-275.\nMinor Comments:\nAn aspect of this problem the authors allude to but never address directly is the topic of strain rates. How it relates to the brittle-ductile transition is implicit in the rupture mechanisms for intermediate-depth earthquakes, but I would like to see it laid out more explicitly (cf., ll. 55-56). One thing in particular that strikes me is that at typical strain rates and temperatures of ~1000ºC, even with a mechanism like thermal runaway to localize strain it is hard to imagine having enough elastic strain energy stored for an earthquake. I understand that this is one of the key points in the manuscript. I think it would be helpful to lay it out for the reader.\nResponse: We agree that understanding how significant elastic strain energy can be stored and released seismically at temperatures approaching 1000°C is important to understanding intermediate and deep earthquakes, and a clearer explanation will benefit the reader.\nThe key is the strong dependence of a material's deformation style on the rate of applied strain. At low strain rate, the slab movement provides the loading and steadily increases the elastic strain. The materials at these temperatures will eventually deform ductilely over geological timescales through processes like dislocation creep. However, mechanisms like the dehydration embrittlement and transformational faulting of metastable olivine can trigger intermediate-depth and deep earthquake ruptures, respectively, which are localized high-strain-rate phenomena, where shear thermal runaway can be initiated. This localization creates a positive feedback loop: localized shear deformation generates heat faster than it can diffuse away, which dramatically weakens the rock (reducing its viscosity) within a narrow zone. This weakening further concentrates the strain, leading to even faster heating and runaway instability that culminates in a rapid and large-scale rupture, both for intermediate-depth and deep focus events.\nWe added discussions in L253-285 to address this point.\nI am aware of several other thermal models for subducting slabs. Qualitatively, the results here look similar, but obviously the details matter for interpreting this event. It would be good to have a comparison (even a brief one) and some references to those other models (e.g., somewhere in ll. 198-205)\nResponse: We have expanded our discussion and included a new supplementary Figure S12 to compare our thermal model with other published models for northern Chile subduction. They contain the D80 model from van Keken and Wilson (2023), and the thermal model from Wada and Wang (2009).\nAs shown in Figure S12, while there are expected differences in slab geometry arising from different model setups and cross-section locations, the key thermal characteristics are quite consistent across all three models. For example, the largest depth of the 400 degree isotherm are all around 100 km, and the thickness of the 650 degree isotherm at 150 km depth are all around or within 30 km, substantially smaller than the depth range of 78 km of along-dip rupture extent of the 2024 Calama earthquake. This means that all these independent thermal models will support our conclusion that the rupture extends beyond the 650 degree isotherm of the serpentinite dehydration and transitioned to thermal runaway.\nNew Figure S12. Comparison of the northern Chile thermal models. (a) Thermal model for the 2024 Mw 7.4 Calama earthquake source region from this study. (b) The D80_new model for northern Chile from van Keken and Wilson (2023). (c) The northern Chile model from Wada and Wang (2009). Numbers on contours show temperature in Celsius degrees.\nDetailed Comments:\nll. 26-34: I think it would be better to disentangle observations and analyses from interpretations about rupture mechanisms here, focusing first on the detailed aspects of the event, then the comparison to thermal models, and finally the connection to dehydration embrittlement or shear thermal instability\nResponse: We have restructured the abstract in L28-41 as requested, which better guides the reader through our evidence and logic.\nl. 42: What is the point the authors are trying to make in saying that “the number of earthquakes drops markedly from shallow to deep”? Could this statement be removed?\nResponse: We have removed this statement accordingly.\nll. 57-58: The references for these intermediate earthquake processes are helpful, but it might also be useful to give a brief (one sentence) description of each process too\nResponse: We added a brief description of these intermediate depth earthquake mechanisms in L65-69.\nl. 60: What exactly do the authors mean by “environmental conditions”? Temperature? Other properties?\nResponse: We have changed the phrases to specify the \"local temperature and mineral compositions\" within the slab to improve the clarity.\nll. 59-62: How would knowing things like rupture extent contribute to constraining the process? Making this explicit will help guide the reader through the rest of the manuscript.\nResponse: We added clarification in L76-78 on how rupture extent is used to constrain the dehydration process: by comparing the spatial extent of the rupture with the predicted region that allows mineral dehydration, we can examine whether the rupture is driven by dehydration embrittlement or other mechanism.\nll. 67-68: For an M 7.8, a few tens of km is not that far away. What is the location quality of the 1950 event?\nResponse: We have revised the manuscript in L83-85 to explicitly discuss the historical location uncertainty by citing the work of Kausel and Campos (1992). Their compilation of published hypocenters for the 1950 event shows a spatial variation of approximately 50-100 km.\nl. 68: I think there is no need for “However” at the beginning of this sentence.\nResponse: We removed it accordingly.\nl. 69: I am unsure about what “constrained only by analog teleseismic data” implies. That it only has a focal mechanism or moment tensor?\nResponse: Yes, we intended to convey that the available data for the 1950 earthquake are only sufficient to determine its overall source properties (i.e., single focal mechanism), not to resolve its detailed rupture process. We have revised the sentence in L86-87 to state this explicitly.\nl. 69: I think it would be good to put a period after “indicating a normal faulting earthquake” and then talk about the 2005 earthquake separately. I initially thought the points were connected, but they are not. Also, just a point to note: the concept of “normal faulting” is a little different for intermediate depth earthquakes, since the slab is dipping and surrounded by rocks (i.e., there is not a requirement for a principal stress axis to be perpendicular to a free surface).\nResponse: We revised accordingly.\nll. 85-88: I would divide this sentence into two, one describing the information obtained from seismic data and the other describing the information obtained from co-seismic displacements\nResponse: We revised accordingly.\nl. 90: Maybe I missed it, but I do not think the “core of the slab” has been defined yet\nResponse: We revised the phrase to \"cold core of the slab\" in the context to clarify that we are referring to the thermal center.\nll. 96ff: I find the description of these results really clear and easy to follow\nResponse: Thanks.\nll. 135ff: I think some references might be useful here (e.g., for EGF approaches, the stretching of the source duration caused by horizontal rupture, or the linearity of the rupture dimension/velocity)\nResponse: We added additional references in L180, 182, and 188.\nll. 153-154: I confused “coherent phases” with seismic phases at first. Is there a better word here?\nResponse: We replaced it with “coherent peaks” for better clarity.\nl. 157: Some numbers would be helpful here: How many pulses? What size asperities?\nResponse: More than ten pulses, and these small scale asperities are superimposed on the five major, lower-frequency subevents ruptured coseismically. We added more detailed and quantitative descriptions accordingly in L195-204.\nll. 174ff: I would make the geodetic analysis its own paragraph. Also, just a comment that at this depth I suspect that it is challenging to resolve a slip distribution, much like it is challenging to resolve fault geometry.\nResponse: We agree that resolving a detailed slip distribution or uniquely determining fault geometry at this depth using only surface geodetic data is inherently challenging, and neither slip distributions nor fault geometry/fault plane selection are very well constrained. However, our point here is mainly that the vertical fault version of the geodetic inversion yields plausible slip distributions and residuals compared to the GNSS that are comparable if not slightly better than those of the horizontal fault. An argument of consistency indeed; the seismological constraints are a much stronger piece of evidence for the rupture geometry.\nWe have revised the main text L223-231 to state these limitations and to clarify that we consider this a supporting consistency argument, with the seismological analysis providing much stronger evidence for the final fault geometry.\nl. 189: What is a “complete half-section”?\nResponse: Our intention was to describe the rupture that penetrates the lower half way to the base of the slab lithosphere. We have revised in L237 by using “if not reached the base of the subducting lithosphere” instead.\nll. 193-197: I am unsure of the argument being made about aftershocks and temperature here. Sorry my comment is not more helpful.\nResponse: We have rephrased L244-250 for better clarity: “Compared with other intraplate earthquakes, 2024 Calama had a significant number of aftershocks, which suggest that its rupture occurred in a colder, brittle regime. This is more analogous to shallower events such as the 2007 Michilla Mw 6.7 and 2020 Calama Mw 6.9 earthquakes, rather than events in high-temperature zones (> 600°C) where aftershocks are typically suppressed.”\nll. 236-237: If this were the case, we might expect to see a temporal relationship between the megathrust earthquake cycle and intermediate depth earthquakes. I know this has been argued for the region further north (where the 2005 Tarapaca occurred), but I am unaware of a similar argument here. Also, in very large megathrust earthquakes, these stresses are recovered, so it might not be a long-term effect. I guess my point is that I am unsure about whether high coupling is a “key factor” or not.\nResponse: Thank you for pointing this out. Our intention is not to argue for a direct, cyclical triggering mechanism. Rather, we aim to frame the high degree of plate coupling as a critical component of the overall tectonic environment that contributes to the long-term stress state within the downgoing slab. We have revised in L309-310 to present this as a potential contributing or background factor rather than a definitive \"key factor”.\nll. 244ff: It seems that there is a little jumping back and forth between intermediate and very deep earthquakes. The authors should be very explicit throughout this paragraph.\nResponse: We have revised the paragraph for better clarity.\nll. 349-350: What is the specific criterion for no longer adding subevents? Is it based on misfit alone? Or is there something like the Akaike Information Criterion applied to select the best number of subevents?\nResponse: This is a good point and similar to question #5 raised by Reviewer #1. To consolidate our response and avoid repetition, we have provided a detailed, comprehensive answer under that response.\nIn brief, our selection of five subevents was based on two criteria: (1) an L-curve trade-off analysis showing that five subevents is the optimal point before misfit reduction becomes marginal (now shown in new supplementary Fig. S2), and (2) a qualitative comparison confirming the 5th subevent is necessary to fit key features in the strong-motion waveforms (Fig. S2).\nFor the full technical explanation and new figures, we respectfully refer you to our detailed response to Reviewer #1, Q5. We have also revised the manuscript to make these details clear to readers.\nl. 371: What is the “width” of the probability density function defined by? Is a Gaussian fit to it?\nResponse: The \"width\" we refer to is the 95% confidence interval of the posterior probability density function (PDF), defined by the 2.5 and 97.5 percentile values of the distribution. For the posterior distribution, we did not assume a Gaussian distribution. For our source inversion problem, some parameters, including the source durations and depths can have naturally asymmetrical PDFs where one side of the distribution is better constrained from data than the other. A simple Gaussian fit would not appropriately capture this.\nWe have added these details in L465-467 for improved clarity.\nll. 407-409: As I said above, I suspect that the fault slip determined by static co-seismic displacements is poorly constrained\nResponse: We have rephrased the context to better clarify its limitation.\nFig. 1: In panel 1a, “Magnitude” is spelled wrong.\nResponse: We corrected accordingly.\nFig. 5: The 1000ºC isotherm is very hard to see. Could the authors make it slightly bolder?\nResponse: We have revised Figures 5 and S11 to improve the visibility of the 1000ºC isotherm.\nFig. 6: Just a note: the scale of the overriding plate makes it seem as if everything is occurring at fairly shallow depths\nResponse: We have added a reference depth on Figure 6."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed my comments and I believe also the other 2 reviewer's comments. I think the paper is in good shape.", "reviewer_2": "I have carefully read the rebuttal letter and the revised manuscript. All of my concerned problems have been resolved well, in addition, the scientific problems mentioned by the other two reviewers have also been answered and resolved. I think the quality of the revised manuscript has been upgraded significantly, so that the paper is suitable for being accepted by NC. I recommend to accept the paper at the present form.", "reviewer_3": "I am pleased with the authors' comprehensive responses and revisions to my comments and those of the other reviewers. I would be glad to see this study published."}, "author_response": ""}]

Decision letter
Editors’ comments:
You will see that while the reviewers are interested in your study, they also raise concerns. These include the interpretation of the rupture of asperities and the rupture propagation from the cold to hot regions of the slab. We find that these concerns limit the strength of the study, and therefore we ask you to address them with additional work.