Most flat-slab subduction regions are marked by an absence of arc volcanism, which is consistent with closure of the hot mantle wedge as the subducting plate flattens below the continent. Farther inland, low surface heat flow is observed, which is generally attributed to cooling of the continent by the underlying flat slab. However, modern flat slabs have only been in place for <20 Ma, and it is unclear whether there has been sufficient time for cooling to occur. We use numerical models to assess temporal variations in continental thermal structure during flat-slab subduction. Our models show that the flat slab leads to continental cooling on timescales of tens of millions of years. Cool slab temperatures must diffuse through the continental lithosphere, resulting in a delay between slab emplacement and surface cooling. Therefore, the timescales primarily depend on the flat-slab depth with shallower slabs resulting in shorter timescales. The magnitude of cooling increases for a shallow or long-lived flat slab, old subducting plate, and fast convergence rates. For regions with flat slabs at 45–70 km depth (e.g., Mexico and Peru), shallow continental cooling initiates 5–10 Ma after slab emplacement, and low surface heat flow in these regions is largely explained by the presence of the flat slab. However, for the Pampean region in Chile, with an ~100-km-deep slab, our models predict that conductive cooling has not yet affected the surface heat flow. The low heat flow observed requires additional processes such as advective cooling from the infiltration of fluids released through dehydration of the flat slab.

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/m2; in contrast, most forearcs have a heat flow of 30–50 mW/m2 (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/m2; 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).

Figure 1.

(A) The Pampean flat-slab area is shown. Large diamonds are heat flow data points from Uyeda and Watanabe (1982), Hamza and Muñoz (1996), and Collo et al. (2018). White circles are the projected Juan Fernandez Ridge track locations following Yáñez et al. (2001). Red triangles show Holocene volcanism (Global Volcanism Program, 2013). Solid black line (X-X’) indicates location of cross section shown in (B) with tick marks every 100 km. Dashed gray lines show terrane boundaries following Ramos et al. (2010) and Martin et al. (2020). Small circles show earthquake locations from the Reviewed International Seismological Centre catalog from 1990 to 2015 for events with depths >50 km and M>4 (International Seismological Centre, 2020). Dotted contours show slab depth at 20 km intervals from Anderson et al. (2007). (B) Cross section along X-X’ profile. Black line at top shows topography/bathymetry. Diamonds on topography show heat flow measurements within 100 km of the profile projected to the closest point on the profile. Red solid line shows the slab surface from Anderson et al. (2007). Dashed red line is inferred slab continuation based on tomography and earthquake locations. Small black dots show the Reviewed International Seismological Centre catalog locations within 20 km of the profile. White circles are projected earthquakes from Anderson et al. (2007) within 20 km of the profile. Colors within dash-dot outlined box are absolute shear wave velocities from the Rayleigh wave tomography of Porter et al. (2012). Background colors are teleseismic P-wave velocity deviations from Portner et al. (2017).

Figure 1.

(A) The Pampean flat-slab area is shown. Large diamonds are heat flow data points from Uyeda and Watanabe (1982), Hamza and Muñoz (1996), and Collo et al. (2018). White circles are the projected Juan Fernandez Ridge track locations following Yáñez et al. (2001). Red triangles show Holocene volcanism (Global Volcanism Program, 2013). Solid black line (X-X’) indicates location of cross section shown in (B) with tick marks every 100 km. Dashed gray lines show terrane boundaries following Ramos et al. (2010) and Martin et al. (2020). Small circles show earthquake locations from the Reviewed International Seismological Centre catalog from 1990 to 2015 for events with depths >50 km and M>4 (International Seismological Centre, 2020). Dotted contours show slab depth at 20 km intervals from Anderson et al. (2007). (B) Cross section along X-X’ profile. Black line at top shows topography/bathymetry. Diamonds on topography show heat flow measurements within 100 km of the profile projected to the closest point on the profile. Red solid line shows the slab surface from Anderson et al. (2007). Dashed red line is inferred slab continuation based on tomography and earthquake locations. Small black dots show the Reviewed International Seismological Centre catalog locations within 20 km of the profile. White circles are projected earthquakes from Anderson et al. (2007) within 20 km of the profile. Colors within dash-dot outlined box are absolute shear wave velocities from the Rayleigh wave tomography of Porter et al. (2012). Background colors are teleseismic P-wave velocity deviations from Portner et al. (2017).

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/m2 (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/m2 (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/m2, and there is little evidence of cooling since the Miocene. This is higher than the heat flow values of <60 mW/m2 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/m2 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/m2, which is lower than that of the back-arc regions above the steep-angle slabs to the north and south (80–120 mW/m2) (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.

Figure 2.

(A) Initial model geometry and boundary conditions are shown. T is temperature and v is velocity. (B) Initial geotherms for the continental (cont.) and the oceanic lithosphere.

Figure 2.

(A) Initial model geometry and boundary conditions are shown. T is temperature and v is velocity. (B) Initial geotherms for the continental (cont.) and the oceanic lithosphere.

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).

TABLE 1.

MODEL PARAMETERS

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/m3 and 0.4 μW/m3, 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/m3 to 3500 kg/m3 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/m3). 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/m2. 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 (Vo in Fig. 2A). The continent has a trenchward motion of 3 cm/yr (Vc), 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 (Vb) is assigned to the side boundaries of the sublithospheric mantle. Models are solved in the continental reference frame, and therefore Vc 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.

Figure 3.

Reference model evolution from 0 m.y. to 20 m.y. with surface elevation. Black arrows show the velocity field of the model. The material colors are shown in the legend and follow those in Figure 2; oceanic crust eclogite is marked by the dark red color.

Figure 3.

Reference model evolution from 0 m.y. to 20 m.y. with surface elevation. Black arrows show the velocity field of the model. The material colors are shown in the legend and follow those in Figure 2; oceanic crust eclogite is marked by the dark red color.

TABLE 2.

LIST OF MODELS AND THE CHANGE IN MOHO TEMPERATURE AND SURFACE HEAT FLOW

Reference Model

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/m2, which is consistent with observations in forearc regions (e.g., Wada and Wang, 2009).

Figure 4.

Thermal evolution of the reference model. (A) Model temperature field at given times with Moho temperature (top plot) and surface heat flow (second plot). Gray lines indicate the initial values for the continental region. In the model plots, white lines are isotherms every 300 °C and red lines are the oceanic Moho. Blue, red, and yellow vertical dotted lines are at 400 km, 600 km, and 800 km. (B) Vertical temperature profile 400 km from the trench; the red line is the initial continental geotherm. (C) Surface geothermal gradient at 400 km from the trench. (D) Moho temperature and (E) surface heat flow over time at three positions from the trench (400 km, 600 km, and 800 km).

Figure 4.

Thermal evolution of the reference model. (A) Model temperature field at given times with Moho temperature (top plot) and surface heat flow (second plot). Gray lines indicate the initial values for the continental region. In the model plots, white lines are isotherms every 300 °C and red lines are the oceanic Moho. Blue, red, and yellow vertical dotted lines are at 400 km, 600 km, and 800 km. (B) Vertical temperature profile 400 km from the trench; the red line is the initial continental geotherm. (C) Surface geothermal gradient at 400 km from the trench. (D) Moho temperature and (E) surface heat flow over time at three positions from the trench (400 km, 600 km, and 800 km).

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/m2 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/m2, and the decreases of surface heat flow are 1.5 mW/m2, 0.4 mW/m2, and 0.2 mW/m2 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).

Effect of Subduction Angle

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/m2 heat flow).

Figure 5.

Normal-angle subduction model. (A) Model temperature field evolution with surface topography, surface heat flow, and Moho temperature; white lines are isotherms every 300 °C; red lines are the oceanic Moho. Temporal changes in (B) Moho temperature and (C) surface heat flow of the normal subduction model (blue line) and flat-slab subduction model (orange line; also see Fig. 4A) at 400 km from the trench.

Figure 5.

Normal-angle subduction model. (A) Model temperature field evolution with surface topography, surface heat flow, and Moho temperature; white lines are isotherms every 300 °C; red lines are the oceanic Moho. Temporal changes in (B) Moho temperature and (C) surface heat flow of the normal subduction model (blue line) and flat-slab subduction model (orange line; also see Fig. 4A) at 400 km from the trench.

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/m2. This shows that the continental cooling in the reference model results from the presence of the flat slab.

Variations in Thickness of Initial Continental Lithosphere and Flat-Slab Depth

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/m2, 57.4 mW/m2, and 54.1 mW/m2 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.

Figure 6.

(A) Temperature fields at 20 m.y. for models with continents that are 80 km, 120 km, and 160 km thick. Note that the flat-slab depth increases with increasing initial continental thickness. The plots on the right show the temporal changes in Moho temperature and surface heat flow of the corresponding model. (B) Geotherm profiles from 0 m.y. to 30 m.y. with an interval of 1 m.y.; red line, green line, and blue lines represent models with continents that are 80 km, 120 km, and 160 km thick, respectively. Temporal changes in (C) Moho temperature and (D) surface heat flow relative to the initial values for the three models at 400 km from the trench.

Figure 6.

(A) Temperature fields at 20 m.y. for models with continents that are 80 km, 120 km, and 160 km thick. Note that the flat-slab depth increases with increasing initial continental thickness. The plots on the right show the temporal changes in Moho temperature and surface heat flow of the corresponding model. (B) Geotherm profiles from 0 m.y. to 30 m.y. with an interval of 1 m.y.; red line, green line, and blue lines represent models with continents that are 80 km, 120 km, and 160 km thick, respectively. Temporal changes in (C) Moho temperature and (D) surface heat flow relative to the initial values for the three models at 400 km from the trench.

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/m2 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/m2 and 0.5 mW/m2 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).

Effects of Other Factors

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

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 7.

Numerical model results with different thermal conductivities are plotted. The green, blue, red, and yellow lines represent models with conductivities of 2.25 W/(m·K) (reference model), 3.25 W/(m·K), 4.25 W/(m·K), and 5.25 W/(m·K), respectively. (A) Moho temperature and Moho temperature change relative to the initial value. (B) Surface heat flow and heat flow change relative to the initial value.

Figure 7.

Numerical model results with different thermal conductivities are plotted. The green, blue, red, and yellow lines represent models with conductivities of 2.25 W/(m·K) (reference model), 3.25 W/(m·K), 4.25 W/(m·K), and 5.25 W/(m·K), respectively. (A) Moho temperature and Moho temperature change relative to the initial value. (B) Surface heat flow and heat flow change relative to the initial value.

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/m2, 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.

Continental Crust Thickness

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/m3 and 0.4 μW/m3 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/m2), and therefore we do not show it here.

Figure 8.

(A) Moho temperature and temperature change relative to the initial value for models with continental crustal thicknesses of 36 km, 42 km, and 48 km. (B) Moho temperature of models with fast (10 cm/yr) and slow (7 cm/yr) convergence rates. (C) Moho temperature of models with old (100 Ma) and young (30 Ma) oceanic plates.

Figure 8.

(A) Moho temperature and temperature change relative to the initial value for models with continental crustal thicknesses of 36 km, 42 km, and 48 km. (B) Moho temperature of models with fast (10 cm/yr) and slow (7 cm/yr) convergence rates. (C) Moho temperature of models with old (100 Ma) and young (30 Ma) oceanic plates.

Convergence Rate

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/m2. These changes are lower than those in the reference model (77 °C and 1.5 mW/m2), 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.

Oceanic Plate Age

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.

Summary of Numerical Models

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. 45), 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.

Limitations of Numerical Models

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.

Timescale for Continental Cooling

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/m2 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/m2 and 479 °C, respectively.

Figure 9.

Time-dependent, 1-D temperature calculations (calc.) are shown for a 120-km-thick continent that is cooled by a flat slab at a depth of 90 km. (A) Vertical temperature profiles after the emplacement of a 600 °C flat slab from 0 m.y. to 60 m.y.; profiles are taken every 5 m.y. The initial continental geotherm is shown in green. (B) Surface heat flow and (C) Moho temperature of the 1-D model as a function of time after the emplacement of a flat slab with different temperatures (blue dotted line: 400 °C; black line: 600 °C; green dotted line: 800 °C). The red line shows the numerical model results for the reference model at 400 km from the trench. The blue and red vertical lines are the lag times and time constants, respectively.

Figure 9.

Time-dependent, 1-D temperature calculations (calc.) are shown for a 120-km-thick continent that is cooled by a flat slab at a depth of 90 km. (A) Vertical temperature profiles after the emplacement of a 600 °C flat slab from 0 m.y. to 60 m.y.; profiles are taken every 5 m.y. The initial continental geotherm is shown in green. (B) Surface heat flow and (C) Moho temperature of the 1-D model as a function of time after the emplacement of a flat slab with different temperatures (blue dotted line: 400 °C; black line: 600 °C; green dotted line: 800 °C). The red line shows the numerical model results for the reference model at 400 km from the trench. The blue and red vertical lines are the lag times and time constants, 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.

The critical times indicate when the thermal effects of the flat slab should be first observed (τl) and the timeframe needed for significant cooling to occur (τc). We can quantitatively determine τl and τc based on the variations in thermal structure relative to the initial and final conditions. The relative change (R) in thermal structure as a function of time is given by:
formula
where X is the measured thermal parameter (Moho temperature or surface heat flow), X0 is the initial value, and XS is the final (steady-state) value. The lag time (τl) is defined as the time at which the temperature field has changed by e% (e = Euler's number):
formula
Following Turcotte and Schubert (2002), the thermal time constant (τc) is the time at which:
formula

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/m2 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).

Figure 10.

(A) Moho temperature and (B) surface heat flow contours are shown for a 120-km-thick continent and 600 °C flat slab with various flat-slab depths and slab emplacement times. The blue line shows the lag time (τl), which is the time at which slab cooling initiates. The red line shows the time constant (τc) for cooling. The red star represents the Pampean flat-slab subduction condition. The plots on the right show the steady-state values of (A) Moho temperature and (B) surface heat flow.

Figure 10.

(A) Moho temperature and (B) surface heat flow contours are shown for a 120-km-thick continent and 600 °C flat slab with various flat-slab depths and slab emplacement times. The blue line shows the lag time (τl), which is the time at which slab cooling initiates. The red line shows the time constant (τc) for cooling. The red star represents the Pampean flat-slab subduction condition. The plots on the right show the steady-state values of (A) Moho temperature and (B) surface heat flow.

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/m2. 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/m2 within 30 m.y.

Figure 11.

Comparison of (A) Moho temperatures and (B) surface heat flows between the 120-km-thick continent (cont.) and the 60-km-thick continent for a 600 °C flat slab at 50 km depth. The blue and red vertical lines are the lag times and time constants, respectively.

Figure 11.

Comparison of (A) Moho temperatures and (B) surface heat flows between the 120-km-thick continent (cont.) and the 60-km-thick continent for a 600 °C flat slab at 50 km depth. The blue and red vertical lines are the lag times and time constants, respectively.

Implications for Modern and Ancient Flat-Slab Regions

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/m2; 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/m2 within 20 m.y. through thermal conductive heat transfer (Fig. 11B). However, 20 mW/ m2 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/m2; 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/m2, which is ~20 mW/m2 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/m2 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.

Science Editor: Shanaka de Silva
Guest Associate Editor: Robert J. Stern
1.
Abers
,
G.A.
, and
Hacker
,
B.R.
,
2016
,
A MATLAB toolbox and Excel workbook for calculating the densities, seismic wave speeds, and major element composition of minerals and rocks at pressure and temperature
:
Geochemistry, Geophysics, Geosystems
 , v.
17
, no.
2
, p.
616
624
. https://doi.org/10.1002/2015GC006171.
2.
Alvarado
,
P.
,
Pardo
,
M.
,
Gilbert
,
H.
,
Miranda
,
S.
,
Anderson
,
M.
,
Saez
,
M.
, and
Beck
,
S.
,
2009
,
Flat-slab subduction and crustal models for the seismically active Sierras Pampeanas region of Argentina
, in
Kay
,
S.M.
,
Ramos
,
V.A.
, and
Dickinson
,
W.R.
, eds.,
Backbone of the Americas: Shallow Subduction, Plateau Uplift, and Ridge and Terrane Collision: Geological Society of America Memoir 204
 , p.
261
278
. https://doi.org/10.1130/2009.1204(12).
3.
Anderson
,
M.
,
Alvarado
,
P.
,
Zandt
,
G.
, and
Beck
,
S.
,
2007
,
Geometry and brittle deformation of the subducting Nazca Plate, Central Chile and Argentina
:
Geophysical Journal International
 , v.
171
, no.
1
, p.
419
434
. https://doi.org/10.1111/j.1365-246X.2007.03483.x.
4.
Antonijevic
,
S.K.
,
Wagner
,
L.S.
,
Kumar
,
A.
,
Beck
,
S.L.
,
Long
,
M.D.
,
Zandt
,
G.
,
Tavera
,
H.
, and
Condori
,
C.
,
2015
,
The role of ridges in the formation and longevity of flat slabs
:
Nature
 , v.
524
, no.
7564
, p.
212
215
. https://doi.org/10.1038/nature14648.
5.
Axen
,
G.J.
,
van Wijk
,
J.W.
, and
Currie
,
C.A.
,
2018
,
Basal continental mantle lithosphere displaced by flat-slab subduction
:
Nature Geoscience
 , v.
11
, https://doi.org/10.1038/s41561-018-0263-9.
6.
Beaumont
,
C.
,
Nguyen
,
M.H.
,
Jamieson
,
R.A.
, and
Ellis
,
S.
,
2006
,
Crustal flow modes in large hot orogens
, in
Law
,
R.D.
,
Searle
,
M.P.
, and
Godin
,
L.
, eds.,
Channel Flow, Ductile Extrusion and Exhumation in Continental Collision Zones
 :
Geological Society
,
London, Special Publication 268
, no.
1
, p.
91
145
. https://doi.org/10.1144/gsl.sp.2006.268.01.05.
7.
Bird
,
P.
,
1988
,
Formation of the Rocky Mountains, western United States: A continuum computer model
:
Science
 , v.
239
, no.
4847
, p.
1501
1507
. https://doi.org/10.1126/science.239.4847.1501.
8.
Bishop
,
B.T.
,
Beck
,
S.L.
,
Zandt
,
G.
,
Wagner
,
L.
,
Long
,
M.
,
Antonijevic
,
S.K.
,
Kumar
,
A.
, and
Tavera
,
H.
,
2017
,
Causes and consequences of flat-slab subduction in southern Peru
:
Geosphere
 , v.
13
, no.
5
, p.
1392
1407
. https://doi.org/10.1130/GES01440.1.
9.
Bloch
,
W.
,
John
,
T.
,
Kummerow
,
J.
,
Salazar
,
P.
,
Krüger
,
O.S.
, and
Shapiro
,
S.A.
,
2018
,
Watching dehydration: Seismic indication for transient fluid pathways in the oceanic mantle of the subducting Nazca slab
:
Geochemistry, Geophysics, Geosystems
 , v.
19
, no.
9
, p.
3189
3207
. https://doi.org/10.1029/2018GC007703.
10.
Burns
,
E.R.
,
Williams
,
C.F.
,
Ingebritsen
,
S.E.
,
Voss
,
C.I.
,
Spane
,
F.A.
, and
DeAngelo
,
J.
,
2015
,
Understanding heat and groundwater flow through continental flood basalt provinces: Insights gained from alternative models of permeability/depth relationships for the
Columbia Plateau, USA
:
Geofluids
, v.
15
, no.
1–2
, p.
120
138
. https://doi.org/10.1111/gfl.12095.
11.
Cai
,
C.
,
Wiens
,
D.A.
,
Shen
,
W.
, and
Eimer
,
M.
,
2018
,
Water input into the Mariana subduction zone estimated from ocean-bottom seismic data
:
Nature
 , v.
563
, p.
389
392
. https://doi.org/10.1038/s41586-018-0655-4.
12.
Cerpa
,
N.G.
,
Wada
,
I.
, and
Wilson
,
C.R.
,
2017
,
Fluid migration in the mantle wedge: Influence of mineral grain size and mantle compaction
:
Journal of Geophysical Research: Solid Earth
 , v.
122
, no.
8
, p.
6247
6268
. https://doi.org/10.1002/2017JB014046.
13.
Clauser
,
C.
, and
Huenges
,
E.
,
1995
,
Thermal conductivity of rocks and minerals
, in
Ahrens
,
T.J.
, ed.,
Rock Physics and Phase Relations: A Handbook of Physical Constants
 :
Washington, D.C.
,
American Geophysical Union
,
Reference Shelf Series
 , v.
3
, p.
105
126
. https://doi.org/10.1029/rf003p0105.
14.
Collo
,
G.
,
Dávila
,
F.M.
,
Nóbile
,
J.
,
Astini
,
R.A.
, and
Gehrels
,
G.
,
2011
,
Clay mineralogy and thermal history of the Neogene Vinchina Basin, central Andes of Argentina: Analysis of factors controlling the heating conditions
:
Tectonics
 , v.
30
, no.
4
, https://doi.org/10.1029/2010TC002841.
15.
Collo
,
G.
,
Dávila
,
F.M.
,
Teixeira
,
W.
,
Nóbile
,
J.C.
,
Sant
,
L.G.
, and
Carter
,
A.
,
2017
,
Isotopic and thermochronologic evidence of extremely cold lithosphere associated with a slab flattening in the Central Andes of Argentina
:
Basin Research
 , v.
29
, no.
1
, p.
16
40
. https://doi.org/10.1111/bre.12163.
16.
Collo
,
G.
,
Ezpeleta
,
M.
,
Dávila
,
F.M.
,
Giménez
,
M.
,
Soler
,
S.
,
Martina
,
F.
,
Ávila
,
P.
,
Sánchez
,
F.
,
Calegari
,
R.
,
Lovecchio
,
J.
, and
Schiuma
,
M.
,
2018
,
Basin thermal structure in the Chilean-Pampean flat subduction zone
, in
Folguera
,
A.
, et al., eds.,
The Evolution of the Chilean-Argentinean Andes
 :
Cham
,
Springer
, p.
537
564
. https://doi.org/10.1007/978-3-319-67774-3_21.
17.
Coney
,
P.J.
, and
Reynolds
,
S.J.
,
1977
,
Cordilleran Benioff zones
:
Nature
 , v.
270
, p.
403
406
. https://doi.org/10.1038/270403a0.
18.
Copeland
,
P.
,
Currie
,
C.A.
, and
Lawton
,
T.
,
2017
,
Location, location, location: The variable lifespan of the Laramide orogeny
:
Geology
 , v.
45
, p.
223
226
. https://doi.org/10.1130/G38810.1.
19.
Currie
,
C.A.
, and
Beaumont
,
C.
,
2011
,
Are diamond-bearing Cretaceous kimberlites related to low-angle subduction beneath western North America?
:
Earth and Planetary Science Letters
 , v.
303
, no.
1–2
, p.
59
70
. https://doi.org/10.1016/j.epsl.2010.12.036.
20.
Currie
,
C.A.
, and
Hyndman
,
R.D.
,
2006
,
The thermal structure of subduction zone back arcs
:
Journal of Geophysical Research: Solid Earth
 , v.
111
, https://doi.org/10.1029/2005jb004024.
21.
Dumitru
,
T.A.
,
1990
,
Subnormal Cenozoic geothermal gradients in the extinct Sierra Nevada magmatic arc: Consequences of Laramide and post‐Laramide shallow‐angle subduction
:
Journal of Geophysical Research: Solid Earth
 , v.
95
, p.
4925
4941
. https://doi.org/10.1029/JB095iB04p04925.
22.
Dumitru
,
T.A.
,
Gans
,
P.B.
,
Foster
,
D.A.
, and
Miller
,
E.L.
,
1991
,
Refrigeration of the western Cordilleran lithosphere during Laramide shallow-angle subduction
:
Geology
 , v.
19
, no.
11
, p.
1145
1148
. https://doi.org/10.1130/0091-7613(1991)019<1145:ROTWCL>2.3.CO;2.
23.
Ferrari
,
L.
,
2004
,
Slab detachment control on mafic volcanic pulse and mantle heterogeneity in central Mexico
:
Geology
 , v.
32
, no.
1
, p.
77
80
. https://doi.org/10.1130/G19887.1.
24.
Ferrari
,
L.
,
López-Martínez
,
M.
,
Aguirre-Díaz
,
G.
, and
Carrasco-Núñez
,
G.
,
1999
,
Space-time patterns of Cenozoic arc volcanism in central Mexico: From the Sierra Madre Occidental to the Mexican Volcanic Belt
:
Geology
 , v.
27
, no.
4
, p.
303
306
. https://doi.org/10.1130/0091-7613(1999)027<0303:STPOCA>2.3.CO;2.
25.
Ferrand
,
T.P.
,
Hilairet
,
N.
,
Incel
,
S.
,
Deldicque
,
D.
,
Labrousse
,
L.
,
Gasc
,
J.
,
Renner
,
J.
,
Wang
,
Y.B.
,
Green
,
H.W.
, II
, and
Schubnel
,
A.
,
2017
,
Dehydration-driven stress transfer triggers intermediate-depth earthquakes
:
Nature Communications
 , v.
8
, p.
1
11
. https://doi.org/10.1038/ncomms15247.
26.
Fromm
,
R.
,
Alvarado
,
P.
,
Beck
,
S.L.
, and
Zandt
,
G.
,
2006
,
The April 9, 2001 Juan Fernández Ridge (Mw 6.7) tensional outer-rise earthquake and its aftershock sequence
:
Journal of Seismology
 , v.
10
, no.
2
, p.
163
170
. https://doi.org/10.1007/s10950-006-9013-3.
27.
Fullsack
,
P.
,
1995
,
An arbitrary Lagrangian‐Eulerian formulation for creeping flows and its application in tectonic models
:
Geophysical Journal International
 , v.
120
, no.
1
, p.
1
23
. https://doi.org/10.1111/j.1365-246X.1995.tb05908.x.
28.
Gans
,
C.R.
,
Beck
,
S.L.
,
Zandt
,
G.
,
Gilbert
,
H.
,
Alvarado
,
P.
,
Anderson
,
M.
, and
Linkimer
,
L.
,
2011
,
Continental and oceanic crustal structure of the Pampean flat slab region, western Argentina, using receiver function analysis: New high-resolution results
:
Geophysical Journal International
 , v.
186
, no.
1
, p.
45
58
. https://doi.org/10.1111/j.1365-246X.2011.05023.x.
29.
Gerya
,
T.
,
2019
,
Numerical solution of the heat conservation equation, in Introduction to Numerical Geodynamic Modelling (second edition)
:
Cambridge, UK
,
Cambridge University Press
, p.
139
155
. https://doi.org/10.1017/9781316534243.
30.
Gleason
,
G.C.
, and
Tullis
,
J.
,
1995
,
A flow law for dislocation creep of quartz aggregates determined with the molten salt cell
:
Tectonophysics
 , v.
247
, no.
1–4
, p.
1
23
. https://doi.org/10.1016/0040-1951(95)00011-B.
31.
Global Volcanism
Program
,
2013
,
Volcanoes of the World
, in
Venzke
,
E.
, eds.,
Smithsonian Institution
 , v
4.10.3
,https://doi.org/10.5479/si.GVP.VOTW4-2013
(accessed 11 November 2019).
32.
Grevemeyer
,
I.
,
Ranero
,
C.R.
, and
Ivandic
,
M.
,
2018
,
Structure of oceanic crust and serpentinization at subduction trenches
:
Geosphere
 , v.
14
, no.
2
, p.
395
418
. https://doi.org/10.1130/GES01537.1.
33.
Gutscher
,
M.A.
,
2002
,
Andean subduction styles and their effect on thermal structure and interplate coupling
:
Journal of South American Earth Sciences
 , v.
15
, no.
1
, p.
3
10
. https://doi.org/10.1016/S0895-9811(02)00002-0.
34.
Gutscher
,
M.A.
,
Spakman
,
W.
,
Bijwaard
,
H.
, and
Engdahl
,
E.R.
,
2000
,
Geodynamics of flat subduction: Seismicity and tomographic constraints from the Andean margin
:
Tectonics
 , v.
19
, no.
5
, p.
814
833
. https://doi.org/10.1029/1999TC001152.
35.
Hacker
,
B.R.
, and
Abers
,
G.A.
,
2004
,
Subduction Factory 3: An Excel worksheet and macro for calculating the densities, seismic wave speeds, and H2O contents of minerals and rocks at pressure and temperature
:
Geochemistry, Geophysics, Geosystems
 , v.
5
, no.
1
, https://doi.org/10.1029/2003GC000614.
36.
Hacker
,
B.R.
,
Abers
,
G.A.
, and
Peacock
,
S.M.
,
2003
,
Subduction factory 1. Theoretical mineralogy, densities, seismic wave speeds, and H2O contents
:
Journal of Geophysical Research: Solid Earth
 , v.
108
, https://doi.org/10.1029/2001JB001127.
37.
Hampel
,
A.
,
2002
,
The migration history of the Nazca Ridge along the Peruvian active margin: A re-evaluation
:
Earth and Planetary Science Letters
 , v.
203
, no.
2
, p.
665
679
. https://doi.org/10.1016/S0012-821X(02)00859-2.
38.
Hamza
,
V.M.
, and
Muñoz
,
M.
,
1996
,
Heat flow map of South America
:
Geothermics
 , v.
25
, no.
6
, p.
599
646
. https://doi.org/10.1016/S0375-6505(96)00025-9.
39.
Henry
,
S.G.
, and
Pollack
,
H.N.
,
1988
,
Terrestrial heat flow above the Andean subduction zone in Bolivia and Peru
:
Journal of Geophysical Research: Solid Earth
 , v.
93
, p.
15153
15162
. https://doi.org/10.1029/JB093iB12p15153.
40.
Hu
,
J.
, and
Liu
,
L.
,
2016
,
Abnormal seismological and magmatic processes controlled by the tearing South American flat slabs
:
Earth and Planetary Science Letters
 , v.
450
, p.
40
51
. https://doi.org/10.1016/j.epsl.2016.06.019.
41.
Humphreys
,
E.
,
Hessler
,
E.
,
Dueker
,
K.
,
Farmer
,
G.L.
,
Erslev
,
E.
, and
Atwater
,
T.
,
2003
,
How Laramide-age hydration of North American lithosphere by the Farallon slab controlled subsequent activity in the western United States
:
International Geology Review
 , v.
45
, no.
7
, p.
575
595
. https://doi.org/10.2747/0020-6814.45.7.575.
42.
Hyndman
,
R.D.
,
Currie
,
C.A.
, and
Mazzotti
,
S.P.
,
2005
,
Subduction zone backarcs, mobile belts, and orogenic heat
:
GSA Today
 , v.
15
, no.
2
, p.
4
10
. https://doi.org/10.1130/1052-5173(2005)015<4:SZBMBA>2.0.CO;2.
43.
International Seismological Centre
,
2020
,
On-line Bulletin
 , https://doi.org/10.31905/D808B830.
44.
Jaupart
,
C.
, and
Mareschal
,
J.-C.
,
2014
,
Constraints on crustal heat production from heat flow data
, in
Holland
,
H.D.
, and
Turekian
,
K.K.
, eds.,
Treatise on Geochemistry (second edition): Reference Module in Earth Systems and Environmental Sciences
 , v.
4
, p.
53
73
. https://doi.org/10.1016/b978-0-08-095975-7.00302-8.
45.
Karato
,
S.I.
, and
Wu
,
P.
,
1993
,
Rheology of the upper mantle: A synthesis
:
Science
 , v.
260
, no.
5109
, p.
771
778
. https://doi.org/10.1126/science.260.5109.771.
46.
Kay
,
S.M.
, and
Abbruzzi
,
J.M.
,
1996
,
Magmatic evidence for Neogene lithospheric evolution of the central Andean “flat-slab” between 30°S and 32°S
:
Tectonophysics
 , v.
259
, p.
15
28
. https://doi.org/10.1016/0040-1951(96)00032-7.
47.
Kay
,
S.M.
, and
Mpodozis
,
C.
,
2002
,
Magmatism as a probe to the Neogene shallowing of the Nazca plate beneath the modern Chilean flat-slab
:
Journal of South American Earth Sciences
 , v.
15
, no.
1
, p.
39
57
. https://doi.org/10.1016/S0895-9811(02)00005-6.
48.
Kim
,
Y.
,
Miller
,
M.S.
,
Pearce
,
F.
, and
Clayton
,
R.W.
,
2012
,
Seismic imaging of the Cocos plate subduction zone system in central Mexico
:
Geochemistry, Geophysics, Geosystems
 , v.
13
, no.
7
, https://doi.org/10.1029/2012GC004033.
49.
Kopp
,
H.
,
Flueh
,
E.R.
,
Papenberg
,
C.
, and
Klaeschen
,
D.
,
2004
,
Seismic investigations of the O’Higgins Seamount Group and Juan Fernández Ridge: Aseismic ridge emplacement and lithosphere hydration
:
Tectonics
 , v.
23
, no.
2
, https://doi.org/10.1029/2003TC001590.
50.
Kumar
,
A.
,
Wagner
,
L.S.
,
Beck
,
S.L.
,
Long
,
M.D.
,
Zandt
,
G.
,
Young
,
B.
,
Tavera
,
H.
, and
Minaya
,
E.
,
2016
,
Seismicity and state of stress in the central and southern Peruvian flat slab
:
Earth and Planetary Science Letters
 , v.
441
, p.
71
80
. https://doi.org/10.1016/j.epsl.2016.02.023.
51.
Lewis
,
T.J.
,
Hyndman
,
R.D.
, and
Flück
,
P.
,
2003
,
Heat flow, heat generation, and crustal temperatures in the northern Canadian Cordillera: Thermal control of tectonics
:
Journal of Geophysical Research: Solid Earth
 , v.
108
, https://doi.org/10.1029/2002jb002090.
52.
Li
,
Z.X.
, and
Li
,
X.H.
,
2007
,
Formation of the 1300-km-wide intracontinental orogen and postorogenic magmatic province in Mesozoic South China: A flat-slab subduction model
:
Geology
 , v.
35
, no.
2
, p.
179
182
. https://doi.org/10.1130/G23193A.1.
53.
Lipman
,
P.W.
,
Prostka
,
H.J.
, and
Christiansen
,
R.L.
,
1972
,
A discussion on volcanism and the structure of the Earth-Cenozoic volcanism and plate-tectonic evolution of the Western United States. I. Early and middle Cenozoic
:
Philosophical Transactions of the Royal Society A: Mathematical, Physical, and Engineering Sciences
 , v.
271
, no.
1213
, p.
217
248
. https://doi.org/10.1098/rsta.1972.0008.
54.
Liu
,
S.
, and
Currie
,
C.A.
,
2016
,
Farallon plate dynamics prior to the Laramide orogeny: Numerical models of flat subduction
:
Tectonophysics
 , v.
666
, p.
33
47
. https://doi.org/10.1016/j.tecto.2015.10.010.
55.
Liu
,
X.
, and
Currie
,
C.A.
,
2019
,
Influence of upper plate structure on flat‐slab depth: Numerical modeling of subduction dynamics
:
Journal of Geophysical Research: Solid Earth
 , v.
124
, no.
12
, p.
13150
13167
. https://doi.org/10.1029/2019JB018653.
56.
Löbens
,
S.
,
Bense
,
F.A.
,
Wemmer
,
K.
,
Dunkl
,
I.
,
Costa
,
C.H.
,
Layer
,
P.
, and
Siegesmund
,
S.
,
2011
,
Exhumation and uplift of the Sierras Pampeanas: Preliminary implications from K–Ar fault gouge dating and low-T thermochronology in the Sierra de Comechingones (Argentina)
:
International Journal of Earth Sciences
 , v.
100
, no.
2–3
, p.
671
694
. https://doi.org/10.1007/s00531-010-0608-0.
57.
Mackwell
,
S.J.
,
Zimmerman
,
M.E.
, and
Kohlstedt
,
D.L.
,
1998
,
High‐temperature deformation of dry diabase with application to tectonics on Venus
:
Journal of Geophysical Research: Solid Earth
 , v.
103
, p.
975
984
. https://doi.org/10.1029/97JB02671.
58.
Manea
,
V.C.
,
Pérez-Gussinyé
,
M.
, and
Manea
,
M.
,
2012
,
Chilean flat slab subduction controlled by overriding plate thickness and trench rollback
:
Geology
 , v.
40
, no.
1
, p.
35
38
. https://doi.org/10.1130/G32543.1.
59.
Manea
,
V.C.
,
Manea
,
M.
,
Ferrari
,
L.
,
Orozco-Esquivel
,
T.
,
Valenzuela
,
R.W.
,
Husker
,
A.
, and
Kostoglodov
,
V.
,
2017
,
A review of the geodynamic evolution of flat slab subduction in Mexico, Peru, and Chile
:
Tectonophysics
 , v.
695
, p.
27
52
. https://doi.org/10.1016/j.tecto.2016.11.037.
60.
Marot
,
M.
,
Monfret
,
T.
,
Gerbault
,
M.
,
Nolet
,
G.
,
Ranalli
,
G.
, and
Pardo
,
M.
,
2014
,
Flat versus normal subduction zones: A comparison based on 3-D regional traveltime tomography and petrological modelling of central Chile and western Argentina (29°–35°S)
:
Geophysical Journal International
 , v.
199
, no.
3
, p.
1633
1654
. https://doi.org/10.1093/gji/ggu355.
61.
Martin
,
E.L.
,
Collins
,
W.J.
, and
Spencer
,
C.J.
,
2020
,
Laurentian origin of the Cuyania suspect terrane, western Argentina, confirmed by Hf isotopes in zircon
:
Geological Society of America Bulletin
 , v.
132
, no.
1–2
, p.
273
290
. https://doi.org/10.1130/B35150.1.
62.
Martinod
,
J.
,
Funiciello
,
F.
,
Faccenna
,
C.
,
Labanieh
,
S.
, and
Regard
,
V.
,
2005
,
Dynamical effects of subducting ridges: Insights from 3-D laboratory models
:
Geophysical Journal International
 , v.
163
, no.
3
, p.
1137
1150
. https://doi.org/10.1111/j.1365-246X.2005.02797.x.
63.
Meade
,
C.
, and
Jeanloz
,
R.
,
1991
,
Deep-focus earthquakes and recycling of water into the Earth’s mantle
:
Science
 , v.
252
, no.
5002
, p.
68
72
. https://doi.org/10.1126/science.252.5002.68.
64.
Miller
,
H.
, and
Söllner
,
F.
,
2005
,
The Famatina complex (NW Argentina): Back-docking of an island arc or terrane accretion? Early Palaeozoic geodynamics at the western Gondwana margin
, in
Vaughan
,
A.P.M.
,
Leat
,
P.T.
, and
Pankhurst
,
R.J.
, eds.,
Terrane Processes at the Margins of Gondwana: Geological Society
 ,
London
,
Special Publication 246
, no.
1
, p.
241
256
. https://doi.org/10.1144/GSL.SP.2005.246.01.08.
65.
Müller
,
R.D.
,
Sdrolias
,
M.
,
Gaina
,
C.
, and
Roest
,
W.R.
,
2008
,
Age, spreading rates, and spreading asymmetry of the world’s ocean crust
:
Geochemistry, Geophysics, Geosystems
 , v.
9
, no.
4
, https://doi.org/10.1029/2007GC001743.
66.
Muñoz
,
M.
,
2005
,
No flat Wadati–Benioff Zone in the central and southern central Andes
:
Tectonophysics
 , v.
395
, no.
1–2
, p.
41
65
. https://doi.org/10.1016/j.tecto.2004.09.002.
67.
Norabuena
,
E.O.
,
Dixon
,
T.H.
,
Stein
,
S.
, and
Harrison
,
C.G.
,
1999
,
Decelerating Nazca‐South America and Nazca‐Pacific plate motions
:
Geophysical Research Letters
 , v.
26
, no.
22
, p.
3405
3408
. https://doi.org/10.1029/1999GL005394.
68.
O’Driscoll
,
L.J.
,
Humphreys
,
E.D.
, and
Saucier
,
F.
,
2009
,
Subduction adjacent to deep continental roots: Enhanced negative pressure in the mantle wedge, mountain building and continental motion
:
Earth and Planetary Science Letters
 , v.
280
, no.
1–4
, p.
61
70
. https://doi.org/10.1016/j.epsl.2009.01.020.
69.
Oyhantçabal
,
P.
,
Siegesmund
,
S.
, and
Wemmer
,
K.
,
2011
,
The Río de la Plata Craton: A review of units, boundaries, ages and isotopic signature
:
International Journal of Earth Sciences
 , v.
100
, no.
2–3
, p.
201
220
. https://doi.org/10.1007/s00531-010-0580-8.
70.
Pankhurst
,
R.J.
,
Rapela
,
C.W.
,
Saavedra
,
J.
,
Baldo
,
E.
,
Dahlquist
,
J.
,
Pascua
,
I.
, and
Fanning
,
C.M.
,
1998
,
The Famatinian magmatic arc in the central Sierras Pampeanas: An Early to Mid-Ordovician continental arc on the Gondwana margin
, in
Pankhurst
,
R.J.
, and
Rapela
,
C.W.
, eds.,
The Proto-Andean Margin of Gondwana: Geological Society
 ,
London
,
Special Publication 142
, no.
1
, p.
343
367
. https://doi.org/10.1144/GSL.SP.1998.142.01.17.
71.
Peacock
,
S.M.
, and
Hyndman
,
R.D.
,
1999
,
Hydrous minerals in the mantle wedge and the maximum depth of subduction thrust earthquakes
:
Geophysical Research Letters
 , v.
26
, no.
16
, p.
2517
2520
. https://doi.org/10.1029/1999GL900558.
72.
Pérez‐Campos
,
X.
,
Kim
,
Y.
,
Husker
,
A.
,
Davis
,
P.M.
,
Clayton
,
R.W.
,
Iglesias
,
A.
,
Pacheco
,
J.F.
,
Singh
,
S.K.
,
Manea
,
V.C.
, and
Gurnis
,
M.
,
2008
,
Horizontal subduction and truncation of the Cocos Plate beneath central Mexico
:
Geophysical Research Letters
 , v.
35
, no.
18
, https://doi.org/10.1029/2008GL035127.
73.
Perry
,
M.
,
Spinelli
,
G.A.
,
Wada
,
I.
, and
He
,
J.
,
2016
,
Modeled temperatures and fluid source distributions for the Mexican subduction zone: Effects of hydrothermal circulation and implications for plate boundary seismic processes
:
Geochemistry, Geophysics, Geosystems
 , v.
17
, no.
2
, p.
550
570
. https://doi.org/10.1002/2015GC006148.
74.
Porter
,
R.
,
Gilbert
,
H.
,
Zandt
,
G.
,
Beck
,
S.
,
Warren
,
L.
,
Calkins
,
J.
,
Alvarado
,
P.
, and
Anderson
,
M.
,
2012
,
Shear wave velocities in the Pampean flat‐slab region from Rayleigh wave tomography: Implications for slab and upper mantle hydration
:
Journal of Geophysical Research: Solid Earth
 , v.
117
, https://doi.org/10.1029/2012jb009350.
75.
Portner
,
D.E.
,
Beck
,
S.
,
Zandt
,
G.
, and
Scire
,
A.
,
2017
,
The nature of subslab slow velocity anomalies beneath South America
:
Geophysical Research Letters
 , v.
44
, no.
10
, p.
4747
4755
. https://doi.org/10.1002/2017GL073106.
76.
Pysklywec
,
R.N.
, and
Beaumont
,
C.
,
2004
,
Intraplate tectonics: Feedback between radioactive thermal weakening and crustal deformation driven by mantle lithosphere instabilities
:
Earth and Planetary Science Letters
 , v.
221
, no.
1–4
, p.
275
292
. https://doi.org/10.1016/S0012-821X(04)00098-6.
77.
Quinteros
,
J.
, and
Sobolev
,
S.V.
,
2013
,
Why has the Nazca plate slowed since the Neogene?
:
Geology
 , v.
41
, no.
1
, p.
31
34
. https://doi.org/10.1130/G33497.1.
78.
Ramos
,
V.A.
,
1988
,
Tectonics of the Late Proterozoic early Paleozoic: A collisional history of southern South America
:
Episodes
 , v.
11
, no.
3
, p.
168
174
. https://doi.org/10.18814/epiiugs/1988/v11i3/003.
79.
Ramos
,
V.A.
,
2004
,
Cuyania, an exotic block to Gondwana: Review of a historical success and the present problems
:
Gondwana Research
 , v.
7
, no.
4
, p.
1009
1026
. https://doi.org/10.1016/S1342-937X(05)71081-9.
80.
Ramos
,
V.A.
, and
Folguera
,
A.
,
2009
,
Andean flat-slab subduction through time
, in
Murphy
,
J.B.
, et al
, eds.,
Ancient Orogens and Modern Analogues: Geological Society
 ,
London
,
Special Publication 327
, no.
1
, p.
31
54
. https://doi.org/10.1144/sp327.3.
81.
Ramos
,
V.A.
,
Cristallini
,
E.O.
, and
Pérez
,
D.J.
,
2002
,
The Pampean flat-slab of the Central Andes
:
Journal of South American Earth Sciences
 , v.
15
, no.
1
, p.
59
78
. https://doi.org/10.1016/S0895-9811(02)00006-8.
82.
Ramos
,
V.A.
,
Vujovich
,
G.
,
Martino
,
R.
, and
Otamendi
,
J.
,
2010
,
Pampia: A large cratonic block missing in the Rodinia supercontinent
:
Journal of Geodynamics
 , v.
50
, no.
3–4
, p.
243
255
. https://doi.org/10.1016/j.jog.2010.01.019.
83.
Ranero
,
C.R.
,
Phipps Morgan
,
J.
,
McIntosh
,
K.
, and
Reichert
,
C.
,
2003
,
Bending-related faulting and mantle serpentinization at the Middle America trench
:
Nature
 , v.
425
, p.
367
373
. https://doi.org/10.1038/nature01961.
84.
Rosenbaum
,
G.
,
Giles
,
D.
,
Saxon
,
M.
,
Betts
,
P.G.
,
Weinberg
,
R.F.
, and
Duboz
,
C.
,
2005
,
Subduction of the Nazca Ridge and the Inca Plateau: Insights into the formation of ore deposits in Peru
:
Earth and Planetary Science Letters
 , v.
239
, no.
1–2
, p.
18
32
. https://doi.org/10.1016/j.epsl.2005.08.003.
85.
Ruiz
,
F.
, and
Introcaso
,
A.
,
2004
,
Curie point depths beneath Precordillera Cuyana and Sierras Pampeanas obtained from spectral analysis of magnetic anomalies
:
Gondwana Research
 , v.
7
, no.
4
, p.
1133
1142
. https://doi.org/10.1016/S1342-937X(05)71089-3.
86.
Sánchez
,
M.A.
,
Ariza
,
J.P.
,
García
,
H.P.
,
Gianni
,
G.M.
,
Weidmann
,
M.C.
,
Folguera
,
A.
,
Lince Klinger
,
F.
, and
Martinez
,
M.P.
,
2018
,
Thermo-mechanical analysis of the Andean lithosphere over the Chilean-Pampean flat-slab region
:
Journal of South American Earth Sciences
 , v.
87
, p.
247
257
. https://doi.org/10.1016/j.jsames.2017.09.036.
87.
Schmidt
,
M.W.
, and
Poli
,
S.
,
1998
,
Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation
:
Earth and Planetary Science Letters
 , v.
163
, no.
1–4
, p.
361
379
. https://doi.org/10.1016/S0012-821X(98)00142-3.
88.
Scire
,
A.
,
Zandt
,
G.
,
Beck
,
S.
,
Long
,
M.
,
Wagner
,
L.
,
Minaya
,
E.
, and
Tavera
,
H.
,
2016
,
Imaging the transition from flat to normal subduction: Variations in the structure of the Nazca slab and upper mantle under southern Peru and northwestern Bolivia
:
Geophysical Journal International
 , v.
204
, no.
1
, p.
457
479
. https://doi.org/10.1093/gji/ggv452.
89.
Smith
,
D.
,
2020
,
Trace elements in Cr-pyrope from the Navajo volcanic field of the Colorado Plateau, SW USA, and implications for the mantle wedge during low-angle subduction
:
Lithos
 , v.
362–363
, no.
105460
, https://doi.org/10.1016/j.lithos.2020.105460.
90.
Stein
,
C.A.
, and
Stein
,
S.
,
1992
,
A model for the global variation in oceanic depth and heat flow with lithospheric age
:
Nature
 , v.
359
, no.
6391
, p.
123
129
. https://doi.org/10.1038/359123a0.
91.
Stevens Goddard
,
A.L.
, and
Carrapa
,
B.
,
2018
,
Using basin thermal history to evaluate the role of Miocene–Pliocene flat‐slab subduction in the southern Central Andes (27° S–30° S)
:
Basin Research
 , v.
30
, no.
3
, p.
564
585
. https://doi.org/10.1111/bre.12265.
92.
Syracuse
,
E.M.
, and
Abers
,
G.A.
,
2006
,
Global compilation of variations in slab depth beneath arc volcanoes and implications
:
Geochemistry, Geophysics, Geosystems
 , v.
7
, no.
5
, https://doi.org/10.1029/2005GC001045.
93.
Syracuse
,
E.M.
,
van Keken
,
P.E.
, and
Abers
,
G.A.
,
2010
,
The global range of subduction zone thermal models
:
Physics of the Earth and Planetary Interiors
 , v.
183
, no.
1–2
, p.
73
90
. https://doi.org/10.1016/j.pepi.2010.02.004.
94.
Thomas
,
W.A.
, and
Astini
,
R.A.
,
1996
,
The Argentine precordillera: A traveler from the Ouachita embayment of North American Laurentia
:
Science
 , v.
273
, no.
5276
, p.
752
757
. https://doi.org/10.1126/science.273.5276.752.
95.
Turcotte
,
D.L.
, and
Schubert
,
G.
,
2002
,
Geodynamics
:
New York
,
Cambridge University Press
,
456
p. https://doi.org/10.1017/CBO9780511807442.
96.
Urbina
,
N.P.
,
Sruoga
,
P.
, and
Malvicini
,
L.
,
1997
,
Late Tertiary gold-bearing volcanic belt in the Sierras Pampeanas of San Luis, Argentina
:
International Geology Review
 , v.
39
, no.
4
, p.
287
306
. https://doi.org/10.1080/00206819709465272.
97.
Usui
,
T.
,
Nakamura
,
E.
,
Kobayashi
,
K.
,
Maruyama
,
S.
, and
Helmstaedt
,
H.
,
2003
,
Fate of the subducted Farallon plate inferred from eclogite xenoliths in the Colorado Plateau
:
Geology
 , v.
31
, no.
7
, p.
589
592
. https://doi.org/10.1130/0091-7613(2003)031<0589:FOTSFP>2.0.CO;2.
98.
Uyeda
,
S.
, and
Watanabe
,
T.
,
1982
,
Terrestrial heat flow in western South America: Tectonophysics
, v.
83
, no.
1–2
, p.
63
70
. https://doi.org/10.1016/0040-1951(82)90007-5.
99.
van Hunen
,
J.
,
van den Berg
,
A.P.
, and
Vlaar
,
N.J.
,
2004
,
Various mechanisms to induce present-day shallow flat subduction and implications for the younger Earth: A numerical parameter study
:
Physics of the Earth and Planetary Interiors
 , v.
146
, no.
1–2
, p.
179
194
. https://doi.org/10.1016/j.pepi.2003.07.027.
100.
von Huene
,
R.
,
Corvalán
,
J.
,
Flueh
,
E.R.
,
Hinz
,
K.
,
Korstgard
,
J.
,
Ranero
,
C.R.
, and
Weinrebe
,
W.
,
1997
,
Tectonic control of the subducting Juan Fernández Ridge on the Andean margin near Valparaiso, Chile
:
Tectonics
 , v.
16
, no.
3
, p.
474
488
. https://doi.org/10.1029/96TC03703.
101.
Wada
,
I.
, and
Wang
,
K.
,
2009
,
Common depth of slab‐mantle decoupling: Reconciling diversity and uniformity of subduction zones
:
Geochemistry, Geophysics, Geosystems
 , v.
10
, no.
10
, https://doi.org/10.1029/2009GC002570.
102.
Wagner
,
L.S.
,
Beck
,
S.
, and
Zandt
,
G.
,
2005
,
Upper mantle structure in the south central Chilean subduction zone (30° to 36°S)
:
Journal of Geophysical Research: Solid Earth
 , v.
110
, https://doi.org/10.1029/2004jb003238.
103.
Wagner
,
L.S.
,
Beck
,
S.
,
Zandt
,
G.
, and
Ducea
,
M.N.
,
2006
,
Depleted lithosphere, cold, trapped asthenosphere, and frozen melt puddles above the flat slab in central Chile and Argentina
:
Earth and Planetary Science Letters
 , v.
245
, no.
1–2
, p.
289
301
. https://doi.org/10.1016/j.epsl.2006.02.014.
104.
Wagner
,
L.S.
,
Anderson
,
M.L.
,
Jackson
,
J.M.
,
Beck
,
S.L.
, and
Zandt
,
G.
,
2008
,
Seismic evidence for orthopyroxene enrichment in the continental lithosphere
:
Geology
 , v.
36
, no.
12
, p.
935
938
. https://doi.org/10.1130/G25108A.1.
105.
Wagner
,
L.S.
,
Caddick
,
M.J.
,
Kumar
,
A.
,
Beck
,
S.L.
, and
Long
,
M.D.
,
2020
,
Effects of oceanic crustal thickness on intermediate depth seismicity
:
Frontiers in Earth Sciences
 , v.
8
, no.
244
, https://doi.org/10.3389/feart.2020.00244.
106.
Warren
,
C.J.
,
Beaumont
,
C.
, and
Jamieson
,
R.A.
,
2008
,
Formation and exhumation of ultra‐high-pressure rocks during continental collision: Role of detachment in the subduction channel
:
Geochemistry, Geophysics, Geosystems
 , v.
9
,
Q04019
, https://doi.org/10.1029/2007GC001839.
107.
Woods
,
M.T.
, and
Okal
,
E.A.
,
1994
,
The structure of the Nazca ridge and Sala y Gomez seamount chain from the dispersion of Rayleigh waves
:
Geophysical Journal International
 , v.
117
, no.
1
, p.
205
222
. https://doi.org/10.1111/j.1365-246X.1994.tb03313.x.
108.
Yáñez
,
G.
,
Cembrano
,
J.
,
Pardo
,
M.
,
Ranero
,
C.
, and
Selles
,
D.
,
2002
,
The Challenger–Juan Fernández–Maipo major tectonic transition of the Nazca–Andean subduction system at 33–34°S: Geodynamic evidence and implications
:
Journal of South American Earth Sciences
 , v.
15
, no.
1
, p.
23
38
. https://doi.org/10.1016/S0895-9811(02)00004-4.
109.
Yáñez
,
G.A.
,
Ranero
,
C.R.
,
von Huene
,
R.
, and
Díaz
,
J.
,
2001
,
Magnetic anomaly interpretation across the southern central Andes (32–34°S): The role of the Juan Fernández Ridge in the late Tertiary evolution of the margin
:
Journal of Geophysical Research: Solid Earth
 , v.
106
, p.
6325
6345
. https://doi.org/10.1029/2000JB900337.
110.
Ziagos
,
J.P.
,
Blackwell
,
D.D.
, and
Mooser
,
F.
,
1985
,
Heat flow in southern Mexico and the thermal effects of subduction
:
Journal of Geophysical Research: Solid Earth
 , v.
90
, p.
5410
5420
. https://doi.org/10.1029/JB090iB07p05410.
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.