Variation in the thermal and dehydration regime below Central America: Insights for the seismogenic plate interface

Summary Slow earthquakes predominant in Costa Rica indicate unstable faulting of segmented Central American megathrusts, but the recurrence of episodic tremors and slips reported to precede a giant earthquake remains still enigmatic. The underlying mechanism is related to the variation in the coupling along the heterogeneous subduction interface which is poorly understood. In this study, we used up-to-date 3D thermal modeling to provide insights into the along-strike variation in the thermal state and hydraulic distribution beneath the Central American subduction zone. Our results show that the subducted Cocos Plate is much warmer than previously estimated, and the slab geometry exhibits remarkable perturbations along the trench. We found that the regions of large dehydration rate along the slab are consistent with the seismicity occurrence depth beneath the Moho. Below the Nicoya Peninsula and the Guatemala-Nicaragua segment of megathrusts, fluids derived from subducted slab result in increased pore fluid pressures and subsequent recurrence of slow slip events and regular earthquakes.


INTRODUCTION
depth as the deep portion of the 2007 SSE. Consequently, Baba et al. 10 analyzed the relationship between slow earthquakes and large earthquakes in the Nicoya Peninsula segment of the Central American subduction zone and proposed that slow earthquakes may occur in the shallower and deeper extents of the rupture regions of large earthquakes. However, the reason for the differential distribution of these slow earthquakes remains enigmatic.
The occurrence of an earthquake in a subduction zone is intimately related to the thermal structure of the subducting slab and the overlying mantle wedge. 22 Thus, to investigate the relationship between fast and slow earthquakes in the Central American subduction zone, the thermal regime of the subducting Cocos Plate must be accurately known. Recent research has demonstrated that the 3D thermal model constrained by the various observations of surface heat flow shows advantages over the existing 2D thermal models in predicting the along-strike variation in the thermal state of the subducting slab (e.g., Qu et al. 23,24 ). Accordingly, in this study, we estimated the thermal state and water content distribution of the Central American subduction zone based on a 3D thermomechanical model. Then, we compared the modeled results with the observed foci of regular interplate earthquakes, large earthquakes, and slow earthquakes, and we analyzed which main mechanisms could be responsible for generating fast and slow earthquakes beneath Central America.

3D temperature structure of the subducting Cocos Plate
Based on the 3D numerical simulation, we calculated the intraplate temperatures of the subducting Cocos Plate at four isodepth contours (measured downward orthogonal to the slab surface) of 0, 4, 8, and 16 km (Figures 4 and S3). In addition, we relocated the fast and slow The background color indicates the surface topography (ETOPO 4 ). The cyan dashed box delineates the study region. The yellow arrows denote the direction (orientation) and velocity (length) of subduction, and the red barbed lines mark convergent plate boundaries (with the teeth on the overriding plate). The red curves represent the isodepth contours of the subducted Pacific plate with a contour interval of 20 km (Slab2.0 5 ). Red triangles indicate arc volcanoes. 3 The solid circles indicate the epicenters of all earthquakes from January 2000 to December 2010 (IRIS 6 ) and M > 5.5 earthquakes from January 1900 to December 2000 (Centennial 7 ). Earthquake depths are indicated by the circle colors. The yellow ellipses with dashed lines represent the distribution of slow earthquakes below the Nicoya Peninsula, including the low-frequency earthquakes (LFEs 8 ), large slip areas of the 2007 slow slip event (SSE 9 ), and very-low-frequency earthquakes (VLFEs 10  iScience Article earthquakes in the Cartesian coordinate system projected from spherical coordinates; these included all regular earthquakes and slow earthquakes. 6,7,19 Identified slow earthquakes have been located mainly offshore El Salvador 19 and below the Nicoya Peninsula. 8,9,10 The slab surface temperatures in the source regions of deep SSE and LFEs vary between 300 C and 900 C, whereas shallow SSE and VLFEs occur at lower surface temperatures varying between 200 C and 250 C ( Figure 4A). Regular interplate earthquakes occur mainly at depths of 20-80 km at the slab surface with surface temperatures ranging from 300 C to 900 C. The surface temperature of the subducting Cocos Plate at a depth of 10 km (from Earth's surface) is mostly approximately 200 C, but the surface temperatures offshore northern Costa Rica are slightly higher than those elsewhere due to the younger age of the slab beneath this region ( Figure 4A). The surface temperature of the incoming plate at a depth of 100 km below Nicaragua is approximately 1,100 C. In contrast, the subducted slab beneath El Salvador is warmer, with surface temperatures exceeding 1,100 C at 100 km depth ( Figure 4A). At the Moho (30 km-depth plate interface), temperatures predominantly in the range of 300 C-500 C control the serpentinization occurring within the mantle wedge. Meanwhile, below central Costa Rica, the surface temperature of the subducted slab is lower, decreasing to 300 C -500 C (mantle wedge) at the depth range of 30-60 km, colder than that below El Salvador. The thermal structure of the slab is associated with the geometry of the subducting Cocos Plate. In other words, the variations in the geometry of the slab, such as its dipping angle variations, occur in association with strong temperature variations across strike ( Figures 4A and 4B).
According to the calculated results ( Figure 4A), we also observed that the subvolcanic slab temperature was as high as >1,100 C on the plate interface. Since the Cocos Plate is warm and thin, the intraslab temperature increases gradually with increasing depth. For example, at the mantle wedge immediately underlying the Moho (30 km) below the central Nicoya Peninsula (Figure 4A), the calculated temperature is 400 C, but at greater depth (16 km into the slab), the slab temperature increases to 550 C -600 C ( Figure 4B). Such elevated intraslab temperatures indicate the absence of a cold center in the slabs of warm subduction environments like the Central American subduction zone. Furthermore, we noticed a sudden drop in temperature within the incoming plate below central Costa Rica compared to other regions at the same depth ( Figure 4B). It indicates that factors apart from slab geometry such as hydrothermal circulation near the trench, slab age, and convergence rate influence the temperature variation in the thermal structure of the slab here.  16 The seafloor ages are derived from EarthByte. 17 Red triangles indicate arc volcanoes. 3 The yellow dashed lines represent the slow earthquakes below the Nicoya Peninsula. 18 Red and yellow polygons indicate the interplate seismogenic zone (10-40 km depth) and intraplate seismogenic zone (>40 km depth), respectively. 14 ll OPEN ACCESS iScience 26, 107936, October 20, 2023 3 iScience Article 3D water content variation along the strike Previous seismic inversion results found that the water content of the subducting Cocos Plate beneath Central America is varying along strike due to differences in crustal thickness, subduction velocity, and dip angle. 26 To estimate the variation in water content of the slab beneath southern Guatemala to central Costa Rica, we assumed that the rocks were water saturated and calculated the 3D distributions of the slab water content (maximum saturation) ( Figure 5) and dehydration ( Figure 6) according to the facies diagrams of mid-ocean ridge basalt (MORB) 27 and ultramafic rocks 28 in the subducting Cocos Plate. Then, similar to the temperature analysis, we computed the slab water content and dehydration distributions on intraslab isodepth contours of 0, 4, 8, and 16 km below the subducting Cocos Plate interface. The subducting Cocos Plate is assumed to be composed of an upper 7 km thick layer of MORB and a lower layer of ultramafic rocks. 28 At the shallowest depths (<20 km, Figure Figure 6B, ultramafic layer). Notably, deep SSE, LFEs, and most fast earthquakes occur in the region where the slab dehydration rate is > 0.01 wt %/km. By contrast, although a handful of shallow SSE and VLFEs are located in the updip region of the dewatering zone (e.g., offshore Nicoya Peninsula), only limited slab dehydration and fluid release occur there according to our calculations.

Comparison with previous thermal structure results
Multiple 2D thermal models have been developed and applied to the Central American subduction zone beneath Nicaragua and Costa Rica, and the slab thermal structure has been estimated along several profiles (e.g., Peacock et al. and Wada et al. 22,29 ). Wada and  3 The white dashed line indicates the potential slow earthquake area predicted by this study. The blue polygon shows the afterslip area of the 2012 Mw 7.3 tsunami earthquake offshore El Salvador. 19 The yellow rectangle indicates the aftershock area of the 1992 Mw 7.6 Nicaragua earthquake. 25 The yellow ellipse with a dashed line shows the distributions of very-low-frequency earthquakes 10 and the shallow area of the 2007 slow slip event. 9 The blue ellipse with a dashed line represents the distribution of lowfrequency earthquakes. 8 23 suggested medium-high heat flows (50-140 mW/m 2 ) beneath northeast Japan and the Cascadia forearc, indicating that the forearc is not as cold as described by the slab-mantle (or interplate) decoupling model. 29 In Central America, we also found high heat flows (>80 mW/m 2 ) at the trenchward forearc (e.g., offshore Costa Rica) (Figures 2 and S1), which may indicate that the subducting slab surface temperature is warmer than previously thought. Therefore, we adopted our model to reevaluate the temperature range of the cold mantle wedge (<500 C) and constrained the bottom of the cold nose to depths of 40-60 km (Figures 7B and 7D).
The thermal structure of the subducting slab depends on many parameters, including the age of the subducting lithosphere, convergence rate, geometry, rate of shear heating, and so on. 31 Maunder et al. 32 found that the slab age and speed have the most important effects on the slab surface temperature. Following the isoviscous mantle-wedge rheology and the olivine mantle-wedge rheology, Peacock et al. 22 predicted slab surface temperatures beneath Nicaragua and Costa Rica of approximately 620 C and 800 C at 3 GPa (100 km depth), respectively. Correspondingly, we calculated the temperature variations along the same three profiles as in Peacock et al. 22 passing through volcanic areas in Nicaragua, northern Costa Rica, and central Costa Rica ( Figure S2A). Our thermal results reveal that at the same depth (100 km), the slab surface temperatures beneath Nicaragua and northern Costa Rica are 1,000 C ( Figures S2B and S2C), while the slab surface  Figure S2D). Compared with previous 2D models, our model is jointly constrained by a variety of additional morphological controlling factors such as the slab curvature, which is relatively sensitive to modest variations in the age of the subducting lithosphere and convergence rate (e.g., [33][34][35] ). Overall, we find that the surface temperature of the subducting Cocos Plate beneath Central America is much warmer than previous studies have suggested, which is consistent with the surface heat flow and seismicity observations. Conductive lithospheric cooling models have predicted heat flows of 95-120 mW/m 2 for 18-24 Myr seafloor. 36 However, Fisher et al. 37 observed particularly low seafloor heat flows (20-40 mW/m 2 ) near the Middle American Trench and inferred the occurrence of vigorous hydrothermal circulation in the shallow basement that extracts heat. Vigorous fluid circulation can affect subduction zone temperatures and thus significantly change the thermal structure of subducting plate. [38][39][40] Hass and Harris 41 incorporated hydrothermal circulation into their thermal models, and their results found that the incoming plate beneath central Costa Rica is warmer near the trench and colder away from the trench. From our results, since the changes in temperature of the slab beneath central Costa Rica shown in Figure 4 are consistent with their findings, we thus infer that hydrothermal circulation also contributes to the variation in thermal structure of the subducting plate beneath Central America along strike.

Distribution of fast-slow earthquakes in the central American subduction zone
We further investigated the ambient conditions required for fast and slow earthquakes to occur at various depths and sought to discover the interconnection among them. We observed that, except for a few earthquakes scattered at shallow depths, most of these earthquakes are iScience Article distributed at absolute depths of 20-80 km, which is related to the dewatering zone of >0.01 wt %/km ( Figure 6). Fluids are typically generated from within the slab by multistage metamorphic reactions requiring high-temperature and high-pressure conditions associated with depth, 42,43 but our results regarding the water content variation suggest that the upper depth limit for the dewatering zone on the upper slab surface is 20 km, which may be attributable mainly to the fact that the shallow slab is characterized by such high temperatures. Moreover, the widely accepted mechanism for the generation of intermediate-depth (50-300 km) earthquakes is dehydration embrittlement from the breakdown of hydrous minerals within the subducting slab. 28,44 We consider that dehydration embrittlement could also have implications for shallow (20-50 km) seismicity in Central America, which is consistent with the observed dewatering zone of the slab surface. Therefore, we propose that warm-plate-environment-controlled dehydration may be the main driving mechanism for fast and slow earthquakes along the Central American subduction zone.
Previous studies have noted that the generation mechanisms of shallow and deep SSEs may be similar. 21 At shallow depths, fluids can be released by the compaction of subducted sediments atop the slab. 43,45 Rü pke et al. 46 modeled the fluid release beneath Nicaragua and Costa Rica, and their results showed that the sediments released 75% of the chemically bound water during shallow (depth <50 km) dewatering. Our results suggest that deep SSE and LFEs (30-50 km) occur on the slab surface in the dewatering zone, while shallow SSE and VLFEs (6-10 km) occur elsewhere. Outerbridge et al. 9 suggested that the occurrence of SSEs at a shallow depth is attributed to the updip transition from stick-slip to stable sliding. However, this transition has difficulty occurring at such shallow depths considering the change in interplate friction conditions. Thus, we stipulate that fluids derived from the compaction of subducted sediments promote the occurrence of shallow SSEs. Overall, beneath the Nicoya Peninsula, fluids were released in different ways and at different depths, ultimately leading to an increase in the pore fluid pressure, which promoted SSEs on the megathrusts atop a young and hot segmented slab beneath Costa Rica with the seafloor age of 16-25 Ma. Furthermore, partially released fluids are conjectured to be stored in the crust beneath the Nicoya Peninsula. 47 Stored fluids accumulated over a long time and triggered ruptures of earthquakes. This may explain why the Nicoya Peninsula lacks regular interplate earthquakes but experienced several large earthquakes (Figure 3, yellow stars and white stars). When the subducting plate interface is in a fluid-rich environment characterized by high pore pressures, the fluid overpressure reduces the effective normal stress and weakens the plate interface, thereby facilitating brittle failure and ultimately causing the occurrence of SSEs. 48,49 Beneath Central America, the presence of the observed dewatering zone on the plate interface from Guatemala to Costa Rica provides similar support for the generation of slow earthquakes. Previous studies located the aftershock area of the 1992 Mw 7.6 Nicaragua earthquake and suggested that it was a slow tsunami earthquake associated with subducted sediments 18,25 ( Figure 4A). Furthermore, Geirsson et al. 19 estimated the seismic deformation from the 2012 Mw 7.3 tsunami earthquake offshore El Salvador and found significant afterslip, causing the event to manifest as exhibiting slow slip over a large region ( Figure 4A). The aftershock area and afterslip area indicate that the plate contact offshore from El Salvador to Nicaragua is weakly coupled and unstable, similar to that of Costa Rica. Consequently, we suggest that the northern-central part of the Central American subduction zone has the potential for slow earthquakes to occur, particularly in the Guatemala-Nicaragua segment of the heterogeneous megathrust ( Figure 4A, surrounded by a dashed white line).
Considering that slow earthquakes require long-term observations over several years to decades, 20 researchers deployed multiple observational networks in the northern-central part of Central America that can provide continuous data. [50][51][52] However, compared to the dense instrumental network located in the Nicoya Peninsula, the lack of near-trench geodetic monitoring in the Guatemala-Nicaragua segment makes detecting slow earthquakes extremely challenging. 53 Moreover, the influence from regular interplate earthquakes ( Figure 4A, M > 3), the noise in GPS data, 20 and instrumental responses 10 also make it difficult to distinguish SSEs and LFEs from the data. In this case, previous studies have been more focused on using models to predict the probability of slow earthquake occurrence. For example, McLellan and Audet 54 predicted the occurrence probability of SSEs in Central America and showed that their occurrence (durations of days to a few weeks) is the highest in northern Central America.
Furthermore, as the subducting plate changes curvature in three dimensions, fluid pore pressure can hence be affected and follow areas where the surrounding stress (including resultant brittle failures) allows them to flow. 55,56 For instance, some 2D and 3D flexural bending studies have found that the stress and deformation of the subducting slab can exert influence on the temperature of the slab and the fluid transport. [57][58][59] In addition, the potential impact from the overriding plate and its thermomechanical controls cannot be ruled out in terms of stress coupling. If the stress and strain evolution of the slab-mantle system is further investigated in the modeling, the composite of subduction regime could be obtained at a higher accuracy, and this work will be left for future studies.

Conclusions
According to 3D thermomechanical modeling of the Cocos Plate subducting beneath Central America from Guatemala to Costa Rica, the following conclusions are drawn.
(1) The surface temperature of the incoming plate beneath Central America is warmer than previously thought, and the cold mantle wedge (<500 C) below the Central American subduction zone is constrained at depths of 40-60 km. (2) Warm-plate-environment-controlled slab dehydration appears to be the mechanism responsible for the recurrence of fast and slow earthquakes along the heterogeneous Central American megathrusts at depths of 20-80 km. (3) Fluids derived from subducted sediments at shallow depths and slab metamorphism at greater depths below the Nicoya Peninsula and the Guatemala-Nicaragua segment of megathrusts result in the increased pore fluid pressure and promote the occurrence of slow earthquakes over a broad depth range.

Limitations of the study
In this study, we found that the surface temperature of the Cocos Plate is warmer than previously estimated. As the plate subducted, warmplate-environment-controlled slab dehydration led to the release of more fluids and drove the occurrence of fast-slow earthquakes in the Central American subduction zone. The limitations of the current study are mainly attributed to the modeling settings; for example, the aquifers are not taken into account, which may result in low surface heat flow at the Nicoya Peninsula. This causes some uncertainty in the calculated temperature of the shallow part (depth <10 km) of the subducted plate. On the other hand, there is some uncertainty in the calculation of plate water content variation. We utilized a water content of 15 wt % for the subducting oceanic plate because harzburgite is the dominant rock type in the uppermost oceanic mantle, which is largely serpentine at lower temperatures. 28 However, the estimates remain uncertain because serpentinites at low temperatures (50 C -300 C) are suggested to be mainly composed of lizardite and magnetite with brucite, but chlorite is not common. 60

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

ACKNOWLEDGMENTS
We thank P. Tackley for sharing the original Stag3d code adapted for this study. We thank IRIS for providing the earthquake catalog at http:// geoserver.iris.edu/. Some figures were created using Generic Mapping Tools software developed by Wessel and Smith. 61 We are grateful to David Yuen, Muriel Gerbault, and two anonymous reviewers for their constructive comments that allowed us to improve this manuscript. We also acknowledge the support of computing facilities provided by the Institute of Tibetan Plateau Research in this study, including the Sugon supercomputing platform and workstations.
Funding: This study is financially supported by the National Natural Science Foundation of China (41988101-0106), the Pioneer Hundred Talents Program of the Chinese Academy of Sciences, and the Second Tibetan Plateau Scientific Expedition and Research Program (2019QZKK0708).

AUTHOR CONTRIBUTIONS
Y.J. conceived the original idea, designed the 3D thermomechanical codes, and performed the numerical experiments. R.Q. and Y.J. expanded the numerical study and wrote the manuscript. L.L., W.Z., Y.Z., C.X., S.Y., and H.F. provided comments that improved the manuscript. All authors discussed the results and interpretations and participated in writing the paper.