3-D thermal structure and dehydration near the Chile Triple Junction and its relation to slab window, tectonic tremors, and volcanoes

Abstract

The southern Chile subduction zone is a complex tectonic environment, where the Chile Ridge, the Nazca (NZ) and Antarctic (AN) plates subduct underneath the South American (SA) plate. The intersection between the NZ, AN and SA plates is referred to as the Chile Triple Junction (CTJ). In this region, a gap, often referred to as a slab window, has been formed between the NZ and AN slabs due to the divergence in their plate motion velocities, with volcanoes existing mainly above the subducted NZ and AN plates. In this study, we constructed a three-dimensional thermomechanical model associated with simultaneous subduction of the NZ and AN plates near the 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 due to the northward migration of the CTJ and the slower convergence rate of the AN plate. Moreover, we calculated the water content and dehydration gradient from the temperature distribution near the upper surface of the slab and discussed their relationship to the distribution of volcanoes. In the northern part of the model domain, high dehydration gradients were obtained below the volcanic chain. Therefore, we suggest that the water released from the slab and the mantle wedge decreased the melting point of the mantle wedge just above the slab and produced melts, which may have contributed to form the overlying volcanoes.

Introduction

In the southern Chile subduction zone, the Nazca (NZ) and the Antarctic (AN) plates are generated by seafloor spreading along the Chile Ridge. The location where the NZ, the AN, and the South American (SA) plates meet is called the Chile Triple Junction (CTJ). The CTJ is considered to have moved northward from the Miocene to the present and is currently located at approximately 46°S (Fig. 1). Within this subduction zone, a gap called a slab window is proposed to have formed where no subducting plate exists due to the subduction of a spreading ridge between NZ and AN plates (Thorkelson and Taylor 1989) that experiences a divergent motion as they subduct beneath the SA plate. The slab window geometry has been estimated by plate reconstruction and observed by seismic tomography (e.g., Breitsprecher and Thorkelson 2009; Russo et al. 2010; Mark et al. 2022). The spatial distribution of the slab window shows that volcanoes exist mainly above the NZ and AN plates, with only a few located above the slab window. Tectonic tremors have also been observed in this region (Ide 2012; Gallego et al. 2013; Sáez et al. 2019). A tectonic tremor is a type of slow earthquake which occurs near the upper surface of the slab, and dehydration near the slab top surface is considered to contribute to its occurrence (e.g., Obara 2002; Iwamoto et al. 2022).

Fig. 1

Tectonic map of the study area in and around the CTJ. The inset shows the location of the study area as a red boxed area within South America (SA). The blue dashed rectangle delineates the model region and the two blue orthogonal arrows indicate the $+x$ and $+y$ axis directions, respectively. The thick black lines delineate the plate boundaries at the Earth’s surface (Bird 2003), and the thin black lines indicate the depth (in km) of the upper surface of the oceanic plate geometry model used in this numerical simulation. The pink dashed line is the approximate boundary of the slab window at a depth of 50 km (Russo et al. 2010). Yellow lines show the current age of the NZ and AN plates (Seton et al. 2020). The colored circles represent the depths and epicenters of all tectonic tremors that occurred from January 2005 to February 2007 (Ide 2012). The light blue dashed line indicates the tectonic tremor-generating area used in Figs. 5 and 6. The volcanoes are represented by solid red triangles. The purple and orange arrows are the convergence rate vectors of the NZ and AN plates, respectively, with respect to the SA plate, calculated using ITRF2014 (Altamimi et al. 2016). The bathymetry and topography data are taken from ETOPO1 (Amante and Eakins 2009)

Previous 2-D thermal modeling of ridge subduction showed elevated temperatures and a significant water flux in the vicinity of the upper surface of the slab, although the model is idealized and assumes that the NZ and AN plates move together after subduction (Iwamori 2000). Groome and Thorkelson (2009) performed 3-D thermal structure modeling in the region, where the slab window exists to investigate the effects of ridge subduction and slab window migration on the thermal structure. They propose a simplified model that only considers thermal conduction, assuming a constant convergence rate and migration velocity of the slab window. Even with such a simplified model, they obtained a high temperature distribution on the forearc and asymmetric temperature structures on both sides of the slab window due to its structure and northward motion.

In and around the CTJ, the presence of the slab window is considered to affect the thermal structure, being highly relevant for heat flow (Villar-Muñoz et al. 2014; 2021), shallow depth melts (Anma and Orihashi 2013) and Patagonian plateau basalts formed by magma rising through the slab window (Gorring et al. 1997; Guivel et al. 2006). To better understand the thermal structure in this region, it is necessary to properly consider various factors, such as a precise subduction geometry, the subduction history of the oceanic plates, and the formation of the slab window. In this study, we use a more realistic 3-D slab geometry than those of previous studies, and consider the spatiotemporal variation in the ages and convergence rates of the NZ and AN plates along the trench axis. Under these conditions, we performed a three-dimensional thermomechanical numerical simulations to predict the thermal structure associated with the simultaneous subduction of these plates including the subduction of the active spreading center. Furthermore, using the obtained thermal structure and phase diagrams of hydrous minerals, we calculated the distribution of water content and the dehydration gradient near the upper surface of the NZ and AN slabs, and investigated the relationship between these distributions and the locations of tectonic tremors and volcanoes.

Results

Comparison between the observed and calculated Curie point depths

Manea and Manea (2011) studied the crustal thermal structure in Mexico, by estimating the distribution of the Curie point depth based on magnetic anomaly data. They found the Curie point depth to be relatively deep in the forearc region, consistent with the position of the flat subducted Cocos plate, and concluded that the Curie point depth varies with the geometry of the subducting plate. In addition, Ji et al. (2019) performed numerical simulations of a 3-D thermal structure, similar to our study, in north-central Chile, and calculated the heat flow from their thermal structure model. Their model is constrained by both the observed heat flow, and also the heat flow converted from the Curie point depth in Li et al. (2017). The Curie point is the temperature at which a material loses its magnetism, and Li et al. (2017) performed the inversion of the magnetic anomaly data to estimate the Curie point depth, defining this depth as the 550 °C isothermal surface based on Mayhew (1982) and Tselentis (1991). In this study, we consider the depth of the 550 °C isotherm of the thermal structure obtained in our numerical simulation as the Curie point depth. We compare such an isotherm spatial distribution with the Curie point depth distribution of Li et al. (2017) to evaluate the validity of the 3-D thermal structure obtained in our numerical simulation. The scarcity of heat flow observations in the region defined by the model domain limits our ability to statistically validate our models against an adequate/representative number of measurements. Instead, we define our optimal model as the one that minimizes the RMS of the residuals between the Curie point depth distribution of our numerical simulation and that of Li et al. (2017), see Eq. (S1) in Text S1. We refer the reader to the Supplementary Information for details of our methodology and the model used to calculate the 3-D thermal structure associated with subduction of the oceanic NZ and AN plates.

Figure 2 shows comparison between the distribution of Curie point depths of Li et al. (2017) and that for the optimal model. The Curie point depth distribution obtained in this study is shallower closer to the ridge due to the young oceanic plates subducting near the trench (Fig. 2b, f). In addition, the calculated Curie point depths were deeper in the northern part of the model domain than those in the southern part around the profile G–G’ in Fig. 2b. This is due to the northern part of the model domain being far from the northward-moving ridge, and the fast convergence rate of the NZ plate. It should be noted that the  calculated Curie point depth along the profile G–G’ is not deep, because this area is strongly affected by frictional heating at the plate boundary (Fig. 2g). On the inland side, the calculated Curie point depths are nearly constant, because the shallow part of the model domain is covered by the stratified upper and lower crusts (Fig. 2h). The Curie point depths calculated in this study show a greater variability than those of Li et al. (2017), especially in terms of the extreme values in the region. However, it should be noted that the positions of these extreme values correspond to seismic velocity anomalies identified in the region, as discussed in the following section.

Fig. 2

a Colors represent the spatial distribution of the Curie point depth from Li et al. (2017), in positions where data exist in the model domain. The black dashed lines are the locations of the profiles shown in (c–h). b Colors represent the spatial distribution of the Curie point depth for the optimal model obtained in this study. Profiles are the same as in (a). c Comparison between observed and calculated Curie point depths along profile C–C’. The black and red lines denote values of Curie point depths in Li et al. (2017) and calculated values obtained from this study, respectively. d–h are the same as (c) except for profiles D–D’ through H–H’, respectively

Thermal structure

Because the grid interval in the depth direction is approximately 3.3 km in this numerical simulation, we performed a linear interpolation of temperature between each grid point in the depth direction and obtained the temperature distribution along the slab surface in the optimal model (Fig. 3). It should be noted that the slab length of AN plate is shorter than that of the NZ plate, because the AN plate has a slower convergence rate meaning the slab did not penetrate across the width of the model domain during the calculation time of 15 Myr. We obtained a tendency for the temperature distribution at the upper surface of the slab to strongly depend on the age of the subducting oceanic plate. Therefore, higher temperatures at the slab surface occurred closer to the ridge, even at the same depth. In addition, the subducted current AN plate region had a higher temperature distribution than the subducted current NZ plate region. This can be attributed to the following two reasons: the CTJ, which is the hottest point along the trench axis, was located at $y=-150\ \mathrm{km}$ at the beginning of the calculation (15 Ma), and subsequently moved northward to its current position of y = 0 km. Therefore, from the beginning of the calculation to the present, the temperature at the slab surface in the current AN plate region is higher than that of the current NZ plate region, because the younger oceanic plate is subducting in the current AN plate region. Another possible cause is that the AN plate is subducted at a slower rate than the NZ plate, so that the subducting plate is heated by the surrounding mantle for a longer period.

Fig. 3

Horizontal projection of the temperature distribution at the upper surface of the slab for the optimal model at present (0 Ma). The contour interval is 50 °C. The thin black box is the model region, and the thick black lines are the plate boundaries at the Earth’s surface (Bird 2003). The pink dashed line is the approximate boundary of the slab window at a depth of 50 km (Russo et al. 2010). The white open circles represent the epicenters of all tectonic tremors that occurred from January 2005 to February 2007 (Ide 2012). The black open triangles are volcanoes

Distributions of water content and dehydration gradient

As previously described, we performed a linear interpolation of temperature in the depth direction, using the thermal structure of the optimal model, to calculate temperature values, at 1 km intervals in the depth direction, throughout the slab surface. We calculated the maximum water content, corresponding to the temperature and depth of each data point, to obtain the distributions of the maximum water content and the dehydration gradient along the slab surface. The slab was assumed to consist of three layers, and different phase diagrams were applied to each layer to determine the water content and dehydration gradient of hydrous minerals near the upper surface of the slab (Fig. 4). The dehydration gradient is defined as the difference in the water content (wt%) per unit length (km) along the subduction direction (Suenaga et al. 2019). We use the phase diagrams of ultramafic rock, turbidite, and MORB  with the following respective layers thickness; the mantle wedge was set up to extend 2 km vertically above the upper surface of the slab, the oceanic sedimentary layer was set from the upper surface of the slab to a depth of 2 km, the oceanic crust was set at depths between 2 and 7 km from the upper surface of the slab, and the slab mantle was set from a depth of 7 km from the upper surface of the slab to the bottom of the oceanic plate (Fig. S1).

Fig. 4

a Phase diagram of ultramafic rock (Tatsumi et al. 2020). The purple and red lines show the p–T paths 2 km above the upper surface of the slab in the vertical cross sections along the profiles E–E’ and F–F’, respectively, in Figs. 5 and 6. b Phase diagram of turbidite (van Keken et al. 2011). The purple and red lines show the p–T paths at the upper surface of the slab in the vertical cross sections along the profiles E–E’ and F–F’, respectively, shown in Figs. 5 and 6. c Phase diagram of MORB (Tatsumi et al. 2020). The purple and red lines show the p–T paths at a depth of 4 km from the upper surface of the slab in the vertical cross sections along the profiles E–E’ and F–F’, respectively, in Figs. 5 and 6. BS: blueschist, LwsBS: lawsonite blueschist, AMP: amphibolite, GR granulite, AmpEC: amphibole eclogite, ZoEC: zoisite eclogite, LwsEC: lawsonite eclogite, DryEC: dry eclogite

The water content distribution at present (0 Ma) is shown in Fig. 5, and the distribution of the dehydration gradient along the current subduction direction is shown in Fig. 6. In all layers, the profiles of water content along the trench axis had lower values near the ridge due to its higher temperature (Fig. 5a–d). At the bottom of the mantle wedge just above the slab, brucite, antigorite, chlorite, and amphibole phases exist, in this order, from the trench to the inland side, and the water content varies greatly at each phase boundary (Figs. 4a and 5a). Therefore, from the trench to the inland side, there are three dehydration zones at the bottom of the mantle wedge, just above the slab, with dehydration gradients of approximately 0.16 wt%/km, 0.22 wt%/km, and 0.12 wt%/km (Fig. 6a). In the oceanic sedimentary layer, the dehydration gradient distribution has small values, with a maximum of 0.04–0.08 wt%/km, due to the gradual change from phengite lawsonite blueschist phase to amphibole phengite zoisite eclogite phase along the subduction direction (Figs. 4b, 5b and 6b). In the oceanic crust, the dehydration reaction from blueschist phase to amphibole phase and amphibole eclogite phase occurs in the southern to central part of the model domain. In contrast, the dehydration reaction from the lawsonite blueschist phase to the lawsonite eclogite phase occurs on the inland side in the northern region (Figs. 4c and 5c). These two different dehydration reactions occur in the oceanic crust, the dehydration zone with 0.1 wt%/km is located near the trench in the southern to central part, whereas the dehydration zone with 0.06 wt%/km is located inland in the northern part (Fig. 6c). The distribution of the dehydration gradient at a depth of 9 km from the upper surface of the slab, which is in the slab mantle layer, is located inside the slab, not directly exposed to the hot surrounding mantle, and so it can maintain lower temperatures on the inland side. Therefore, the dehydration reactions from antigorite phase to chlorite phase and from chlorite phase to amphibole phase have moved farther inland, rather than above the slab, which uses the same phase diagram (Figs. 5d and 6d).

Fig. 5

a Horizontal projection of the water content distribution on the plane 2 km vertically above the upper surface of the slab at present (0 Ma) obtained in this study. The water content distribution is shown only in the region, where the temperature is higher than 200 °C, for which phase diagram values exist. The thin black box is the model region. The two black dashed lines are the profiles E–E’ and F–F’ shown in (e) and (f), respectively. The black open triangles are volcanoes. The encircled green solid line indicates the tectonic tremor-generating area. b Same as (a) except for the upper surface of the slab. c Same as (a) except for the plane at a depth of 4 km from the upper surface of the slab. d Same as a except for the plane at a depth of 9 km from the upper surface of the slab. e Distribution of water content in the vertical cross section of the slab along the current subduction direction E–E' of the NZ plate with respect to the SA plate. The two solid black lines indicate the upper and lower surfaces of the oceanic plate. The green open circles indicate tectonic tremors with residual travel time errors of less than 0.5 s from Ide (2012). Tectonic tremors are plotted only within a width of 10 km (one-sided width of 5 km) along the profile E–E’. The red solid triangles indicate volcanoes within a width of 50 km (one-sided width of 25 km) along the profile E–E’. f Same as e except for the profile along the current subduction direction F–F' of the AN plate with respect to the SA plate

Fig. 6

a Horizontal projection of the dehydration gradient distribution on the plane 2 km vertically above the upper surface of the slab at present (0 Ma) obtained in this study. The dehydration gradient distribution is shown only in the region, where the temperature is higher than 200 °C, for which phase diagram values exist. The thin black box is the model region. The two black dashed lines are the profiles E–E’ and F–F’ shown in e and f, respectively. The black open triangles are volcanoes. The encircled green solid line indicates the tectonic tremor-generating area. b Same as a except for the upper surface of the slab. c Same as a except for the plane at a depth of 4 km from the upper surface of the slab. d Same as a except for the plane at a depth of 9 km from the upper surface of the slab. e Distribution of dehydration gradient in the vertical cross section of the slab along the current subduction direction E–E’ of the NZ plate with respect to the SA plate. The two solid black lines indicate the upper and lower surfaces of the oceanic plate. The green open circles indicate tectonic tremors with residual travel time errors of less than 0.5 s from Ide (2012). Tectonic tremors are plotted only within a width of 10 km (one-sided width of 5 km) along the profile E–E’. The red solid triangles indicate volcanoes within a width of 50 km (one-sided width of 25 km) along the profile E–E’. f Same as e except for the profile along the current subduction direction F–F' of the AN plate with respect to the SA plate

Discussion

Thermal structure

The 3-D model presented in this study is an idealized model with a constant northwards migration velocity of the CTJ. The chosen value of 1.0 cm/year is more-closely aligned to the recent tectonic situation, and, therefore, is most representative of the regions, where the spreading center has not yet been subducted (to the north of the present-day CTJ) and where the ridge has recently been subducted (in the vicinity of the present-day CTJ). These regions coincide with the positions of observed tectonic tremors (Ide 2012), young NZ plate subduction (Russo et al. 2010), recent ridge subduction (Gallego et al. 2010) and the southernmost part of the Southern Volcanic Zone in Chile (Stern 2004). The locations of older slab windows, located southwards of the present-day CTJ, will be positioned slightly to the east of their actual positions in this model, due to the assumption of AN plate subduction starting in the model at 15 Ma, and, therefore, making further eastwards progress than that calculated by Breitsprecher and Thorkelson (2009). However, the thermal conditions on the upper surface of the slab for these older slab windows will be largely representative of the actual situation.

We compared the temperature distribution at the upper surface of the slab obtained in this study with the boundary of the slab window at a depth of 50 km in Russo et al. (2010) (Fig. 3). It should be noted that the horizontal distance from the trench to the location where the depth of the upper surface of the slab is 50 km is approximately x = 170 km with the slab geometry used in this study. The isotherms at the upper surface of the slab obtained in this study align approximately with the boundary of the slab window inferred by seismic tomography. The calculated temperature inland from the slab window boundary exceeds 1000 °C (Fig. 3). In addition, the seismic tomography of Russo et al. (2010), at depths of 100 and 200 km showed a fast Vp anomaly associated with the subduction of the NZ plate to the north of 46°S. This is consistent with the results of this study which shows lower temperatures north of the CTJ and an increased Curie point depth in the position of this fast tomographic anomaly. A second seismic tomography study of Gallego et al. (2010) observes a slow Rayleigh wave group velocity anomaly corresponding to the most-recently subducted ridge underneath the Taitao Peninsula (Fig. 1), which is consistent with the westward extension of the higher (> 800 °C) temperatures observed on the upper surface of the slab and shallowest Curie point depths obtained in this study. Therefore, the thermal structure of the upper surface of the slab is consistent with the locations of the slab window and the oceanic plates observed in previous studies. Although we do not account for the materiality of the slab window in the used continuous 3-D slab geometry, our analysis recognizes the slab window as a thermal boundary.

Distributions of water content and dehydration gradient

Hydrous minerals, transported by the subduction of the oceanic plate, produce water through dehydration reactions, which decrease the effective normal stress near the upper surface of the slab and generate tectonic tremors (e.g., Obara 2002). It is also shown that dehydration from the hydrous minerals decreases the melting point of the mantle wedge to produce melts and contribute to the formation of volcanoes (e.g., Iwamori 1998; Kawakatsu and Watada 2007). Therefore, we compared the distributions of water content and the dehydration gradient obtained in this study with the source region of tectonic tremors and the distribution of volcanoes (Figs. 5 and 6).

At the bottom of the mantle wedge, just above the slab, no high dehydration gradient was observed within the tectonic tremor-generating area: the dehydration from antigorite to chlorite was approximately 0.22 wt%/km on the updip side of the tremor area, and the dehydration from chlorite to amphibole was approximately 0.12 wt%/km on the downdip side (Fig. 6a). In the oceanic sedimentary layer and crust, the water content gradually decreases going from the trench to the updip limit of the tectonic tremor-generating area and reaches nearly 0 wt% in the tremor area (Fig. 5b, c). Thus, the high dehydration gradient values are located near the trench, on the updip side of the tectonic tremor-generating area (Fig. 6b, c). In the slab mantle, the maximum water content is about 2–4 wt% near the tectonic tremor-generating area, but a high dehydration gradient was not identified near the tremor area due to the relationship between the maximum water content distribution and the subduction direction (Figs. 5d and 6d). Because of these results, it is difficult to explain the occurrence of tectonic tremors by the dehydration gradient distribution of this optimal model, as previous studies have shown that dehydration contributes to the occurrence of tectonic tremors in southwest Japan and southern Alaska (e.g., Obara 2002; Iwamoto et al. 2022). This is because the tectonic tremor-generating area is near the CTJ, where the upper surface of the slab experiences significantly elevated temperatures, and consequently the water is dehydrated before it can be transported to the tremor area. Therefore, to associate the dehydration with the occurrence of tectonic tremors, it is necessary to shift the location of dehydration in the downdip direction by decreasing the temperature of the slab. Incidentally, when the frictional heating at the plate boundary was not considered and the northward ridge migration rate was 4.0 cm/year, the temperature at the upper surface of the plate boundary was decreased. Thus, the hydrous minerals can exist in the tectonic tremor-generating area and dehydration was observed in the mantle wedge just above the slab, oceanic sedimentary layers, and oceanic crust near the tremor area (Additional file 1: Figs. S3 and S4). Therefore, we believe that a parameter setting that results in lower temperatures near the upper surface of the slab when compared to the optimal model is preferable to explain the relationship between dehydration and the tectonic tremor-generating area.

In the northern part of the model domain, where volcanoes are densely aligned in the north–south direction, the dehydration gradients along the current subduction direction are large near the volcanic chain in each layer (Fig. 6). Specifically, there are two dehydration zones of 0.22 wt%/km and 0.12 wt%/km in the mantle wedge just above the slab. In addition, approximately 0.03 wt%/km in the oceanic sedimentary layer, 0.06 wt%/km in the oceanic crust, and 0.22 wt%/km in the slab mantle. Therefore, a considerable amount of water was supplied near the upper surface of the slab and in the bottom of the mantle wedge just above the slab, decreasing the melting point and producing magma to form volcanoes above. In addition, it is considered that little water is transported to the downdip side of the tectonic tremor-generating area, cutting off the Southern Volcanic Zone, because the tectonic tremor-generating area is located closer to the ridge and has a higher temperature distribution compared to the northern part of the model domain.

Conclusions

In this study, we constructed a three-dimensional thermomechanical model, in and around the CTJ, associated with subduction of the NZ and AN plates, underneath the SA plate. From the relationship between temperature and depth obtained in our numerical simulation, we calculated the water content and dehydration gradient near the upper surface of the slab. The significant results obtained in this study are summarized as follows:

• A) The temperature distribution at the upper surface of the slab is higher closer to the ridge, and the current subducted AN plate has a higher temperature than the current subducted NZ plate. This is because the CTJ has moved northward to reach its current position from the beginning of the calculation (15 Ma), and the AN plate is subducting at a slower convergence rate than the NZ plate.

• B) At the mantle wedge just above the slab, no high dehydration gradient was observed within the tectonic tremor-generating area. The dehydration from antigorite to chlorite was approximately 0.22 wt%/km on the updip side, and the dehydration from chlorite to amphibole was approximately 0.12 wt%/km on the downdip side of the tremor area. The water content in the oceanic sedimentary layer and oceanic crust gradually decreased on the updip side of the tectonic tremor-generating area, because the ridge is subducting and the temperature is high near the current CTJ. Therefore, it is difficult to explain the relationship between the occurrence of tectonic tremors and dehydration with the optimal model in this study. To explain this relationship, the model with lower slab surface temparature with no frictional heating at the plate boundary and high ridge migration rate is preferred.

• C) In the northern part of the model domain, the relatively low temperature distribution allows the transportation of hydrous minerals to the deeper inland side. The maximum dehydration gradients are 0.12–0.22 wt%/km in the mantle wedge just above the slab, 0.03 wt%/km in the oceanic sedimentary layer, 0.06 wt%/km in the oceanic crust, and 0.22 wt%/km in the slab mantle. We propose that water released from these layers decreased the melting point, producing a melt which forms the volcanoes above.