In this paper, numerical models were built by the 2D finite element method “ASPECT” (Kronbichler et al., 2012), which solves the three conservation equations of, respectively, momentum, mass, and heat.

1) Stokes equation:

ASPECT regards the mantle as a high-viscosity fluid. The common rheological criteria used in geodynamics are adopted in this paper (Billen and Gurnis, 2001; Lallemand et al., 2005; Jadamec et al., 2013; Holt et al., 2015; Leng and Gurnis, 2015). The viscosity for dislocation or diffusion creep is defined as: η = 1 2 A − 1 n d m n ε ˙ l 1 − n n ⁡ exp ( E + P V n R T ) , (4) where η is the viscosity, A is the prefactor, n is the stress exponent, ε l ˙ is the second invariant of the deviatoric strain rate tensor, d is the grain size, m is the grain size exponent, E is the activation energy, V is the activation volume, R is the gas exponent, and T is the absolute temperature.

In case of diffusion creep, m ≠ 0 and n = 1, the viscosity coefficient in 4) is η diff . While dislocation creep m = 0 and n > 1, the viscosity coefficient in 4) is η disl . So the contributions of diffusion and dislocation creep to the effective viscosity are harmonically averaged into a composite viscosity: η comp = ( 1 η diff + 1 η disl ) − 1 , (5)

Thus, the viscous stress is defined as: σ visco = 2 × η diff , disl , comp  × ε l ˙ , (6) where η diff ,  disl , comp  is one of η diff , η disl , and η comp according to different rheologies (Glerum et al., 2018; Liu et al., 2019).

The ductile rheology is combined with a brittle/plastic rheology to yield an effective visco-plastic rheology. In our implementation, it is defined by the Drucker–Prager criterion (Davis and Selvadurai, 2002): η Drucker − Prager = { η diff , disl , comp  , σ visco < σ yield σ yield 2 ε l ˙ , σ visco ≥ σ yield , (7) σ yield is the yield value: σ yield = C ⁡ cos ⁡ ϕ + P ⁡ sin ⁡ ϕ , (8) where ϕ is the internal friction angle, P is the pressure, and C is the cohesion.

The final effective viscosity η eff is capped by the defined minimum viscosity η m i n (10²⁰ Pas) and maximum viscosity η m a x (10²⁵ Pas) to avoid extremely low or high viscosity values (Glerum et al., 2018): η eff = min ( max ( η Drucker − Prager , η m i n ) , η m a x ) . (9)

Based on the geological background of the CTJ, we built a series of 2D geodynamic numerical models, which fall into Antarctic models and Nazca models, respectively, corresponding to wedge and slab subduction.

Model 1, which simulates the wedge subduction of the Antarctic Plate, corresponds to profile 1 with the convergence rate set as 1.8 cm/year (Figure 1) (Cande and Leslie, 1986; Maksymowicz et al., 2012). The other four models, i.e., Models 2–5, were built to investigate the effect of the slab age and wedge subduction on the lithospheric thermal structure with a set convergence rate of 6.6 cm/year. Models 2 and 3 correspond to wedge subduction, and Models 4 and 5 represent slab subduction. The slab of Model 3 is older than that of Model 2, and the slab of Model 5 is older than that of Model 4. When the calculation finishes, the slab age of Model 4 is the same as that of Model 2 and the slab age of Model 5 is the same as that of Model 3. The slab ages of Models 1–5 are shown in Table 1. Model 4 is a simplification based on profile 2.

The computational domain is 2,500 × 660 km; the continental and oceanic plates are 1,500- and 1,000-km long, respectively. To our knowledge, there is a lack of direct observation of the slab dip of Antarctic plate, which is only inferred to be smaller than that of the Nazca Plate (Breitsprecher and Thorkelson, 2009; Agurto-Detzel et al., 2014). Since there are rarely earthquakes south of CTJ, there have been many different tomographic models depending on the various cases in reality. These models differ greatly for this area. On the basis of the tomography (Figure 2) and previous numerical model settings for the CTJ (van der Hilst and de Hoop, 2005; Scherwath et al., 2006, 2009; Simmons et al., 2010; Völker et al., 2011; Maksymowicz et al., 2012) and the comparison of different initial slab dips ( Surface Heat Flow and Seismogenic Zone and Seismicity Data ), the initial slab dips of Models 1–5 were set in a way in which it gradually increases with the depth from 0° at the trench to 30° at about 120-km deep.

The initial temperature field and the boundary conditions were set as follows: The upper boundary was maintained at a constant temperature of 0°C. The temperature at the lithosphere–asthenosphere boundary was set as 1,300°C (Tetreault and Buiter, 2012). An initial geothermal gradient of 0.25°C/km was set for the asthenosphere; thus, the temperature at the lower boundary was set as 1,440°C (Tetreault and Buiter, 2012). A plate cooling model was used for the oceanic plate with a maximum plate thickness of 95 km. The thickness and initial temperature field of the oceanic lithosphere were calculated according to the plate age (Stein and Stein 1992; Richards et al., 2018): T ( x , y ) = T m [ y a + ∑ n = 1 ∞ 2 n π exp ( − β n x a ) sin ( n π y a ) ] , (10) where a is the asymptotic thermal plate thickness, T _m is the basal temperature, β n = R 2 + n 2 π 2 − R , R = v a / ( 2 κ ) , κ is the thermal diffusivity, and v is the half-spreading rate. x and y denote the horizontal and vertical coordinates. The lithosphere thickness of the continental plate was set as 100 km, and the thickness of the upper and lower crust was both 20 km (Assumpção et al., 2013; Bagherbandi et al., 2017; Liu et al., 2019). A linear interpolation method was used for the initial temperature field of the continental plate (Li and Shi, 2016). The thickness of each layer is listed in Table 2 .

A constant velocity was applied to the upper part of the left boundary to simulate the convergence of the plates and also applied to the lower part in the opposite direction (Liu et al., 2019). The upper boundary was a free surface and the right and lower boundaries were set to be free slip (Figure 3). The parameters obtained in relevant studies are listed in Table 2 (Tetreault and Buiter, 2012; Liu et al., 2019). For the radioactive heat production, this study refers to Xu et al. (2019) and, for the viscosity values, Glerum et al. (2018) and Liu et al. (2019).

According to the initial slab length, a calculation of 5.5 Ma was needed to obtain the initial temperature field of the slab. The calculation process is as follows: First, the boundary velocity was set as 0, and a calculation of 5.5 Ma was carried out. Subsequently, a boundary velocity corresponding to the actual convergence rate was applied to the left boundary, and another calculation of 5.5 Ma was done. In total, the model was run for 11 Ma. After that, the plate ages of Models 1 and 4 match those of profiles 1 and 2, respectively, in Figure 1. The slab ages of the models and the corresponding relationship between the models and profiles are shown in Table 1.

The Nazca model in which the old part subducted first corresponded to slab subduction, whereas the Antarctic model in which the young part subducted first corresponded to wedge subduction. Typically, the slab dip gradually increases with the depth and stops growing at about 80–150 km deep. Beneath this depth, it remains almost constant down to the boundary between the upper and lower mantle (Lallemand et al., 2005). According to the definition of the shallow dip and the deep dip of Lallemand et al. (2005), the mean shallow dips and mean deep dips of Models 2, 3, 4, and 5 are given in Figure 4.

The solid blue lines in Figures 4 and 5 are the sediment contours. Figure 5 shows the evolution of the sediment outlines. The upper parts of the outlines are the upper surface of the slabs. As is shown in Figure 4, the older the slab, the greater the slab dip. In the Antarctic models, for example, the slab of Model 2 is younger than that of Model 3 and the mean shallow dip and mean deep dip of Model 2 are about 2° and 26° smaller than those of Model 3 (Figure 4A, B). In the Nazca models, the slab of Model 4 is younger than that of Model 5, and the mean shallow dip and mean deep dip of Model 4 are around 1° and 4° smaller than those of Model 5 (Figure 4C, D).

When the slab ages of the wedge subduction and the slab subduction are the same, the slab dip of wedge subduction is smaller than that of slab subduction. For example, the mean shallow dip and mean deep dip of Model 2 (Antarctic model, Figure 4A) are, respectively, about 3° and 22° smaller than that of Model 4 (Nazca model, Figure 4C). The slabs of Model 3 and Model 5 are too old for their thickness to increase with age, so they have similar slap dips (Figure 4B, D). Among Models 2, 3, 4 and 5, Model 2 has the smallest slap dip. Based on the results above, the slab dip differences of the areas north and south of the CTJ are considered a consequence of the wedge subduction.

In the effort to obtain a more reasonable initial slab dip, an initial slab dip of 15° was applied to Models 6 and 7. Except for the initial slab dip, the other settings of Models 6 and 7 were the same as those of Models 1 and 4. Figure 6 shows the observed surface heat flow (Lucazeau, 2019) and computed surface heat flow of the different models. In the Antarctic models, the observed heat flow data is too limited to constrain the initial slab dip (Figure 6A). In the Nazca models, the heat flow of Model 4 with an initial slab dip of 30°is more consistent with observed data in the vicinity of the trench than that of Model 7 (Figure 6B).

The larger the convergence rate, the higher the maximum heat flow in the vicinity of the trench. The convergence rate of Model 1 is smaller than that of Model 2, and the maximum heat flow of Model 1 is 133.7°mW/m² lower than that of Model 2 in the oceanic plate and 23.7°mW/m² lower in the forearc (Figure 6C).

The surface heat flow, including the oceanic surface heat flow and forearc heat flow, is generally high in the vicinity of the trench (Figure 6). The surface heat flow in the forearc is complex and depends on several factors such as the topography of the overlying plate or the different slab dips. As is shown in Figure 6D, the forearc heat flow of Model 2 is high, and the heat flow values of Models 3, 4, and 5 gradually stabilize and are almost the same far away from the trench.

Generally, the younger the plate is, the higher the maximum oceanic surface heat flow is in the vicinity of the trench (Figure 6D). In the Antarctic models, the oceanic plate of Model 2 is younger than that of Model 3, and the maximum oceanic heat flow of Model 2 (300.8 mW/m²) is higher. The maximum forearc heat flow of Model 2 (70.8mW/m²) is higher than that of Model 3 (59.2°mW/m²; Figure 6D), which is related to the young slab and small slab dip. In the Nazca models, the oceanic plate of Model 4 is younger than that of Model 5, and the maximum oceanic heat flow of Model 4 is higher (427.0 mW/m²) in the vicinity of the trench, which is related to the young plate and adjacent ridge. However, the maximum forearc heat flow of Model 5 (94.7 mW/m²) is higher than that of Model 4 (73.2 mW/m²).

The position of the interplate seismogenic zone is related to the interplate seismicity and primarily controlled by temperature (Tichelaar and Ruff, 1991; Oleskevich et al., 1999; Völker et al., 2011). The earthquake distribution in Chile (Tichelaar and Ruff, 1991; Oleskevich et al., 1999; Völker et al., 2011) and in the Nazca subduction zone is consistent with this hypothesis. Klotz et al. (2006) generally defines 1) the interseismic locked zone as the region in which the plate contact is at temperatures from 100°C to 350°C, 2) the temperature range from 350°C to 450°C as the transition zone from full locking to full slip, and 3) the seismogenic zone as the width of the locked zone plus half the width of the transition zone (Völker et al., 2011).

The seismicity data (1906–2021) in the vicinity of the profiles from the ISC (International Seismological Centre (2021), On-line Bulletin, https://doi.org/10.31905/D808B830) (Bondár and Storchak, 2011; Storchak et al., 2017; 2020) were projected to profile 1 and profile 2 to contrast with seismogenic zones of our models (Figure 7).

The seismicity data of profile 1 is too little to constrain the seismogenic zone south of the CTJ. In the north of the CTJ, the seismicity data of profile 2 are more consistent with Model 4 than with Model 7 in the seismogenic zone. However, the observed deep seismic events within the slab of profile 2 are more consistent with Model 7 than with Model 4. This is due to the shallow and deep dip of the slab increase during the subduction in Model 4. In the Slab dip section, it is found that the different slab dips of the areas north and south of the CTJ are due to the wedge subduction. Therefore, the initial slab dip of profile 1 and profile 2 were set to be the same. Based on the comparison of the seismogenic zone and heat flow, the numerical models with an initial slab dip that gradually increases from 0° to 30° (Models 1–5) better reflect the interplate seismogenic zone and surface heat flow than those with a dip of 15° (Models 6–7). The widths of the seismogenic zones of Models 1–5 according to the definitions of Klotz et al. (2006) are shown in Table 1.

When the slab subduction and wedge subduction are of the same, old slab age and convergence rate, the width of the seismogenic zone of the former is smaller than that of the latter. As Table 1 shows, the width of the seismogenic zone of Model 5 is 3 km narrower than that of Model 3. The opposite results were obtained for slabs with the same young age and convergence rate, as is illustrated by the 4-km-wider seismogenic zone in Model 4 compared with Model 2. Generally, the widths of the locked zone, transition zone, and seismogenic zone of Model 2 are smaller than those of the other models when the convergence rate is 6.6 cm/year (Table 1).

The older the slab, the wider the seismogenic zone. For instance, the slab of Model 5 is older than that of Model 4, and the seismogenic zone of Model 5 is 2 km wider than that of Model 4. The slab of Model 3 is older than that of Model 2, and the seismogenic zone of Model 3 is 9 km wider than that of Model 2.

The seismogenic zone related to wedge subduction widens with the increase in the convergence rate. The convergence rate of Model 2 is larger than that of Model 1, and the seismogenic zone of Model 2 is 11 km wider than that of Model 1 (Table 1). Besides, as is shown in Table 1, the intersections of the slab surface with the 350°C and 450°C isotherms of Model 2 are much deeper than those of Model 1. The depth of the downdip limit of the seismogenic zone of Model 2 is about 10 km deeper than that of Model 1.

Figure 8 shows the temperature field within 350 km ahead of the trench. Affected by the subduction of the cold slab, the temperature in the vicinity of the upper surface of the slab is low. In contrast, the overlying plate far away from the slab is only slightly affected. The x- and y-coordinates at the intersections of the upper surface of the slab with the isotherms of 100°C, 150°C, 350°C, and 450°C (Figure 8) are shown in Table 1.

At the same age, the slab of wedge subduction tends to be hotter than that of slab subduction horizontally. As Table 1 shows, the slab ages of Models 4 and 2 are the same, and the x-coordinates of Model 4 at intersections of the slab surface with the 100, 350, and 450°C isotherms are larger than those of Model 2. Compared with the differences between their x-coordinates, the differences between the y-coordinates of Models 2 and 4 are much smaller.

Generally, the older the slab is, the lower the temperature is in the vicinity of the upper surface of the slab in the horizontal direction. In the Antarctic models, for example, the slab of Model 2 is younger than that of Model 3. The x-coordinates of Model 2 at the intersections of the slab surface with the isotherms of 100, 150, 350, and 450 °C are smaller than those of Model 3. For Models 2–5, the difference in the y-coordinates is much smaller than the difference in the x-coordinates (Table 1). For example, the maximum difference in the x-coordinates between Models 4 and 5 is 13 km at the 150°C isotherm, whereas the maximum difference in the y-coordinates is 3 km at the 450°C isotherm.

The temperature in the vicinity of the slab surface lowers with a decrease in the convergence rate horizontally. In the Antarctic models, the convergence rate of Model 1 is smaller than that of Model 2. The x-coordinates of Model one at the intersections of the slab surface with the isotherms of 100, 150, 350, and 450°C are larger than or equal to those of Model 2 (Table 1). In addition, the y-coordinates of Model 1 at the intersections of the slab boundary with the 100°C and 150°C isotherms are 8 and 3 km, respectively, deeper than those of Model 2, whereas at the intersections of the slab boundary with the 350 and 450°C isotherms, the y-coordinates of Model 1 are, respectively, 7 and 12 km shallower than those of Model 2. This means that the convergence rate has a more significant effect on the temperature field than the slab age in the vertical direction.

Our simulations show that the slab dip of the Antarctic Plate, which increases southward from the CTJ, is smaller than that of the Nazca Plate when the same slab age is chosen. To the south of the CTJ, the slab between the CTJ and profile 1 is younger than that of profile 1 and, thus, its slab dip is smaller than that of profile 1. It is also smaller than the dip of the Nazca slab (profile 2) because the slab dip of profile 2 is almost the smallest to the north of the CTJ ( Slab Dip ). This is consistent with the tomography of the area south of the CTJ in the vicinity of 47°30S, which shows a smaller slab dip compared with the Nazca Plate (Maksymowicz et al., 2012). However, there is a lack of evidence farther south, But the results are consistent with the inference of Breitsprecher and Thorkelson (2009) based on kinematics. In addition, the slab age is generally considered to be the main factor in the change of the slab dip. The older the oceanic plate age, the greater the slab dip (Jarrard, 1986; Stern, 2002). Tomographic evidence, the tomography in the vicinities of 45°30′S, 44°30′S (Maksymowicz et al., 2012), 43°S (Scherwath et al., 2006; 2009), 38°S (Contreras-Reyes, et al., 2008), and farther north (van der Hilst and de Hoop, 2005), shows that the slab dip increases northward from the CTJ. Our results show that, with the slab age, the slab dip increases when the slabs are young but remains almost unchanged when the slabs are old. The results in Slab Dip show that the slab age and slab dip increase northward from profile 2, which is also consistent with the tomographic evidence above.

If we adopt the definitions of Klotz et al. (2006), the simulation results of Models 2 and 3 show that the temperature of the Nazca slab decreases, and the width of the seismogenic zone increases northward from the CTJ ( Seismogenic Zone and Seismicity Data , Temperature Field ). This is consistent with the numerical simulation results of Völker et al. (2011), the width of the seismogenic zone Wang et al. (2007) inferred from GPS measurements, and the earthquake distribution to the north of the CTJ (Figure 1).

Wedge subduction significantly affects the seismogenic zone. Our simulation results show that the seismogenic zone of the Antarctic slab is narrower than that of the Nazca slab when the slabs are young. The seismogenic zone of profile 1 is 15 km narrower than that of profile 2 (Table 1). Then the seismogenic zone of profile 1 is narrower than that of the Nazca slab because the seismogenic zone of profile 2 is almost the narrowest to the north of the CTJ. The slab between the CTJ and profile 1 is younger than that of profile 1, so its seismogenic zone is narrower than that of profile 1 and much narrower than that of the Nazca Plate. Our simulations by considering exclusively the slab age and convergence rate show that the width of the seismogenic zone increases southward from the CTJ to profile 1. In addition, a few earthquakes occur to the south of the CTJ according to statistics (Murdie et al., 1993). South of 43°S, the WB plane is not observable for the young age of the Antarctic slab, as proposed by Kirby et al. (1996), Lange et al. (2007), and Agurto-Detzel et al. (2014). This hypothesis is also supported by our finding that the younger slab corresponds to a narrower seismogenic zone. Besides, the oldest slab age of the Antarctic Plate (at the trench) is around 30 Ma, as Figure 1 shows. Therefore, old slabs do not exist south of the CTJ. In fact, the spreading ridge raises the temperature of the South America Plate to a much higher degree and leads to a much narrower seismogenic zone in the vicinity of the CTJ than our simulation results indicate. This may cause difficulties in observing the WB plane, while wedge subduction is not conducive to this phenomenon. Thus, the wedge subduction, slab age, and convergence rate have weaker effects than ridge subduction on the width of the seismogenic zone.

The focal depth is generally shallower south of, than north of, the CTJ (Agurto-Detzel et al., 2014). Because of the wedge subduction, the slab dip of the Antarctic Plate is smaller than that of the Nazca Plate (Slab dip section). Although the Antarctic slab is hotter and has a smaller dip than the Nazca slab, our results suggest that there is little difference between the depths of the seismogenic zones of wedge subduction and slab subduction (Table 1). The downdip limit depth of the seismogenic zone of Model 2 is approximately 1 km deeper than that of profile 2 (Model 4), but 10 km deeper than that of profile 1 (Model 1). Therefore, the effects of slab age and wedge subduction on the focal depth are smaller than that of convergence rate, and the difference between the convergence rate of the Antarctic Plate and that of the Nazca Plate is a main cause of the different focal depths to the south and north of the CTJ. To the north of the CTJ, slip distributions of a 900-km-long rupture zone terminate at a similar depth were inversed by Barrientos and Ward (1990) and Moreno et al. (2009) (Völker et al., 2011). Our results show that the downdip limit of the seismogenic zone increases by only about 2 km, while the slab age increases by 30 Ma, which is consistent with the results of Barrientos and Ward (1990) and Moreno et al. (2009). However, this is inconsistent with the statistics related to the focal depths of the area between 36°S and 40°S (Bohm et al., 2002). In conclusion, the CTJ is not only a wedge subduction boundary but also an important factor controlling the lithospheric thermal structure of the overriding plate. Without considering the spreading ridge beneath the South American plate, the differences in the slab dip, seismic distribution, and focal depth between the northern and southern sides of the CTJ are consistent with the results of our simulation controlled by the plate age, the wedge subduction, and the convergence rate.

Several processes that occur in the subduction zone, such as shear heating, partial melting, hydrothermal cooling, and dehydration, were not covered in this study. Dehydration and partial melting are related to element recycling, intermediate and deep earthquakes, and double seismic zones (Seno et al., 2001; Li and Ni, 2020). Dehydration is essentially a result of the slab responding to the temperature and pressure change during subduction. Dehydration has significantly different effects in cold and hot subduction zones (Li and Ni, 2020). The subduction zone to the north of the CTJ is cold and that to the south is hot. The slab north of the CTJ can, thus, retain more water at greater depths, whereas the slab south of the CTJ dehydrates extensively at shallower depths. Shear heating and hydrothermal cooling are related to the thermal structure of the subduction zone. These processes and their occurrences north of the CTJ are discussed in Völker et al. (2011). The shear heating may give rise to a higher heat flow and a hotter slab south of the CTJ, corresponding to a narrower seismogenic zone; in contrast, the effect of hydrothermal circulation is the opposite of it, leading to a lower heat flow and a colder slab.

The observed data of the area to the south of the CTJ are scant. In this paper, effects of the wedge subduction of Antarctica plate are determined based on 2D numerical simulations. However, south of the CTJ, the thermal structure and seismicity near the CTJ are mainly affected by the spreading ridge subduction, which is a very complex 3D process. In this study, we discuss the effect of wedge subduction on the vicinity of the CTJ only based on the results obtained for the area far away from the CTJ.

The geological conditions in the CTJ area are complex, particularly the lithospheric thermal structure, and the distribution of seismogenic zones are affected by a series of factors such as the ridge subduction, the slab window, the plate age, the slab dip, the hydrothermal circulation, volcanic activities, and the uneven thickness of the lithosphere. With the ridge subduction and other factors ruled out, the effects of the slab age and convergence rate on the lithospheric thermal structure and seismogenic zone were investigated in this paper. The results of our simulations, however, are still complex, especially those obtained for the forearc heat flow. In fact, compared with the plate age and convergence rate, the ridge subduction and the slab window have greater effects of on the thermal structure of the lithosphere of the area south of the CTJ, where the WB plane cannot be observed. Although a series of studies showed that the definition of the seismogenic zone in Klotz et al. (2006) is in good agreement with geological observations on the north of the CTJ, several limitations and uncertainties remain with respect to the mechanism controlling interplate earthquakes. For future research, we will consider the effects of the different factors during the ridge subduction in the CTJ area.