3D thermal structural and dehydration modeling in the southern Chile subduction zone and its relationship to interplate earthquakes and the volcanic chain

Abstract

In southern Chile, the Nazca plate is subducting beneath the South American plate. This region was struck by megathrust earthquakes in 1960 and 2010 and is characterized by the existence of a volcanic chain. In this region, we modeled a three-dimensional thermal structure associated with the subduction of the Nazca plate by using numerical simulations. Based on the obtained temperature distribution, we determined the updip and downdip limit temperatures for the region ruptured by these two megathrust earthquakes. In addition, the distributions of water content and dehydration gradient were calculated by using appropriate phase diagrams and compared with the location of the volcanic chain. As a result, we infer that the coseismic slip of the 2010 Mw8.8 Maule earthquake occurred only at temperatures lower than and around the 350 °C isotherm that resembles the beginning of the brittle‒ductile transition. We also deduce that the rupture of the 1960 Mw9.5 Valdivia earthquake propagated up to the 450 °C isotherm because the magnitude was considerably large and the young hot plate subducted near the Chile Ridge. In addition, the hydrous minerals in the turbidites, MORB and ultramafic rocks released fluids via dehydration reactions, and dehydrated water migrated upward almost vertically, decreasing the melting point of the mantle wedge and contributing to the formation of the volcanic chain.

Introduction

In the southern Chile subduction zone, between 33°S and 44°S, the Nazca plate is subducting beneath the South American plate in an east‒northeastward orientation at ~ 6.2 cm/yr (Altamimi et al. 2016) (Fig. 1). This region has historically been struck by great earthquakes (e.g., Kanamori 1977; Contreras-Reyes and Carrizo 2011). For example, the great 1960 Mw 9.5 Valdivia earthquake was the world's largest instrumentally recorded earthquake; during this earthquake, the generated tsunamis reached great distances within the region around the Pacific Ocean (e.g., Ho et al. 2019). The 2010 Mw 8.8 Maule earthquake also occurred in this region and generated a devastating tsunami (e.g., Moreno et al. 2012; Lin et al. 2013). In addition, a volcanic chain called the Southern Volcanic Zone extends from 33°S to 46°S and is almost parallel to the trench. The top of the slab surface is located ~ 250–350 km horizontally from the trench at depths greater than 80 km below the volcanic chain.

Fig. 1

Tectonic map of the southern Chile subduction zone. The map shows the area within the black dashed boxed area in the inset. The area surrounded by blue dashed lines is the horizontal projection of the 3-D model region. The thick black barbed line represents the plate boundary at the Earth’s surface (Bird 2003), and the thin black line denotes the upper surface of the subducted Nazca plate with a contour interval of 20 km (Hayes et al. 2018). The yellow lines indicate the current seafloor age of the Nazca plate (Seton et al. 2020). The red solid line and orange solid line represent the coseismic slip distributions of the 1960 Valdivia earthquake (Mw9.5) and 2010 Maule earthquake (Mw8.8) with contour intervals of 10 m and 5 m, respectively (Ho et al. 2019; Moreno et al. 2012). The red solid triangles denote volcanoes. The purple arrows represent the convergence rates of the Nazca plate with respect to the South American plate calculated based on ITRF2014 (Altamimi et al. 2016). The bathymetry and topography data are obtained from ETOPO1 (Amante and Eakins 2009). The map was created by using Generic Mapping Tools (GMT) (Wessel and Smith 2016)

Herein, we briefly summarize previous studies of thermal structural models in the Chile subduction zone. Wada and Wang (2009) constructed a 2-D thermal model, considering slab and mantle wedge decoupling to explain the surface heat flow. Spinelli et al. (2016) constructed thermal structural models within 2-D vertical cross-sections at two profiles along the subduction of the Nazca plate with seafloor ages of 6 Myr and 14 Myr and showed that hydrous minerals in the oceanic crust contribute to temperature, dehydration and melt. Valdenegro et al. (2019) constructed a 2-D thermal structural model at 33°S and indicated that most of the earthquakes in the past 20 years have occurred at temperatures lower than the 350 °C isotherm, which corresponds to the brittle‒ductile transition region. Ji et al. (2019) constructed a 3-D thermal structural model in the north–central Chile subduction zone (19°S–33°S) to investigate the temperature range in which megathrust and slow earthquakes occur. They concluded that the temperature threshold between megathrust and slow earthquakes is 400–600 °C, the threshold for the dehydration gradient in the subduction direction is 0.03–0.05 wt.%/km for slow earthquakes, and megathrust earthquakes occur at temperature and dehydration values that are less than these parameters. Rodriguez Piceda et al. (2022) calculated the 3-D thermal distribution in the south central Andes (29–39°S). They numerically solved the steady-state heat conduction equation at shallow depths and calculated temperatures at high depths from the results of S-wave seismic tomography. By combining these two temperatures, they calculated the surface heat flow and identified the differences between the values in the orogenic zone and the values in the forearc due to the shallow cold slab. Iwamoto et al. (2023) constructed a 3-D thermomechanical model associated with simultaneous subduction of the Nazca (NZ) and Antarctic (AN) plates near the Chile Triple Junction (CTJ). The results show that the current temperature distributions on the upper surface of the slabs are higher closer to the Chile Ridge, and the AN plate has a distribution of elevated temperatures relative to the NZ plate at the same depth. They also calculated the water content and dehydration gradient from the temperature distribution near the upper surface of the slab and discussed their relationship to the distributions of volcanoes and tectonic tremors.

Studies that compare the 3-D thermal structure and dehydration distribution in the southern Chile subduction zone with megathrust earthquakes, interplate earthquakes and the volcanic chain are lacking. A 3-D model can address differences in mantle flow and temperature in the along-arc direction, which cannot be reproduced by the 2-D model, indicating that the former is a much more realistic model than the latter. In addition, a 3-D model can capture the temperature distribution and dehydration gradient distribution at the plate boundary where earthquakes occur in an areal manner, which is more suitable for capturing phenomena that extend over the horizontal plane in the along-arc direction than the 2-D models employed in previous studies. Therefore, in this study, we constructed a 3-D thermomechanical model in the southern Chile subduction zone, which is associated with the subduction of the Nazca plate, and determined the temperature range of the slip areas of the 1960 Valdivia earthquake and 2010 Maule earthquake. Furthermore, we calculated the water content near the plate boundary based on the obtained thermal structure and phase diagrams of hydrous minerals and discussed the relationship between the distribution of dehydration and the Southern Volcanic Zone.

Methods and model

We modeled a thermal structure in the Chile subduction zone that has been associated with the subduction of the Nazca plate over the last 20 Myr using convergence rates and obliquity angles from tectonic reconstruction models and modern space geodesy calculations (Müller et al. 2019; Young et al. 2019; Cao et al. 2020). The grid intervals were approximately 10 km × 10 km × 3.3 km (x × y × z), which were reasonably spaced for the target earthquakes and the volcanic chain (Additional file 1: Fig. S1). The Nazca plate was subducted along a prescribed guide based on the Slab2 geometry (Hayes et al. 2018). The models reached a steady state after 20 Myr of subduction. In the energy equation, thermal conduction, viscous dissipation, radioactive heating, and frictional heating were considered (Ji and Yoshioka 2015). In accordance with previous studies (Honda 1985; van Keken et al. 2002; Wada and Wang 2009), we examined a case of incorporating a rigid mantle wedge into the model to reproduce the decoupling at the plate interface. The effective friction coefficient (μ’) at the plate boundary and the presence or absence of a rigid mantle wedge were tuned to match surface heat flow observations (Additional file 1: Fig. S2) (Abbott 2008; Delisle and Ladage 2002; Grevemeyer et al. 2003, 2005, 2006; Hamza et al. 2005; Hamza and Muñoz 1996; Mas et al. 2000; Muñoz and Hamza 1993). For each model, we calculated the weighted root mean square (RMS) values of residuals between the calculated and observed heat flow (Additional file 1: Table S2). Based on the result minimizing RMS of heat flow, we obtained the optimal model when μ’ was zero without the rigid mantle wedge. The parameters of this model were approximately consistent with the distribution of the Curie point depths by Li et al. (2017), which covered a wider area than the heat flow values. Therefore, we selected this model as the optimal model in this study. We refer the reader to Supplementary Information for further details on data, methods and numerical model (Additional file 1: Texts S1 to S3) as well as to the distributions of heat flow, Curie point depth and temperature at the plate boundaries for other models (Additional file 1: Texts S4 and S5).

Based on the optimal thermal structural model (0 Ma) that was recently obtained from our numerical simulation, we applied the temperatures and depths along the geometry of the slab surface to the appropriate phase diagrams and calculated the water content distribution. When calculating the water content, the slab was divided into three layers (Additional file 1: Fig. S1), and a different phase diagram was applied to each layer (Additional file 1: Fig. S2). The oceanic sedimentary layer was defined as the vertical layer from the plate boundary to a depth of 2 km, for which the turbidite phase diagram was applied (van Keken et al. 2011). The oceanic crust was defined at a depth range of 2 km to 7 km below the plate boundary, for which the MORB phase diagram was applied (Tatsumi et al. 2020). The thickness of the sedimentary layer and the oceanic crust described above were determined by seismic reflection images in the Chile subduction zone (Olsen et al. 2020; Ramos et al. 2018). The slab mantle was defined as the layer from a depth of 7 km below the plate boundary to the base of the oceanic plate. The ultramafic rock phase diagram was employed for the slab mantle and the bottom of the mantle wedge to a height of 2 km above the plate interface (Tatsumi et al. 2020). Based on the water content distribution, we calculated the dehydration gradient, which was the difference in the water content per unit length (km) along the subduction direction (Suenaga et al. 2019). Here, the area where the dehydration gradient was greater than 0.1 wt.%/km was referred to as a dehydration zone.

Results

Thermal structure

The comparison between the observed heat flow and calculated heat flow for the optimal model (μ’ = 0.000 and without rigid mantle wedge) is shown in Fig. 2. The distribution of the calculated heat flow near the trench is higher to the south because the age of the Nazca plate is younger to the south (Fig. 2a and 2e). Regarding the change in the calculated heat flow values from the trench to the inland, the values at approximately x = 0–200 km are small due to the subduction of the cold oceanic plate (Fig. 2b, c and d). Incidentally, the value at approximately x = 0–150 km in the optimal model is lower than those in the models with μ’ = 0.005 because frictional heating at the plate boundary does not act in the optimal model. The value at approximately x = 150–250 km in the optimal model is higher than those in the models with rigid mantle wedge because hot mantle material enters the corner of the mantle wedge. On the inland side at approximately x = 250–450 km, the heat flow was nearly constant as the shallow part of the model domain is covered by the stratified upper and lower crusts (Fig. 2f and g). We compared the Curie point depth distributions calculated from our optimal model and from those in Li et al. (2017), which resulted in approximately the same trend. The details are shown in Additional file 1: Text S5 in the Supplementary Information.

Fig. 2

a Comparison between the observed and calculated heat flows for the optimal model. The color contours are the spatial distributions of the calculated heat flow obtained in this study with a contour interval of 10 mW/km². The colored squares represent the average heat flow and average locations of the observation points, which are located within 1 degree of each latitude and longitude in the same reference papers (Abbott 2008; Delisle and Ladage 2002; Grevemeyer et al. 2003, 2005, 2006; Hamza et al. 2005; Hamza and Muñoz 1996; Mas et al. 2000; Muñoz and Hamza 1993). The white solid lines are the locations of the profiles shown in (b–g). The map was created by using generic mapping tools (GMT) (Wessel and Smith 2016). a Comparison between observed and calculated heat flows along profile B-B′ in (a). The solid black line is the optimal model. The green, blue and pink dashed lines are the model with μ’ = 0.005 without rigid mantle wedge, the model with μ’ = 0.000 with rigid mantle wedge and the model with μ’ = 0.005 with rigid mantle wedge, respectively. The colored diamonds are all observations (not average) within a one-sided width of 30 km. c Same as (b), except for profile C–C′. d Same as (b), except for profile D-D′. e Same as (b), except for profile E‒E’. f Same as (b), except for profile F-F′. g Same as (b), except for profile G-G’

As the grid interval in the depth direction was approximately 3.3 km in this numerical simulation, we performed a linear interpolation of the temperature between each grid point in the depth direction and obtained the temperature distribution along the plate boundary in the optimal model (Fig. 3). The temperature distribution at the plate boundary tended to be strongly dependent on the geometry of the subducting Nazca plate. Therefore, the temperature changes along the subduction direction, namely, the dip direction, were larger in the southern part of the model domain than in the northern part because the former dip angle was relatively large. Furthermore, the age of the oceanic plate at the trench decreased toward the southern part (Fig. 1), resulting in slightly higher temperatures at the plate boundary in the south than in the north at the same depth.

Fig. 3

Horizontal projection of the temperature distribution at the plate boundary in the present day (0 Ma) obtained in this study. The contour interval is 100 °C. The temperature distribution is shown only in the region where the upper surface of the oceanic plate is shallower than the bottom of the model (200 km depth). The thin black box represents the horizontal projection of the model region. The red solid line and orange solid line denote the coseismic slip distributions of the 1960 Valdivia (Mw9.5) earthquake (Ho et al. 2019) and 2010 Maule (Mw8.8) earthquake (Moreno et al. 2012), respectively, with contour intervals of 5 m. The map was created by using Generic Mapping Tools (GMT) (Wessel and Smith 2016)

Distributions of the water content and dehydration gradient

The water content distribution at present (0 Ma) for the optimal thermal model is shown in Fig. 4.The distribution of the dehydration gradient is shown in Fig. 5. Notably, Fig. 4 shows the maximum water content contained in the hydrous minerals that can be calculated from the phase diagram. Therefore, the actual water content and dehydration values may be less than the values obtained in this study. However, it is likely that the dehydration reaction occurs at the proposed location because the maximum water content significantly changes at the temperature and pressure under which the dehydration reaction occurs. Therefore, we suggest that the location of dehydration in Fig. 5 is reasonable.

Fig. 4

a Horizontal projection of the current water content distribution on the plane 2 km vertically above the plate boundary obtained in this study, for which the ultramafic rock phase diagram was applied. The water content distribution is shown only in the region where the upper surface of the oceanic plate is shallower than the bottom of the model (depth of 200 km) and where the temperature is higher than 200 °C, for which phase diagram data exist. The thin black box represents the horizontal projection of the model region. The two black dashed lines (e) and (f) denote the profiles along the current subduction direction passing through the positions $(x, y)$ = (5 km, 200 km) and (5 km, -510 km), respectively. The black open triangles are volcanoes. The map was created by using Generic Mapping Tools (GMT) (Wessel and Smith 2016). b Same as (a), with the exception of the plate boundary, for which the turbidite phase diagram was applied. c Same as (a), with the exception of the plane at a depth of 4 km from the plate boundary, for which the MORB phase diagram was applied. d Same as (a), with the exception of the plane at a depth of 9 km from the plate boundary, for which the ultramafic rock phase diagram was applied. e Depth distribution of the water content in the slab in the vertical cross-section along the current subduction direction (e) in (a), passing through the position $(x, y)$ = (5 km, 200 km). The two solid black lines indicate the upper and lower surfaces of the oceanic plate. The red solid triangles indicate volcanoes within a one-sided width of 30 km. f Same as (e), with the exception of the profile (f) in (a), passing through the position $(x, y)$ = (5 km, -510 km)

Fig. 5

a Horizontal projection of the current dehydration gradient distribution on the plane 2 km vertically above the plate boundary obtained in this study, for which the ultramafic rock phase diagram was applied. The dehydration gradient distribution is shown only in the region where the upper surface of the oceanic plate is shallower than the bottom of the model (depth of 200 km) and where the temperature is higher than 200 °C, for which phase diagram data exist. The thin black box represents the horizontal projection of the model region. The two black dashed lines (e) and (f) denote the profiles along the current subduction direction passing through the positions $(x, y)$ = (5 km, 200 km) and (5 km, − 510 km), respectively. The black open triangles are volcanoes. The map was created by using Generic Mapping Tools (GMT) (Wessel and Smith 2016). b Same as (a), with the exception of the plate boundary, for which the turbidite phase diagram was applied. c Same as (a), with the exception of the plane at a depth of 4 km from the plate boundary, for which MORB phase diagram was applied. d Same as (a), with the exception of the plane at a depth of 9 km from the plate boundary, for which the ultramafic rock phase diagram was applied. e Depth distribution of the dehydration gradient in the slab in the vertical cross-section along the current subduction direction (e) in (a), passing through the position $(x,y)$ = (5 km, 200 km). The two solid black lines indicate the upper and lower surfaces of the oceanic plate. The red solid triangles indicate volcanoes within a one-sided width of 30 km. The blue histogram indicates the vertical sum of the dehydration gradients at every 10 km in the horizontal distance along the profile. f Same as (e), with the exception of the profile (f) in (a), passing through the position $(x, y)$ = (5 km, − 510 km)

Above the slab, brucite, antigorite, chlorite, and amphibole phases exist in this order from the trench to inland, and the water content varies greatly at each phase boundary (Fig. 4a). Therefore, there are three dehydration zones above the slab with dehydration gradients of approximately 0.16 wt%/km, 0.22 wt%/km, and 0.12 wt%/km from the trench side (Fig. 5a). In the oceanic sedimentary layer, the water contents are 3–4 wt% on the trench side due to the dominance of the lawsonite blueschist phase, but the dehydration reaction in the eclogite phase gradually occurs in inland locations (Fig. 4b). Because of the gradual change in the water content, the distribution of the dehydration gradient is approximately 0.02–0.06 wt%/km (Fig. 5b). In the oceanic crust, the dehydration reactions from the lawsonite blueschist phase to the lawsonite eclogite phase and from the blueschist phase to the amphibole phase rapidly change the water content (Fig. 4c). This change results in a dehydration gradient of approximately 0.10 wt%/km at this location in the dehydration zone (Fig. 5c). The water content and dehydration gradient of the slab mantle, which are calculated using the same hydrous mineral, are distributed like those above the slab (Figs. 4d and 5d). However, the dehydration zones are located slightly inland because the temperature distribution of the slab mantle is colder than that above the slab. The values of the dehydration zones are 0.12 wt%/km, 0.20 wt%/km, and 0.12 wt%/km from the trench side (Fig. 5d).

Discussion

Temperature distribution

Rodriguez Piceda et al. (2022) calculated the 3D thermal distribution in the south central Andes (29–39°S). They calculated the surface heat flow and discovered that the values reached 80–105 mW/m² in the orogenic zone and 40–75 mW/m² in the forearc, which were caused by the shallow cold slab. The model reproduced heat flow in the orogenic zone by lateral changes in the radiogenic continental crust. However, since we have not considered such lateral variations and volcanic heat sources in our numerical simulation, the inland heat flow is almost constant at approximately 80 mW/m², and there is no heat flow higher than that in the surrounding area near the volcanic chain. Conversely, we obtain a relatively low heat flow on the forearc, which is similar to the values in Rodriguez Piceda et al. (2022) because of the effect of the subducting cold oceanic plate. It is possible that their model did not consider the flow velocity, resulting in a different thermal structure from that of our model, which obtaines the temperatures and flow velocities at each time step, by solving a coupled problem on mass, momentum and energy equations.

We compared the temperature distribution of the source regions of the two megathrust earthquakes that occurred at the plate boundary (Fig. 3). The average temperatures along the strike were determined at the updip and downdip limits, with slips larger than 5 m for the coseismic slip distribution of the 1960 Valdivia earthquake (Ho et al. 2019) and that of the 2010 Maule earthquake (Moreno et al. 2012). We only determined the average temperature at the downdip limit of the 1960 Valdivia earthquake because the coseismic slip reached the trench in all inverted slip distributions by Barrientos and Ward (1990), Moreno et al. (2009), and Ho et al. (2019). As a result, the average temperature at the downdip limit of the 1960 Valdivia earthquake was estimated to be 432 ± 57 °C. The average temperatures of the 2010 Maule earthquake were estimated to be 167 ± 52 °C at the updip limit and 289 ± 27 °C at the downdip limit.

Hyndman et al. (1997) suggested that the temperatures at the updip limit and downdip limit of the seismogenic zone at the plate boundary are 100–150 °C and 350 °C, respectively, and that megathrust earthquakes occurring in this temperature range can propagate to the 450 °C isotherm as the slip decreases. Our results are consistent with the viewpoint of the rate-and-state friction law that controls the type of slip occurring at the plate interface (Lay and Kanamori 1981; Rice and Gu 1983). Here, the location of the 350 °C isotherm coincides with the transition from a shallower region of the plate interface dominated by brittle (velocity weakening) frictional behavior to a deeper region dominated by ductile (velocity strengthening) frictional behavior, which has been imaged in the region by several interplate coupling models in the Chile subduction zone (Moreno et al. 2010; Métois et al. 2012). In addition, Valdenegro et al. (2019) compared the 2-D thermal model along 33°S with the location of hypocenters with moment magnitudes greater than 5.0 that occurred between 33 and 34°S along the Chile subduction zone over the past 20 years. Their results showed that most of the earthquakes occurred in the forearc and arc domains within the upper crust located on the colder side of the 350 °C isotherm. For the 2010 Maule earthquake, the average temperature of the downdip limit was consistent with that of previous studies, suggesting that temperature conditions are important for the occurrence of interplate earthquakes in the southern Chile subduction zone and that the slip behavior changes around the 350 °C isotherm. Conversely, for the 1960 Valdivia earthquake, the average temperature of the downdip limit is between 350 °C and 450 °C, which appears to be exceptional. This trend arises because the 1960 (Mw 9.5) Valdivia earthquake was the instrumentally recorded largest earthquake in the world, and a young hot plate is subducting near the Chile Ridge, located at approximately 46°S.

Distributions of the water content and dehydration gradient

Spinelli et al. (2016) constructed 2-D thermal structure models of a vertical cross-section in the southern Chile subduction zone to investigate the location of dehydration from a turbidite and MORB. In their models, similar to the dehydration gradient distribution of the oceanic sedimentary layer obtained in our study, the dehydration from the turbidite in their models was located close to the trench side of the volcanic chain and did not contribute to melt production. In contrast, the oceanic crust was estimated to be dehydrated just below the volcanic chain and contribute to the melt production and the formation of volcanoes. In the dehydration gradient distribution of this study, the dehydration zone of MORB was located slightly trenchward of the volcanic chain south of 42°S (Fig. 4c and f), which is near the profiles of Spinelli et al. (2016). This may be because the temperature range at the slab surface below the volcanic chain was estimated to be 565–775 °C in their model, whereas we estimated a slightly higher temperature range of approximately 825 ± 115 °C.

Fluid generation and its movement are important for understanding igneous activity associated with volcanism in subduction zones. The dehydrated water from the slab rises upward and decreases the melting point of the mantle wedge, contributing to the production of magma. In addition, previous studies that perform seismic waveform analyses and numerical simulations indicate that most of the water released from the subducting slab is captured at the bottom of the mantle wedge above the slab in the shallow part of the subduction zone. Then, the water is transported as serpentinite and chlorite along the top of the slab surface to the deep portion (e.g., Iwamori 1998; Kawakatsu and Watada 2007; Wada et al. 2012, Tastumi et al. 2020). At great depths, serpentinite and chlorite undergo dehydration reactions, and the released water rises vertically upward, resulting in melt production and forming volcanoes.

We compare the locations of the dehydration zones and the distribution of the volcanoes. The volcanic chain is nearly parallel to adjacent dehydration zones in the horizontal projection of all sedimentary layers, oceanic crust, and slab mantle and overlaps some dehydration zones (Fig. 5a–d). Furthermore, the vertical sum of the dehydration gradients in the slab is high below the volcanoes (Fig. 5e and f). Therefore, we propose that the water dehydrated from the slab rises to the mantle wedge, decreases the melting point, and contributes to the production of magma. Assuming that the fluid rises vertically after dehydration, the dehydration zone should be concentrated immediately below the volcanic chain rather than spread out around the volcanic chain. We propose a supplementary process that includes dehydration from the bottom of the mantle wedge above the slab, as in the previous studies described above, and the dehydration process from the inside the slab, as described below.

At the bottom of the mantle wedge immediately above the dehydration zone, the water content is greater than 0 wt% in the forearc (Fig. 4a), and the vertical sum of the dehydration gradients in the slab is high in the forearc (Fig. 5e and f). Therefore, water released from the subducting slab is likely to react with the bottom of the mantle wedge along the plate interface, partially becoming hydrous minerals. Subsequently, the water is dragged to a deep portion by the subduction of the oceanic plate. The fluid above the slab mantle is released again immediately below the volcanic chain, rises to a shallow depth and decreases the melting point to produce magma. However, based on the amount of released water, this process is supplementary, as described above.

In addition, a relationship between dehydration from the slab and the source distribution of regular earthquakes has been suggested in Cascadia, where a young oceanic plate is subducting (e.g., Ji et al. 2017), which is the same as in the southern Chile subduction zone. Therefore, we compare the distribution of the dehydration gradient with the slip distribution of the 1960 Valdivia and the 2010 Maule earthquakes. The dehydration zones of the sedimentary layer and the slab mantle overlap with the distributions of these coseismic slips, suggesting that dehydration from the slab might have influenced the occurrence of megathrust earthquakes in the southern Chile subduction zone. Although the distributions of dehydration zones and coseismic slips overlap, there is little correlation between the region with the maximum value of the dehydration gradient and that with the most coseismic slips. Therefore, further study is needed to determine how much dehydration may have affected the occurrence of these earthquakes.

Conclusions

In this study, we constructed a three-dimensional thermomechanical structure model associated with the subduction of the Nazca plate in the southern Chile subduction zone. From the relationship between temperature and depth obtained from this numerical simulation, we calculated the water content and dehydration gradient near the plate boundary. The significant results obtained in this study are summerized as follows:

• 1. For the 2010 Maule earthquake (Mw8.8), the temperature range for the source region with coseismic slips larger than 5 m was determined to range from 167 ± 52 °C to 289 ± 27 °C, suggesting that most of the coseismic slip occurred on the lower side of the 350 °C isotherm at the plate boundary, dominated by a brittle frictional regime. However, seismic slip did not occur in the ductile regime. In contrast, for the 1960 Valdivia earthquake (Mw9.5), we determined temperature range with coseismic slips larger than 5 m as a maximum of 432 ± 57 °C at the downdip limit of the rupture, which appeared to be exceptional. This result arose because the magnitude of the 1960 Valdivia earthquake was considerably large and the young hot oceanic plate is subducting near the Chile Ridge.

• 2. Below the volcanic chain, the dehydration gradients in the oceanic sedimentary layer, the oceanic crust, the slab mantle, and at the bottom of the mantle wedge above the subducted slab, were approximately 0.04 wt%/km, 0.10 wt%/km, 0.12 wt%/km, and 0.12 wt%/km, respectively. Based on the results of the vertical sum of the former three processes, we suggest that water that dehydrated from the slab rose vertically and formed the volcanic chain. In addition to this effect, we propose the following supplementary process: dehydration from the bottom of the mantle wedge. Fluids released from the subducted slab were captured at the bottom of the mantle wedge above the slab and subducted to depth. The fluids at the bottom of the mantle wedge were released below the volcanic chain, rose upward vertically, and decreased the melting point, contributing to magma production to form volcanoes.