Cooling of the continental plate during flat-slab subduction

The thermal structure of continental lithosphere is significantly modified by a nearby subduction zone. In areas where the oceanic plate descends steeply into the mantle, partial melting in the mantle wedge generates a volcanic arc on the overlying continental plate, which is located 150–300 km from the trench (Turcotte and Schubert, 2002; Syracuse and Abers, 2006). The volcanic arc divides the thermal regime of the continent into a cool forearc region and a hot back-arc region (Dumitru et al., 1991). The forearc temperatures are primarily controlled by conductive cooling by the underlying subducting plate (Wada and Wang, 2009). Conversely, in the arc and back-arc, the presence of a hot mantle wedge and the thermal effects of magmatism lead to high continental temperatures (e.g., Currie and Hyndman, 2006). As a result, back-arcs have a high surface heat flow of 80 ± 20 mW/m²; in contrast, most forearcs have a heat flow of 30–50 mW/m² (Gutscher, 2002; Hyndman et al., 2005; Currie and Hyndman, 2006; Wada and Wang, 2009; Marot et al., 2014).

In some subduction zones, the subducting oceanic plate bends in the shallow mantle to become subhorizontal below the continental plate instead of continuing to descend steeply into the mantle. In these regions, the flat slab is in contact with the overlying continental lithosphere and displaces the hot mantle wedge (Lipman et al., 1972; Dumitru et al., 1991; Humphreys et al., 2003; Axen et al., 2018). Flat-slab regions are generally characterized by the absence of an active volcanic arc (e.g., Kay and Abbruzzi, 1996; Gutscher et al., 2000; Ramos et al., 2002) and relatively low surface heat flow (20–70 mW/m²; Henry and Pollack, 1988; Hamza and Muñoz, 1996; Marot et al., 2014; Sánchez et al., 2018). The low heat flow is usually explained by cooling of the continental plate by the cool flat slab at its base (e.g., Dumitru et al., 1991; van Hunen et al., 2004; Currie and Hyndman, 2006; Manea et al., 2017).

Several episodes of flat-slab subduction have occurred in both North and South America since the Late Cretaceous (e.g., Lipman et al., 1972; Coney and Reynolds, 1977; Gutscher et al., 2000; Ramos and Folguera, 2009). In the Pampean flat-slab region in central Chile (~27–33°S; Fig. 1), seismic data shows that the subducted Nazca plate becomes flat at a depth of 90–110 km at distances of 300 km to 600 km from the trench before it bends downward into the mantle (Anderson et al., 2007; Gans et al., 2011; Porter et al., 2012; Marot et al., 2014). This area is marked by an absence of present-day arc magmatism (Kay and Abbruzzi, 1996). Based on the magmatic record, it is inferred that the onset of flat-slab formation at current latitudes began at ca. 11 Ma (Yáñez et al., 2001) and reached its current geometry by ca. 5 Ma (e.g., Gutscher et al., 2000; Kay and Mpodozis, 2002; Löbens et al., 2011). The development of flat-slab subduction coincides with an eastward expansion of magmatism and the full cessation of volcanism in the last 5–2 Ma (e.g., Urbina et al., 1997; Kay and Mpodozis, 2002). The formation of the flat slab is partially attributed to the subduction of the Juan Fernandez Ridge, which is assumed to be buoyant due to its comparatively thick oceanic crust (von Huene et al., 1997; Ramos et al., 2002; Kopp et al., 2004). Trenchward motion of the continent and increasing suction force due to a thick cratonic root nearby likely also contribute to the formation of flat-slab subduction (e.g., O’Driscoll et al., 2009; Manea et al., 2012).

There are limited thermal data for the Pampean flat-slab region (Fig. 1). Surface heat flow has been determined directly through 32 borehole temperature measurements and indirectly through Curie depth estimates from magnetic data (e.g., Uyeda and Watanabe, 1982; Hamza and Muñoz, 1996; Collo et al., 2018; Sánchez et al., 2018). Above the flat slab, surface heat flow is 40–60 mW/m² (Fig. 1B; e.g., Hamza and Muñoz, 1996; Collo et al., 2018), which is lower than the typical back-arc values of 60–100 mW/m² (e.g., Gutscher, 2002; Hyndman et al., 2005; Currie and Hyndman, 2006). The low heat flow suggests a cool Pampean lithosphere, but there is debate about the extent of cooling. Heat flow measurements have uncertainties of at least ± 15% (Lewis et al., 2003). Furthermore, basin thermochronology measurements by Stevens Goddard and Carrapa (2018) indicate that the region above the Pampean flat-slab has a relatively high geothermal gradient that corresponds to a surface heat flow of 58–80 mW/m², and there is little evidence of cooling since the Miocene. This is higher than the heat flow values of <60 mW/m² from Collo et al. (2011, 2017, 2018). Stevens Goddard and Carrapa (2018) propose that the lack of significant cooling is due to the short time since the formation of the flat slab.

Other flat-slab regions also appear to have relatively cool continental lithosphere. Two well-studied modern flat slabs are in central Mexico and Peru. In both places, the subducted plates are inferred to have developed flat geometries in the last 20 Ma, and the flat slabs are found at depths of ~45 km in Mexico and 60–70 km in Peru (e.g., Ferrari et al., 1999; Pérez‐Campos et al., 2008; Hampel, 2002; Antonijevic et al., 2015; Kumar et al., 2016; Bishop et al., 2017). Neither area has an active volcanic arc above the flat portions (e.g., Ferrari, 2004; Rosenbaum et al., 2005). In central Mexico, the surface heat flow is 13–40 mW/m² between the trench and the Trans-Mexico Volcanic Belt (~350 km from the trench), where the slab resumes a steep angle (Ziagos et al., 1985). In Peru, the heat flow above the Peruvian flat-slab is 20–60 mW/m², which is lower than that of the back-arc regions above the steep-angle slabs to the north and south (80–120 mW/m²) (Henry and Pollack, 1988; Muñoz, 2005).

In the western United States, Farallon flat-slab subduction is inferred to have induced widespread cooling of the continental lithosphere from the Late Cretaceous to Eocene (Dumitru et al., 1991; Smith, 2020). The subducted Farallon plate started to develop a flat geometry at ca. 90 Ma, which eventually created a flat-slab segment at 120–150 km depth that extended more than 1500 km from the plate margin and persisted until ca. 50 Ma (e.g., Coney and Reynolds, 1977; Copeland et al., 2017). Thermochronological data in the region of the volcanic arc show that cooling started at 5–10 million years (m.y.) after the local cessation of magmatism (Dumitru et al., 1991). Fission-track data indicate that upper crust temperature decreased by ~200 °C over 10 m.y. after the onset of cooling (Dumitru, 1990). Farther into the continental interior, thermochronological data also show cooling of the shallow continent, but the magnitude is reduced (Dumitru et al., 1991). In addition, cooling in the continental interior started later and occurred more gradually than in the former arc region. This is interpreted to reflect both the extra time it took for the Farallon flat-slab to reach the continental interior and the fact that the flat slab was at a greater depth here; therefore, it took longer for the thermal effects of the cool slab to reach the shallow crust (Bird, 1988; Dumitru et al., 1991).

The observations above support the idea that flat-slab subduction regions are associated with a relatively cool lithosphere compared to that of normal-slab regions. This is consistent with the idea that the flat slab displaces the hot mantle wedge, which leads to the cessation of arc volcanism and conductive cooling of the continental lithosphere due to the presence of the cool oceanic plate (Dumitru et al., 1991; Wada and Wang, 2009; Currie and Hyndman, 2006). However, cooling through thermal conduction is a relatively slow process (Dumitru, 1990). The present-day flat slabs only achieved the horizontal geometries within the last 20 Ma (Ferrari et al., 1999; Yáñez et al., 2001; Hampel, 2002; Ramos et al., 2002; Rosenbaum et al., 2005). Therefore, it is unclear whether there has been sufficient time for the continent to cool. In this study, we use numerical models to assess temporal variations in the continental thermal structure during the development of a flat slab.

The two-dimensional numerical models are based on our previous work on the development of a flat slab and the factors controlling its depth (Liu and Currie, 2019). The thermal-mechanical evolution of the model is calculated using the finite-element code SOPALE (Fullsack, 1995). It solves the equations of conservation of mass, force balance, and energy balance assuming plane strain and incompressibility (Fullsack, 1995).

The model domain is 3000 km wide and 660 km deep (Fig. 2) and is oriented parallel to the convergence direction. The Eulerian calculation mesh has elements with a horizontal width of 10 km and a height of 2 km for the top 90 km, 5 km for the middle 100 km, and 10 km for the bottom 470 km. Within the domain, we model subduction of an oceanic plate below a continental plate, where the trench position is located 400 km from the left boundary. In the following sections, we report all distances with respect to the trench.

The geometry of the reference model is based on the present-day subduction zones in South America. The oceanic plate has a 70-km-thick lithosphere and a thermal structure consistent with a 30 Ma oceanic plate, which is based on the average oceanic plate age in the flat subduction regions in Mexico and South America (e.g., Müller et al., 2008; Quinteros and Sobolev, 2013). We include a 600-km-wide oceanic plateau. Our modeled plateau has an 18-km-thick crust that corresponds to the average thicknesses of the ridges in the Peruvian (the Nazca Ridge) and Pampean (the Juan Fernandez Ridge) flat-slab regions (e.g., Woods and Okal, 1994; Yáñez et al., 2002; Gans et al., 2011). The adjacent normal oceanic crust is 8 km thick. Below the oceanic crust, there is a harzburgite layer that is twice the thickness of the crust to represent melt-depleted mantle (van Hunen et al., 2004). For simplicity, we use a laterally uniform continental plate. In our reference model, we use a 48 km crustal thickness and a total lithospheric thickness of 120 km. In later models, we vary both the crustal and lithospheric thickness.

The material properties are given in Table 1 and follow values used in earlier studies (e.g., Beaumont et al., 2006; Liu and Currie, 2019). All materials undergo viscous-plastic deformation. The rheologies are based on the laboratory experiments, and we apply a linear scaling factor (f) to vary the effective viscosity to allow for strength variations relative to the laboratory materials (Beaumont et al., 2006). We also include frictional-plastic strain softening and viscous strain weakening to represent processes that cause a decrease in material strength during deformation such as changes in grain size (e.g., Beaumont et al., 2006). The exception is the continental crust, which has a dry, strong continental crust without strain softening or weakening (Table 1) so that deformation is minimized, and the resulting thermal structure only depends on the effects of the flat slab. The continental upper crust has the viscous rheology of wet quartzite (Gleason and Tullis, 1995) with f = 50, and the lower crust has the viscous rheology of dry Maryland diabase (Mackwell et al., 1998) with f = 1. The oceanic crust uses the dry Maryland diabase rheology with f = 0.1. The oceanic and continental mantle lithosphere have a wet olivine rheology (Karato and Wu, 1993) with f = 10 to approximate dry conditions. The sublithospheric mantle has a wet olivine rheology without any scaling (f = 1).

The thermal conductivity of all materials is 2.25 W/(m·K), and the upper and lower continental crusts have radiogenic heat production of 1.0 μW/m³ and 0.4 μW/m³, respectively; there is no heat production for other materials (Jaupart and Mareschal, 2014). The density of all materials is temperature-dependent. The oceanic crust undergoes a phase change when the temperature and pressure conditions reach the eclogite stability field (Hacker et al., 2003). The phase change is included as an increase of the reference density from 2950 kg/m³ to 3500 kg/m³ with no other changes (e.g., rheology). Normal oceanic plate fully transforms to high-density eclogite, whereas the oceanic plateau crust undergoes 30% eclogitization (reference density of 3125 kg/m³). This allows the plateau to remain buoyant enough to develop a flat slab; the assumption is that this region is relatively dry and therefore does not fully transform to eclogite (Liu and Currie, 2019). The amount of eclogitization does not change the flat-slab subduction dynamics as long as the plateau remains less dense than the shallow sublithospheric mantle (Liu and Currie, 2019).

Figure 2A shows the boundary conditions used in the models. The temperatures of the top and bottom boundaries are 0 °C and 1564 °C, respectively. The bottom temperature is calculated based on a 1300 °C mantle adiabat (vertical gradient of 0.4 °C/km). At the left boundary, the oceanic lithosphere has a geotherm that is consistent with a 30 Ma plate (Stein and Stein, 1992; Fig. 2B). The side boundaries are insulated boundaries with no horizontal heat flux. The initial structure of the model is calculated based on these boundary conditions and the thermal properties of each material. This results in a continental geotherm (Fig. 2B) that is consistent with its 120 km thickness and a surface heat flow of 57.4 mW/m². In later sections, we vary the thickness and the associated thermal structure of the continental lithosphere. The top surface is a stress-free boundary that allows for topography development. The bottom boundary is a closed, free-slip boundary. The side boundaries are no-slip boundaries, and they have prescribed horizontal velocities to generate plate convergence.

The average convergence rate between the Nazca and South American plates over the last 15 Ma was ~10 cm/yr (Norabuena et al., 1999). In the model, the oceanic plate enters the model domain through the left boundary at 7 cm/yr (V_o in Fig. 2A). The continent has a trenchward motion of 3 cm/yr (V_c), which is consistent with motion of South America relative to the deep mantle (Quinteros and Sobolev, 2013). To maintain mass balance of the model, a small outflux velocity (V_b) is assigned to the side boundaries of the sublithospheric mantle. Models are solved in the continental reference frame, and therefore V_c is added to all side boundaries (Liu and Currie, 2016).

Models are run in three phases. In the first two phases, the models undergo isostatic adjustment and subduction initiation for 500 km plate convergence. At the end of phase two, the oceanic plateau is located at the trench with a normal-angle subduction zone. Phase three starts at this point (the 0 m.y. panel of Fig. 3), and all times are reported relative to this time. Models are run for 30 m.y. in phase three to study flat-slab subduction over this timeframe, but we note that the Pampean flat-slab at its current latitude only initiated in the last 11 Ma (Yáñez et al., 2001), and the development of the Peruvian flat-slab began at ca. 11 Ma (Hampel, 2002). In phase three, we examine different parameters to test their effects on the continental thermal structure. All models are listed in Table 2.

We first present the reference model, which has an initial continental lithospheric thickness of 120 km and continental crustal thickness of 48 km. Figure 3 shows the model evolution from 0 m.y. to 20 m.y. The models are run for an additional 10 m.y. in which the slab maintains a flat geometry with the edge of the slab reaching ~2800 km. At 0 m.y., the oceanic plateau enters the trench. At this time, the subduction zone has a normal-angle slab and a hot mantle wedge. We expect that this would lead to a volcanic arc ~200 km from the trench (Fig. 4A) based on where the mantle wedge temperature is >1200 °C (e.g., Schmidt and Poli, 1998). The Moho temperature and surface heat flow are high near the putative arc due to mantle wedge corner flow (Fig. 4A). In the forearc region (<100 km from trench), the cool oceanic plate is in contact with the continent at shallow depth (<50 km), and the surface heat flow is <40 mW/m², which is consistent with observations in forearc regions (e.g., Wada and Wang, 2009).

As the oceanic plateau enters the trench, it remains buoyant relative to the mantle. This causes a decrease in subduction angle. Owing to the density contrast between the oceanic plateau and the oceanic plate that has already subducted, the slab undergoes a break-off at ~5 m.y. With continued plate convergence, the buoyant oceanic plateau causes the slab to form a flat-slab geometry beneath the continental lithosphere at 90–100 km depth.

As the slab flattens, we expect a cessation of arc volcanism as the hot mantle wedge is displaced. The presence of the flat slab results in rapid cooling in the regions near the predicted arc (~200 km from the trench). The Moho temperature decreases by >200 °C and the surface heat flow decreases by >15 mW/m² by 20 m.y. (Fig. 4A). The magnitude of cooling becomes smaller farther into the continental interior.

Figure 4B shows the vertical temperature profile at 400 km from the trench at 0 m.y., 10 m.y., 20 m.y., and 30 m.y. The flat slab reaches this location at ~5 m.y., and the slab depth is ~90 km. After the slab is in contact with the continental lithosphere, the continent cools through thermal conduction. The lowermost continent cools first and has a greater amount of cooling during the model run than the shallower regions. Figure 4C shows the geothermal gradient near the ground surface (upper 2 km). For the first 20 m.y., the surface geothermal gradient has a relatively constant value of 24.1 °C/km and then decreases to ~23.5 °C/km over the next 10 m.y. with small fluctuations that appear to be due to minor crustal deformation and topographic changes. The near-surface thermal structure is not immediately affected by the presence of a flat slab, and there is an ~15 m.y. delay after the flat-slab emplacement at 5 m.y. This delay reflects the time needed for conductive cooling to propagate from the slab to the shallow crust.

To compare the lateral variations in continental thermal structure over time, we plot the Moho temperature (Fig. 4D) at three locations (400 km, 600 km, and 800 km from the trench). The flat slab reaches the three locations at 5 m.y., 8 m.y., and 11 m.y., and the Moho temperature starts to decrease from its initial value of 745 °C at 13 m.y., 16 m.y., and 19 m.y. with later times for positions farther from the trench. This shows that the onset of cooling at the Moho has a delay of ~8 m.y. following flat-slab emplacement. By the end of the 30 m.y. model run, the Moho temperatures have decreased by ~77 °C, 28 °C, and 15 °C at 400 km, 600 km, and 800 km, respectively. In all locations, the flat-slab depth is approximately the same, and therefore the greater amount of cooling closer to the trench is primarily due to the earlier flat-slab emplacement and longer cooling time. A secondary effect is that the slab temperature is lower below this region. The slab surface temperature increases by ~100 °C every 100 km inboard (Fig. 4A).

Figure 4E shows the temporal changes in surface heat flow at 400 km, 600 km, and 800 km from the trench. This is determined based on the near-surface thermal gradient and the thermal conductivity (2.25 W/(m·K)). The initial surface heat flow is 57.4 mW/m², and the decreases of surface heat flow are 1.5 mW/m², 0.4 mW/m², and 0.2 mW/m² over 30 m.y. at 400 km, 600 km, and 800 km from the trench, respectively (Fig. 4E). This change in heat flow is well below the uncertainties in heat flow measurements (Lewis et al., 2003), and it is unlikely that this could be resolved in observations.

In the following sections, we vary model parameters to examine the effect on the continental thermal structure. To compare the results, we report the Moho temperature and surface heat flow at 400 km from the trench (blue dotted line in Fig. 4A), as this position is in the continental interior where observations show low surface heat flow (Fig. 1).

To demonstrate the effect of subduction geometry on continental thermal structure, we present a model that maintains a normal (steep) subduction angle throughout the 30 m.y. model run. This model uses the same parameters as the reference model except that the oceanic plateau crust undergoes 100% eclogitization, which makes it negatively buoyant, and no flat slab develops.

Figure 5A shows the model evolution. The normal-angle slab results in the persistent presence of mantle wedge with a temperature of >1200 °C above a slab depth of 100 km. Thus, we expect that there would be a continuous magmatic arc ~200 km from the trench. In the forearc (<150 km), the Moho temperature and surface heat flow are below the initial values due to cooling from the subducting plate. At distances of 150–350 km, the continent is heated by mantle wedge flow, and the Moho temperature and surface heat are higher than the initial values; the peaks are at 200 km (950 °C Moho and ~65 mW/m² heat flow).

Figures 5B and 5C compare the models with normal-angle and reference flat-slab subduction at 400 km. The normal-angle subduction model has negligible changes in the continental Moho temperature and surface heat flow during the 30 m.y. model run. This contrasts with the flat-slab model in which the Moho temperature decreases by 77 °C and surface heat flow decreases by 1.5 mW/m². This shows that the continental cooling in the reference model results from the presence of the flat slab.

In this section, we examine variations in initial continental thickness. As we showed in an earlier study (Liu and Currie, 2019), the flat-slab depth is primarily controlled by the initial continental thermal structures owing to the temperature-dependent rheology of mantle lithosphere. As a slab flattens, it is able to displace the lowermost continental mantle lithosphere with temperatures >900 °C. The flat-slab depth increases for thicker continents, as these regions are cool and thus there is a thicker portion of the lithosphere that is too strong to be displaced.

We present models with initial continental thicknesses of 80 km, 120 km (reference model), and 160 km. These variations result in differences in the initial thermal structure, with Moho temperatures of 930 °C, 745 °C, and 680 °C and surface heat flows of 65 mW/m², 57.4 mW/m², and 54.1 mW/m² for models with continental thickness of 80 km, 120 km, and 160 km, respectively. The general evolution of all models is the same; slab flattening is induced by subduction of the oceanic plateau. The main difference is the flat-slab depth. At 400 km from the trench, the flat slabs are at depths of 60 km, 90 km, and 115 km with a deeper slab for a thicker continental plate (Fig. 6A). Figure 6B shows the geotherms of the three models at 400 km over the 30 m.y. model run.

To compare the three models, we plot the Moho temperature change and surface heat flow change at 400 km (Figs. 6C–6D) relative to the initial values for each model. The flat slab reaches this position at 5 m.y. in all models. For the thinnest continent (80 km thick), the Moho temperature starts to decrease at ~8 m.y. and drops by ~450 °C by 30 m.y. In contrast, the continents with thicknesses of 120 km and 160 km experience a Moho temperature decrease of 77 °C and 31 °C starting at 13 m.y. and 16 m.y., respectively. Therefore, the times between flat slab emplacement (5 m.y.) and onset of cooling are 3 m.y., 8 m.y., and 11 m.y. for the three models with the shortest time for the thinnest continent (shallowest flat slab). For the 80 km continent, the surface heat decreases by ~18 mW/m² by 30 m.y. starting at ~12 m.y. (7 m.y. after flat-slab emplacement), whereas the surface heat flow drops by 1.5 mW/m² and 0.5 mW/m² for the 120 km and 160 km continents, respectively, starting at >20 m.y.

These models demonstrate that the flat-slab depth strongly affects the continental cooling. A shallower flat slab leads to both an earlier onset and a greater amount of cooling of the overlying continent because conductive cooling due to the slab transfers through a thinner layer of lithosphere (Dumitru et al., 1991).

This section presents models that examine variations in thermal conductivity, continental crustal thickness, convergence rate, and oceanic plate age. The models are based on the reference model, and only the tested parameter is different from the reference value. In all models, the development of a flat slab is similar to the evolution of the reference model (Fig. 3).

Thermal conductivity controls the rate of conductive heat transfer (Turcotte and Schubert, 2002). The thermal conductivity of lithospheric rocks ranges from 1.5 W/(m·K) to 7 W/(m·K) depending on the rock composition and temperature, and values are higher for cool and mafic rocks (Clauser and Huenges, 1995). We tested thermal conductivities of 2.25 W/(m·K) (reference value), 3.25 W/(m·K), 4.25 W/(m·K), and 5.25 W/(m·K), where this value is applied to the entire continental lithosphere (crust and mantle). All models have the same lithospheric thickness (120 km), but the variation in conductivity results in different initial thermal structures (Fig. 7). Models with higher thermal conductivities have lower lithospheric temperatures so that all models have the same temperature at 120 km depth (i.e., 1348 °C, the mantle adiabat at this depth). As a result, these models have different flat-slab depths. The slab depths are 90 km, 97 km, 105 km, and 110 km, with deeper slabs having a higher thermal conductivity due to the lower temperature and thus stronger lithosphere.

Figure 7A shows the temporal variation in absolute Moho temperature at 400 km from the trench as well as the difference relative to the initial value. The initial slab emplacement at this location is the same (5 m.y.). The model with highest thermal conductivity has a Moho temperature decrease of 103 °C by 30 m.y. compared to a decrease of 77 °C for the reference model with the lowest conductivity. The onset of Moho cooling occurs up to 5 m.y. earlier as the thermal conductivity increases. Figure 7B shows that the surface heat flow changes of the four models decreases by as much as 3.5 mW/m², and the largest decrease is in models with higher thermal conductivity. These results show that continental cooling is enhanced by a higher thermal conductivity, as the cool slab temperatures are more efficiently transferred through the continent even though the slab is at a greater depth.

We test models with continental crustal thicknesses of 36 km, 42 km, and 48 km (reference model) in which the upper crust thickness is twice that of the lower crust. The radiogenic heat productions are 1.0 μW/m³ and 0.4 μW/m³ for the upper crust and lower crust, respectively. All models have the same continental lithospheric thickness of 120 km. Figure 8A shows the Moho temperature and temperature changes over the 30 m.y. model run at 400 km from the trench. The model with the thinnest crust (36 km) has the coolest initial Moho temperature (538 °C), and the flat slab causes the Moho temperature to decrease starting at 17 m.y. (12 m.y. after slab emplacement), with an overall decrease of ~13 °C by 30 m.y. In contrast, the reference models with a 48 km crust experience Moho cooling starting at 13 m.y. (8 m.y. after slab emplacement), and it reaches 77 °C by 30 m.y. The earlier onset and greater magnitude of cooling for a thicker crust occur because the Moho is closer to the cool flat slab. The surface heat flow decrease for all crustal thicknesses is similar to that seen in the reference model (~1.5 mW/m²), and therefore we do not show it here.

In the Pampean flat-slab region, the convergence rate between the Nazca plate and South America plate has slowed from 10 cm/yr to 7 cm/yr in the Neogene (Quinteros and Sobolev, 2013). Here, we present a model with a convergence rate of 7 cm/yr for comparison to the reference model (10 cm/yr). With the slower convergence rate, the flat slab reaches 400 km from the trench at ~7 m.y., compared to 5 m.y. for the reference model, but the overall slab dynamics are similar. Figure 8B shows the temporal changes in Moho temperature of the two models at 400 km from the trench. Both models experience Moho cooling starting ~8 m.y. after flat-slab emplacement. By 30 m.y., the model with slower convergence has a Moho temperature decrease of 50 °C and a surface heat flow decrease of 0.7 mW/m². These changes are lower than those in the reference model (77 °C and 1.5 mW/m²), which shows that the cooling rate is enhanced when the convergence rate is higher. As both models have a similar flat-slab depth (90 km), the cooling appears to be caused by the fact that the faster slab is cooler.

The age of the oceanic plate controls the thermal state of the oceanic lithosphere; an older plate has a cool and thick lithosphere (e.g., Stein and Stein, 1992). The previous models use a young oceanic plate (30 Ma) to represent the young, flat-subducting plates in Mexico and South America (e.g., Müller et al., 2008). The Farallon flat-slab in the Cretaceous western United States had a 100 Ma oceanic plate (Usui et al., 2003). To show the effect of subducting plate age on the continental thermal structure, we test a model with an oceanic plate age of 100 Ma. This oceanic plate has a thickness of 90 km based on the GDH1 plate cooling model (Stein and Stein, 1992). Figure 8C shows the continental Moho temperature over the 30 m.y. model run for these models. In both models, the Moho temperature starts to decrease ~8 m.y. after flat slab emplacement. With an old oceanic plate, the Moho temperature decreases by 105 °C at 30 m.y., which is greater than that for the young oceanic plate (77 °C). Both models have the same convergence rate and flat-slab depth, and therefore the difference is due to the oceanic plate temperatures. With an older plate, the flat slab is cooler, and this results in a higher cooling rate as seen by the greater slope in Figure 8C.

The thermal-mechanical models allow us to examine the temporal changes in continental thermal structure during flat-slab subduction. Table 2 summarizes the results. The model with normal-angle subduction does not show a significant change in back-arc thermal structure over the 30 m.y. model run. All models with a flat slab show some continental cooling over time. In these models, the emplacement of a cool flat slab beneath the continent displaces the hot mantle wedge. The cold slab is in contact with the continent and cools the overlying lithosphere through thermal conduction. The cooling has a greater impact on areas closer to the trench (Figs. 4–5), as the flat slab reaches these areas earlier, which allows for a longer cooling time. A secondary effect is that the slab has a lower temperature in this region. The bottom of the continental lithosphere is most affected by the flat slab because it is closest to the cool slab surface and cools first. Similarly, the continental Moho temperature starts to decrease before near-surface temperatures, and the overall magnitude of cooling is greater at the Moho; this is especially true for models with a thick crust (Fig. 8A).

The main factor that controls the cooling of the continent is the flat-slab depth (Fig. 6). With a shallower flat slab, the continental crust experiences an earlier onset and greater magnitude of cooling compared to deeper slabs. Thus, less time is needed for reduced temperatures to be observed at the surface (Dumitru et al., 1991). Other parameters have a second-order influence on the continental thermal structure. Higher thermal conductivity, faster convergence rates, and older oceanic lithosphere result in slightly greater cooling by either increasing the heat transfer rate or decreasing the slab temperature.

All of the models presented in this work are in 2-D, and therefore they do not allow an assessment of along strike variation in slab dip and thermal structure. In a 3-D flat-slab subduction setup, if a tear occurs at the edge of the flat segments, mantle flow around the slab edge may bring hot asthenosphere into contact with the continent, which leads to heating (e.g., Hu and Liu, 2016; Martinod et al., 2005). The heat may also conductively transfer laterally into the flat-slab region and may slow the cooling process of this region near the flat slab.

We also note that in our models, a flat slab develops after the break-off of the dense slab in front of the oceanic plateau (Fig. 3). Other models of flat-slab subduction show that a flat slab can develop without a break-off, such that there is a continuous slab with a steep section ahead of the flat segment (e.g., van Hunen et al., 2004; Liu and Currie, 2016). Seismic tomography studies of modern flat-slab regions (Mexico, Peru, and Pampea) show that the geometry may be continuous (Kim et al., 2012; Scire et al., 2016; Portner et al., 2017), although a slab gap is observed below the southern part of the Pampean flat-slab (Portner et al., 2017). With a continuous slab, slab flattening is accompanied by an inboard migration of the hot mantle wedge corner. This may result in a migration of the volcanic arc, where magmatism could cause enhanced heating of the shallow continent for a few million years following flat-slab emplacement. This would make the continent hotter than we have assumed here, which would delay slab cooling. Alternately, as the flat slab advances, the mantle wedge may be filled by continental mantle lithosphere that is displaced by the flat slab, which would fully terminate magmatism (Axen et al., 2018). Future models that include magmatism are needed to assess these effects.

To assess the timescale for continental cooling and to explore a greater range of parameters, we present 1-D temperature calculations. These solve the time-dependent heat conduction equation using the conservative finite difference method (Gerya, 2019). The initial 1-D temperature structure of the continent is given by a geotherm based on the assumed continental thickness using the same thermal parameters as those of the numerical models (Table 1). A cool flat slab is emplaced at a prescribed depth and with a constant temperature. The resulting thermal structure for the overlying material is then calculated for a longer time of 60 m.y. compared to the 30 m.y. in the 2-D numerical models. This encompasses the long-lived flat slabs below the Late Cretaceous–Eocene western United States (40 m.y.; e.g., Copeland et al., 2017) and Mesozoic South China (60 m.y.; e.g., Li and Li, 2007).

Figure 9A shows the evolving 1-D thermal structure for a 120-km-thick continent with a 600 °C slab at 90 km depth based on the reference numerical model at 400 km from the trench (Fig. 3). Owing to the low-temperature slab, the continent starts to cool, and the greatest cooling occurs in the deep lithosphere. Figures 9B–9C show the temporal changes in surface heat flow and Moho temperature. The Moho temperature starts to decrease 8 m.y. after emplacement of the flat slab and cools by 200 °C at 60 m.y. In contrast, the surface heat flow starts to decrease at 17 m.y. and drops by ~10 mW/m² at 60 m.y. The Moho temperature and surface heat flow variations from the 1-D calculations agree well with those of the reference numerical model. If the 1-D calculations were run for longer, the geotherm would reach steady-state conditions (dashed lines in Fig. 9) with surface heat flow and Moho temperature of 45 mW/m² and 479 °C, respectively.

Based on Figure 9, continental cooling can be divided into three stages that are defined by two critical times: τ_l (the lag time) and τ_c (the cooling time or thermal time constant). In the early stages of the flat slab (at times less than τ_l), the thermal structure of the continent retains its initial values. At time τ_l, conductive cooling starts to decrease the temperature. This marks the start of the second stage, where temperatures are determined by both the initial continental structure and the flat-slab depth and temperature. In the third stage (at times greater than τ_c), the continental cooling approaches the steady-state values, and the temperatures are dominated by the effect of the flat slab.

For the reference model (Fig. 9), the lag time (τ_l) occurs at 8.3 m.y. for the Moho temperature and 16.7 m.y. for the surface heat flow. The thermal time constant (τ_c) is 46.8 m.y. for the Moho temperature and 59.3 m.y. for the surface heat flow. The longer timescales for the surface heat flow are due to the fact that the ground surface is further from the slab.

For comparison, we show the thermal evolution for a 90-km-deep slab with slab temperatures of 400 °C and 800 °C (blue and green dotted lines in Figs. 9B–9C). Note that τ_l and τ_c are the same as in the model with a 600 °C slab because the slab is at the same depth in all models. However, at times greater than τ_l, the temperatures depend on the slab temperature. When the slab is 400 °C, the Moho cools by an extra ~80 °C, and the surface heat flow decreases by an extra 3 mW/m² over 60 m.y. compared to the reference 600 °C slab. This is consistent with the greater cooling observed in models that had an older oceanic plate or faster convergence rate (Table 2).

We also examine the effects of different slab depths, from 50 km to 120 km, for a 120 km continental thickness and 600 °C slab temperature. The presence of the cool slab causes the continent to cool from its initial conditions to steady-state values that depend on the slab depth (Fig. 10). The timescale for cooling also depends on the slab depth; the surface heat flow lag time increases from τ_l = 4 m.y. (50 km slab) to τ_l = 32 m.y. (120 km slab). Similarly, the thermal time constant (τ_c) increases with slab depth, and where the slab is shallower than 60 km, the Moho temperature and surface heat flow reach steady-state within ~30 m.y. Note that both τ_l and τ_c vary nonlinearly with slab depth, as the timescales for thermal conduction are proportional to the square of the thickness of the lithosphere above the slab surface (Dumitru, 1990; Turcotte and Schubert, 2002).

In the final calculations, we demonstrate the effects of initial continental structure on the cooling when the slab is at a very shallow depth (short cooling timescale). We test two continental thicknesses of 60 km and 120 km with a 600 °C slab at a depth of 50 km (2 km below the Moho). The thinner continent has a hotter initial structure than the reference 120 km continent with initial Moho temperature of 1129 °C and surface heat flow of 75.3 mW/m². The timescales for cooling are the same as for the two models as they have the same slab depth (Fig. 11). At the Moho, τ_l = 0.026 m.y. and τ_c = 0.031 m.y., and in both models, the Moho temperature rapidly decreases to the steady-state value of 590 °C (Fig. 11A). The surface heat flow starts to decrease at τ_l = 4.1 m.y. (Fig. 11B). At times greater than τ_l, the heat flow is greater for the 60 km continent owing to its initially hotter structure, but both models are close to the steady-state value of 50 mW/m² within 30 m.y.

Our numerical models provide insight into the interpretation of surface heat flow observations for modern flat-slab regions. The critical variables are (1) the flat-slab depth and (2) the length of time the flat slab has been emplaced below the heat flow measurement point. For shallow slabs like the Mexican and Peruvian flat-slabs that have existed for <20 Ma, our results indicate that conductive cooling due to the underlying slab is expected to reduce surface heat flow measurements. In central Mexico, observations show very low surface heat flow (13–40 mW/m²; Ziagos et al., 1985) above a flat slab located at ~45 km depth (Manea et al., 2017). Our results indicate that a flat slab at ~50 km depth would reduce surface heat flow by over 20 mW/m² within 20 m.y. through thermal conductive heat transfer (Fig. 11B). However, 20 mW/ m² may be insufficient to explain the very low heat flow measurements even when error is considered. One possibility is that the flat slab is cooler than the 600 °C assumed in our calculations despite the young slab age and relatively slow convergence rate. Perry et al. (2016) suggest that the hydrothermal circulation within the oceanic crust of the Cocos plate could cool the subduction zone by up to 180 °C, which could account for this difference.

Similar to the Mexican flat-slab, the flat slab in southern Peru also has a shallow slab depth of 60–70 km and is immediately beneath the continental Moho (Kumar et al., 2016; Bishop et al., 2017). Again, the shallow slab decreases the cooling lag time and increases the cooling magnitude. However, the Peruvian flat-slab is migrating southward along the trench in addition to its trench-perpendicular growth direction (Hampel, 2002; Rosenbaum et al., 2005). This suggests that a thermal gradient north-to-south should be observable, and the continent above the older flat slab to the north should be colder than the continent above the younger flat slab to the south. Unfortunately, there are insufficient heat flow measurements to fully explore this hypothesis. The existing heat flow measurements (20–60 mW/m²; Currie and Hyndman, 2006) are generally consistent with the cooling that we predict, but more work is needed to fully incorporate the along-strike evolution of this complicated flat slab along with any other factors that contribute to continental heat flow.

In contrast with shallow flat slab regions, low surface heat flow measurements above deeper flat slabs are not well explained by our models. In the Pampean flat-slab region, the flat slab lies at 90–110 km depth and was emplaced at <11 Ma (e.g., Yáñez et al., 2001; Kay and Mpodozis, 2002). Both our reference numerical model and our 1-D calculations show that within the ~10 m.y. existence time of the Pampean flat-slab, the surface heat flow should not decrease through conductive heat transfer (Figs. 4, 9, and 10). However, observed surface heat flow is 40–60 mW/m², which is ~20 mW/m² less than typical back-arc values (e.g., Hamza and Muñoz, 1996; Collo et al., 2018; Fig. 1) despite the limited heat flow measurements. We suggest that more heat flow measurements are needed in this region for future studies.

Some recent work suggests that heat flow measurements could be artificially suppressed using water during borehole drilling and that a new assessment of thermochronological data suggests much more average heat flow above the Pampean flat-slab (Stevens Goddard and Carrapa, 2018). If that is the case, this would be consistent with our results. However, the low heat flow measurements cited above are consistent with other previous studies using Curie point temperatures (Ruiz and Introcaso, 2004; Sánchez et al., 2018), crustal seismicity depths (Alvarado et al., 2009), and thermochronological data (e.g., Collo et al., 2011, 2017). Collectively, these studies suggest relatively cool temperatures. One possibility is that these cool temperatures reflect processes beyond conductive cooling. For example, near surface hydrological processes, such a changing groundwater flow pattern and local permeability could affect the surface heat flow significantly; low surface heat flow is observed in high recharge upland regions (Burns et al., 2015).

The surface heat flow is also affected by the initial continental structure prior to flat-slab subduction. Both initial thickness of the lithosphere and crustal radiogenic heat production affect the thermal structure. The low heat flow observations near the Pampean flat-slab (Fig. 1) are relatively uniform from the longitude of the paleo-arc to the east across the flat-slab region and into the Paleoproterozoic Rio de la Plata craton (Oyhantçabal et al., 2011). The flat slab traverses a series of north-south–trending, accreted terranes and associated arcs that include the Cuyania, Famatina, and Pampian terranes. The Cuyania terrane is likely of Laurentian origin and may have drifted from the Ouachita Embayment prior to its accretion to South America in the mid-Ordovician (Thomas and Astini, 1996; Ramos, 2004; Ramos et al., 2010). The Famatinian arc developed on continental crust during the early to mid-Ordovician with the accretion of the Cuyania terrane onto the Pampian basement (Pankhurst et al., 1998; Miller and Söllner, 2005; Ramos et al., 2010). The Pampian terrane is believed to have been part of Rodinia and likely accreted onto the western margin of the Rio de la Plata craton in the late Proterozoic (Ramos, 1988; Ramos et al., 2010). Given the diverse origins of these terranes, it is unlikely that the largely uniform heat flow observations are due solely to uniformly low radiogenic heat production. While a thick cratonic root may explain the low heat flow measurements, the ~100 km depth of the modern flat slab makes it unlikely that the overlying terranes had a much thicker lithosphere prior to flat-slab emplacement (Liu and Currie, 2019). We therefore do not ascribe the low heat flow measurements to pre-flat-slab tectonic structures.

Our models only consider conductive heat transfer from the cool flat slab. One possible additional cooling mechanism is flat-slab dehydration. Downgoing oceanic plates are believed to be extensively hydrated to depths of up to 30 km along outer-rise faults (Ranero et al., 2003; Fromm et al., 2006; Bloch et al., 2018; Cai et al., 2018; Grevemeyer et al., 2018; Wagner et al., 2020). The release of this water is believed to play a role in the generation of intermediate depth seismicity (e.g., Meade and Jeanloz, 1991; Ferrand et al., 2017) and flux melting in normal subduction zones. The depth of this release in normal subduction zones depends dominantly on the temperature of the downgoing plate, which in turn depends on its age, convergence rate, and dip angle (e.g., Syracuse et al., 2010). In the case of flat slabs, a steady-state thermal profile is not achieved within the lifetime of known flat slabs, and therefore we have the added variable of the duration of flat-slab subduction beneath a given portion of the overriding plate. In the case of the Pampean flat-slab, there is copious intermediate depth seismicity along the shallowest portions of the flat slab (e.g., Anderson et al., 2007), which suggests that the horizontally lying slab is actively dehydrating. However, numerous tomographic studies (Wagner et al., 2005, 2006, 2008; Porter et al., 2012; Marot et al., 2014) have shown evidence of high shear wave velocities (>4.7 km/s) above the seismically active portion of the flat slab. These high shear wave velocities are difficult to reconcile with the presence of any significant amount of serpentinite (Hacker and Abers, 2004; Abers and Hacker, 2016). This suggests that the mantle above the flat slab may be too hot for hydrous phases to be stable. In this case, water released from the flat slab would pass directly through the mantle above it and into the overriding crust. Previous studies have shown that aqueous fluids are able to migrate through the mantle lithosphere above ~110 km depth (e.g., Peacock and Hyndman, 1999; Humphreys et al., 2003; Cerpa et al., 2017). This flux of hydrous fluids may produce an additional mechanism for heat advection that could hasten the cooling of regions of deep flat slab subduction. The details of any such advective cooling would depend heavily on the rates of fluid migration and the fluid pathway, which are topics for future research.

In the western United States, it is believed the Farallon plate developed a flat geometry in the Late Cretaceous (Coney and Reynolds, 1977; Copeland et al., 2017). Thermochronology data show there was widespread cooling of the Cordilleran lithosphere during Farallon flat subduction (Dumitru et al., 1991). Within 200 km of the trench, thermochronology data indicate that the geothermal gradient decreased to one-fifth of the initial value of 50 °C/km and cooling was delayed by 5–10 m.y. after the cessation of arc magmatism (Dumitru, 1990). Dumitru's (1990) model shows that the delay time is dependent on the thickness of the continental mantle lithosphere above the slab surface and predicts that the slab here was at a depth of 35–50 km. This agrees with our results of significant surface cooling 5–10 m.y. after flat-slab emplacement for a slab at ~50 km depth (Fig. 10). Thermochronology data also show that the amount of cooling in the western United States decreased from the trench to the continental interior with a later onset of cooling further from the plate margin (Dumitru et al., 1991). Our model results are consistent with this observation that the region near the trench has the earliest and greatest amount of cooling and that the cooling decreases with distance landward from the trench (Fig. 3). The widespread cooling in the western United States was likely caused by the long duration of the Farallon flat-slab (~40 m.y.), 20–30 m.y. longer than present-day flat slabs. As shown in Figure 9, if a 90-km-deep slab existed for 40 m.y., it could decrease the Moho temperature and surface heat flow by an extra ~130 °C and 5 mW/m² compared to a slab with a 10 m.y. duration. The cooling could also have been enhanced by the old (100 Ma) and cool Farallon plate that had a fast convergence rate (12 cm/yr). In addition, the old Farallon plate may also have allowed slab dehydration to occur farther inboard (Humphreys et al., 2003; Currie and Beaumont, 2011), which would increase the possibility of fluid-induced cooling.

The relatively low surface heat flow observed in continental areas above a flat slab is commonly attributed to cooling of the continent by the presence of the cold oceanic plate at the base of the continental lithosphere (e.g., Dumitru et al., 1991; van Hunen et al., 2004; Currie and Hyndman, 2006; Manea et al., 2017). The numerical models in this study allow an assessment of the evolution of continental thermal structure during flat-slab subduction. Our models show that once a flat slab forms, the continent undergoes conductive cooling over timescales of tens of millions of years. Cool temperatures diffuse from the slab upward through the continental lithosphere, which leads to a delay between the emplacement of the flat slab and the onset of cooling of the shallow continent (Fig. 10). This lag time depends primarily on the depth of the flat slab, with surface cooling initiating at 4 m.y. and 32 m.y. for slabs at 50 km and 120 km depths, respectively. At a given time after flat-slab emplacement, the magnitude of cooling is enhanced by a high thermal conductivity, an old and/or shallow oceanic plate, and a fast convergence rate.

Present-day flat slabs have only been in place for <20 Ma. The Mexico and Peru flat slabs are at shallow depths (45–70 km), and therefore the observed low surface heat flow in these areas is likely due to slab cooling. In contrast, the Pampean flat-slab is at 90–100 km depth, and conductive cooling cannot explain a reduction of surface heat flow in this area. We suggest that additional cooling may arise from advective heat transfer as cool fluids are released through dehydration of the flat-slab and flux through the continental lithosphere. Our models also have implications for ancient flat-slab regions. Widespread cooling in the western United States in the Late Cretaceous is consistent with the presence of a long-lived old Farallon flat-slab with a fast convergence rate and slab dehydration across a wide area.

This work was supported by an NSERC Discovery Grant. Numerical models use the SOPALE finite element code (Fullsack, 1995; http://geodynamics.oceanography.dal.ca/sopaledoc.html) developed by the Dalhousie Geodynamics Group under the direction of Christopher Beaumont, Dalhousie University, Halifax, Canada. All of the model parameters used in this study are available in Table 1. We thank Shanaka de Silva, Robert Stern, and two anonymous reviewers for constructive comments and suggestions on the manuscript.