Reconciling Geophysical and Petrological Estimates of the Thermal Structure of Southern Tibet

The thermal structure of the Tibetan plateau—the largest orogenic system on Earth—remains largely unknown. Numerous avenues provide fragmentary pressure/temperature information, both at the present (predominantly informed though geophysical observation) and on the evolution of the thermal structure over the recent past (combining petrological, geochemical, and geophysical observables). However, these individual constraints have proven hard to reconcile with each other. Here, we show that models for the simple underthrusting of India beneath southern Tibet are capable of matching all available constraints on its thermal structure, both at the present day and since the Miocene. Many parameters in such models remain poorly constrained, and we explore the various trade‐offs among the competing influences these parameters may have. However, three consistent features to such models emerge: (i) that present‐day geophysical observations require the presence of relatively cold underthrust Indian lithosphere beneath southern Tibet; (ii) that geochemical constraints require the removal of Indian mantle from beneath southern Tibet at some point during the early Miocene, although the mechanism of this removal, and whether it includes the removal of any crustal material, is not constrained by our models; and (iii) that the combination of the southern extent of Miocene mantle‐derived magmatism and the present‐day geophysical structure and earthquake distribution of southern Tibet require that the time‐averaged rate of underthrusting of India relative to central Tibet since the middle Miocene has been faster than it is at present.


Introduction
As the largest modern continental orogeny, the Tibetan plateau represents an ideal modern analog for major collisional orogenic systems since the onset of modern-style tectonics. The temperature structure and evolution of the Tibetan plateau is of vital importance in understanding not only the seismicity, volcanism, rheology, and deformation of the modern collision zone but also the evolution of orogenic systems in general. Significant areas of the continental interiors are dominated by the geological products of such collision zones, with their chemical and rheological history governed by their thermal evolution, and the resultant metamorphic and igneous processes they were subjected to, during orogenic growth and decay. Despite this importance, the thermal structure of the lower crust and upper mantle of Tibet remains a controversial topic.
Thermochronology of rocks exposed along the Himalayan front offers a window into the kinematics and shallow thermal structure of the orogenic front and Himalayan prism (e.g., Bollinger et al., 2006;Godin et al., 2001;Herman et al., 2010). However, the incorporation of material into the Himalayas only represents a relatively small volume of the total lithosphere involved in the collision zone. Constraining the thermal structure of the deeper Tibetan plateau is more difficult, and despite extensive work, there are relatively few available constraints on either its contemporary structure or its evolution through time. Many of the constraints that do exist are hard to reconcile into a consistent thermal structure, particularly at depth beneath the plateau. Similarly, and in part due to the paucity of observational constraints on the deeper rheological structure of Tibet, debate continues about the geodynamic evolution of the plateau. Processes including large-scale underthrusting, crustal channel flow, gravitational spreading, convective delamination and slab tearing, break-off, and removal are all invoked to explain particular geological, geochemical, and geophysical observations.
Here, we focus on the deep thermal evolution of the southern half of the Tibetan plateau, limited to within ∼500 km of the Himalayan Front, hereafter referred to as "south/southern Tibet." We summarize existing proxies for the thermal and rheological structure of southern Tibet and present a series of two-dimensional models for the thermal evolution of the region from the Miocene to present, to investigate what models are compatible with available geophysical and geochemical constraints. As we shall discuss, these models do not provide a unique solution for the structure and evolution of southern Tibet but instead demonstrate that reconciling available geophysical and geochemical constraints is possible within a coherent framework and highlight some consistent features that are common to all successful models.

Geophysical Constraints
Receiver functions, when combined with co-located surface-wave dispersion studies, provide vital information on the present-day crustal geometry of the collision zone (Figures 1 and 2). In general, the crustal thickness varies from 35-40 km across much of peninsular India to 75-80 km under southern Tibet, with this increase being roughly linear over the first 200-250 km of the collision zone, north of the Himalayan front, where incoming Indian material is underthrust beneath Tibet. While there is significant local variability in this trend (e.g., Subedi et al., 2018), particularly in the region of most rapid change (and most variable geometry) beneath the high Himalaya (e.g., Hubbard et al., 2016;Ingleby et al., 2020), the available information supports a broadly two-dimensional crustal-scale pattern in southern Tibet, perpendicular to the southern edge of the Plateau.
This simple picture, of Indian lithosphere thrusting northwards beneath the Tibetan plateau in a manner that, in terms of the overall crustal thickness, varies little along strike, is also consistent with tomographic imaging of the deeper seismological structure of the plateau. Surface-wave tomography at frequencies appropriate to both crustal (Gilligan & Priestley, 2018) and mantle (M. Chen et al., 2017) levels shows a gradually thickening crustal pile beneath the Himalayas with low seismic wave speeds, an approximately uniform crustal thickness beneath the Tibetan plateau itself, and a seismically fast uppermost mantle, consistent with the presence of cold India-derived lithosphere beneath the southern plateau. Due to its different sensitivity, body-wave tomography (e.g., Li & Song, 2018;Liang et al., 2016;Ren & Shen, 2008) shows wider lateral variability in this seismological picture. Such studies commonly detect intraplateau features, such as the graben systems of southern Tibet, but match the broad image of a seismologically fast region of Indian lithosphere underlying the thickened crust of southern Tibet.
The underthrusting of Indian crustal material, along with progressive mineralogical and rheological evolution northwards, is also suggested by modeling the variation in terrestrial gravity observations across the Himalayas and onto southern Tibet (Cattin et al., 2001;Hétenyi et al., 2007). While such models only provide a relatively coarse image of the density distribution at depth beneath the plateau, and their petrological implications are subject to the poorly known thermal structure, they indicate a progressive densification of Indian crustal material over the first ∼400 km of the plateau, most probably driven by fluid release during dehydration of underthrust Indian material as it heats up.
We now consider geophysical constraints on the temperature distribution of India and southern Tibet. With the well-established dependence of seismic slip on the temperature and composition of lithospheric material, the distribution of earthquakes within the collision zone provides an initial constraint on the thermal structure of the plateau (Figures 1a and 2b). Across much of the Tibetan plateau itself, the depth extent of small-and moderate-magnitude earthquakes (Craig et al., 2012;Langin et al., 2003;Sloan et al., 2011), the rupture extent of large earthquakes (Elliott et al., 2010;Wright et al., 2013), and the geodetically determined locking depth of major crustal faults (Wright et al., 2013) all terminate at ∼12to 13-km depths, with this transition in crustal rocks typically occurring at ∼350°C (W.-P. Jackson et al., 2008).
In peninsular India, seismicity is observed to occur throughout the thickness of the crust in all areas where major seismic activity has been observed (northwestern peninsular India, around the Shillong plateau, and beneath the Himalayan foreland) and in some cases into the top of the uppermost lithospheric mantle (Craig et al., 2012;Monsalve et al., 2006;Schulte-Pelkum et al., 2019). Where Indian material is underthrust beneath the Tibetan plateau, a narrowing tongue of seismicity extends beneath the Himalayas and southern Tibet following the Indian lower crust and Moho (Monsalve et al., 2006;Schulte-Pelkum et al., 2019), gradually dying out at 400-500 km north of the Himalayan front (Craig et al., 2012;Priestley et al., 2008). These  Craig et al. (2012), updated with Sippl et al. (2013), Mitra et al. (2014), Paul et al. (2015), Reynolds et al. (2015), Kufner et al. (2016), Mendoza et al. (2016), Ainscoe et al. (2017), Diehl et al. (2017), Gibbons and Kvaerna (2017), Huang et al. (2017), Negi et al. (2017), Kanna et al. (2018), and Parija et al. (2018). Dashed lines highlight areas of deep seismicity associated with the Hindu-Kush deep seismic zone and the Indo-Burman subduction zone and also the northern boundary of the Tibetan plateau. (b) Receiver-function observations of crustal thickness (after Gilligan & Priestley, 2018). (c) Shear wave splitting measurements (after Agius & Lebedev, 2017). (d) Volcanism across Tibet, where both chemistry and age are well constrained (after Yakovlev et al., 2019). (e) Locations of xenolith suites (Chan et al., 2009;Ding et al., 2007;Hacker et al., 2000;Liu et al., 2011), seismological observations of the inferred ABQT (Mechie et al., 2004;Sheehan et al., 2014), and the region of low S n propagation in northern Tibet (McNamara & Owens, 1995). Dashed white line indicates ∼500 km from the Himalayan front along the India-Asia convergence vector. earthquakes represent the internal deformation of underthrust Indian lithosphere. The rate of earthquake occurrence is low and is concentrated into two regional groups in northwestern and eastern South Tibet. This patchy distribution may either reflect the slow strain rates in the host material and hence may not be a representative seismic catalog (similar to that seen in peninsular India) or a further control on seismicity based on spatial variations in the mineralogical structure of underthrust India (Alvizuri & Hetényi, 2019;Jamtveit et al., 2019). In general, seismicity associated with either the upper mantle (Boettcher et al., 2007;Craig et al., 2014) or a mafic and anhydrous lower crust (Jackson et al., 2008;Mackwell, 1998) can persist to temperatures of at least 600-650°C, before, at standard tectonic strain rates, deformation becomes ductile. Beneath the Himalaya, the pattern in seismicity is more complex, as shown on Figure 2, accommodating this bifurcation into shallow upper crustal seismicity, and deeper seismicity roughly along the Moho Satellite magnetic measurements (Alsdorf & Nelson, 1999), although lacking in spatial resolution, show a magnetic low across the plateau associated with a shallow depth (∼15 km) of the crustal Curie temperature (∼550°C at 15 km, from which they infer that temperatures of 600-650°C are reached regionally at 16-to 18-km depth; Figure 2d).
Limited surface heat flow measurements from across Tibet offer a direct estimate of the shallow geotherm (International Heat Flow Commission, 2011). However, heat flow measurements are logistically complex and expensive, are often taken in nonrepresentative sections of the plateau (rift basins, geothermal Xenolith data for Namring (Chan et al., 2009). Dashed boxes show P/T uncertainties. (f) Xenolith data for Sailipu (Liu et al., 2011). Inset histogram shows the temperature population of the xenolith suite. (g) Xenolith data for Yibuchaka (Ding et al., 2007). Inset histogram shows the temperature population of the xenoliths suite, with lithology indicated by color (gray for mafic granulite, green for olivine-bearing, and black for websterite). (h) Xenolith data for Taipinghu (Hacker et al., 2000). Red line is the temperature profile used in thermal models for the northern boundary condition. Inset dates on panels (e)-(h) show the age of emplacement, with the age of peak metamorphism, if known, given in brackets. See Figure 1 for the locations of xenolith suites shown in panels (e)-(h). areas, etc.), and may be strongly influenced by shallow radiogenic structure or recent magmatism. These factors result in a wide distribution of values, spanning a range of 200 mW m −2 . As a result of this uncertainty over which values may be representative of the deeper thermal structure of the plateau, we exclude surface heat flow measurements from within the plateau from the rest of our analysis.
Seismological data provide a range of constraints on the present-day temperature structure of the plateau (Figure 2d). Densely spaced seismic acquisitions in the northern Himalayas (HIMNT; Sheehan et al., 2014) and southern Tibet (InDEPTH III; Mechie et al., 2004) detected high reflectivity contrasts in V p /V s ratios linked to high seismic attenuation, inferred to represent the transition from α to β quartz within the crust (hereafter referred to as the ABQT). While questions remain about why this transition has not been detected more commonly across the plateau, laboratory constraints on the pressure-temperature controls of this transition (Angel et al., 2017) allow an estimate of the temperature at the depths where this interface is inferred, at depths from ∼45-20 km, shallowing northwards (Figures 1e and 2d).
The variation in seismic anisotropy revealed by both shear-wave tomography and shear-wave splitting across the plateau provides critical information on the geometry of strain. Recent shear-wave tomography (Gilligan & Priestley, 2018) shows the development of significant radial anisotropy in the midcrust of the central Tibetan plateau, in contrast to the seismic structure of both India and southernmost Tibet, where anisotropy is low. Similarly, shear-wave splitting (Figures 1c and 2c; after Agius & Lebedev, 2017), while showing wide variation across the collision zone (which may partly reflect the range of data selection and analysis methods used), shows a consistent northwards increase in amplitude and, more importantly, a shift toward a consistent preferential alignment in the fast direction, such that more than ∼300 km north of the Himalayan front, the fast direction is consistently 070-090°E. The development of a coherent anisotropic structure within the central plateau is suggestive of a process in which gradual heating of underthrust Indian lithosphere, as it moves northwards beneath Tibet, leads to its progressive rheological weakening and softening. This weakening allows increasing volumes of crustal and mantle material (starting with the midcrust) to deform and flow, resulting in a situation where a weaker, ductile material from central and northern Tibet flows out of the collision zone to the east, leading to the development of a preferred mineral orientation, resulting in both a coherent shear-wave splitting fabric and midcrustal anisotropy (e.g., Gilligan & Priestley, 2018;Pandey et al., 2015).
Beneath northern Tibet, frequency-dependent inefficient propagation of the seismic phases P n and S n (head waves that propagate along the Moho) indicates an inverted gradient in seismic velocities beneath the Moho, over a depth extent that increases northwards ( Figure 2d; Barron & Priestley, 2009;McNamara & Owens, 1995;Ni & Barazangi, 1983). While this observation does not provide an absolute estimate of the temperature at this point, it does suggest an inversion of the vertical temperature gradient across the Moho (ie., the upper mantle just below the Moho may be hotter than the mantle beneath it).

Petrological and Geochemical Constraints
While geophysical observations are limited to direct constraints on the present-day structure of the plateau, geochemical information offers the opportunity to constrain the past temperature at depth within the plateau. In all cases, the present-day surface location of the samples providing these geochemical constraints is displaced relative to the deep structure of the plateau at the time they were emplaced, due to subsequent motion within the collision zone. Hence, sample locations must be corrected for the movement of material within the collision zone in the intervening time since their age of chemical equilibration.
Southern Tibet has experienced numerous phases of volcanism over the last 50 Myr. Here, we draw on a recent compilation by Yakovlev et al. (2019) for volcanic rocks within the Tibetan plateau for which both detailed geochemical and chronological data are available (Figures 1d and 2d). Erupted volcanics with Mg# ≪ 0.7, particularly those in northern Tibet, likely reflect either intracrustal melting (for more recent samples), the effects of arc volcanism (in older samples), or highly fractionated mantle-derived melts. While intracrustal melting has the potential to add some further constraints on the recent temperature at depth within the plateau (particularly at midcrustal depths), this would require a level of detailed analysis (and knowledge of the source rock composition) that is beyond the scope of this study. Instead, we exclude samples with Mg# < 0.7 and focus on those samples with Mg# > 0.7 (in southern Tibet, these are typically potassic or ultrapotassic rocks), as a method for identifying samples which likely have equilibrated with refractory peridotite (Mg# > 0.89), and hence derive from the melting of a source at subcrustal levels. The generation of such melts requires temperatures sufficient to partially melt a source rock-either mantle peridotite or various metasomatic or subducted components-and to equilibrate that melt with mantle peridotite. Obviously, such processes must take place at depths where mantle peridotite is present, below the Moho. We then assume that the duration of intracrustal storage for such rocks is negligible compared to the timespan of the India-Asia collision (e.g., Mutch et al., 2019), and thus that eruption ages are representative of the age of melting. This argument is supported by the requirement for a high Mg#, indicating that the magma has undergone limited fractionation and/or crustal assimilation. The spatial distribution of such volcanics hence constrains the age and extent of widespread (but predominantly small melt-fraction) mantle melting beneath the plateau. The detailed data set, including a breakdown of Mg#, rather than the binary division used here, is shown in the supporting information Figure S1 and fully described in Yakovlev et al. (2019).
Much of the older volcanism across the plateau prior to 50 Ma is likely associated with precollisional arc-type volcanism in the accreted terranes of central Tibet. From 50 to ∼30 Ma, central and southern Tibet show little magmatic activity. After 30 Ma, there is a progressive spread of magmatism, including magmatism with Mg# > 0.7, from central Tibet both northwards and southwards. In northern Tibet, magmatic activity with Mg# < 0.7 continues to the present day. In southern and central Tibet, the majority of magmatic activity (including that with Mg# > 0.7) terminates at or before ∼12 Ma (and predominantly by 15 Ma). This general trend in magmatic activity is mirrored along the length of Himalayan arc, although there is some variability in the time at which various phases of this behavior occurred, suggesting that the underlying geodynamic processes were diachronous along strike (Webb et al., 2017). Crucially, the extent of Mg# > 0.7 magmatism in south and central Tibet into the mid-Miocene indicates that the mantle beneath this region at this time was hot enough, and shallow enough, to melt-a feature that ends by 12 Ma.
More accurate constraints on temperature at depth can be obtained from detailed analysis of xenoliths contained within some of these volcanics. Here, we draw on the four published sets of xenolith data for Tibet (excluding the Pamir): from Taipinghu in the northern Qiangtang terrane (Hacker et al., 2000); Yibuchaka, southern Qiangtang terrane (Ding et al., 2007); Namring, just north of the Indus-Yarlung suture (Chan et al., 2009); and Sailipu, western Lhasa terrane (Liu et al., 2011), all shown on Figure 1e.
In the case of Sailipu (Liu et al., 2011), the xenoliths are peridotite, entrained in an ultrapotassic host rock. The ultrapotassic lava most probably formed via melting of mantle material. For Yibuchaka (Ding et al., 2007), the xenoliths are a mixture of mafic granulites, lherzolites, and websterites and therefore span the Moho, from lower crust to uppermost mantle. While accurate temperature estimates exist for samples from these two locations (see Figures 2f and 2g), the pressure constraints for the Sailipu and Yibuchaka xenolith suites are based on the absence of garnet, based on which they are inferred to originate above the garnet-spinel transition and are therefore only an upper limit. For the xenoliths from Yibuchaka and Sailipu, only eruption ages are known (28 and 17 Ma, respectively), not the age at which the xenoliths were in chemical equilibrium with their surroundings. As a result, it is unclear whether these xenoliths represent the thermal structure of the plateau at, or shortly before, their time of emplacement, or instead reflect prior processes, (and potentially arc-related volcanism), that affected the plateau. We therefore do not consider them further.
The best-constrained set of xenolith data for southern Tibet comes from an andesitic/potassic dyke near Namring (Chan et al., 2009). These comprise micaceous, felsic, and mafic granulite, and ultramafic xenoliths. All have well-constrained temperature estimates, and three have well-constrained pressure estimates (shown on Figure 2e). Crucially, both the age of emplacement (12.7 Ma) and the age of peak metamorphism (16.8, 15.6, and 14.4 Ma for the three separate xenoliths shown in Figure 2e) are known. Hence, the pressure/temperature estimates from this xenolith suite can be used to constrain the thermal structure of the plateau in that time range.
Recently erupted crustal (metasedimentary and mafic) xenoliths from Taipinghu (Hacker et al., 2000) with robust thermobarometry provide information on the near-present-day geotherm of north-central Tibet ( Figure 2h). Additionally, the near-anhydrous nature of these xenoliths helps to explain why, despite the high inferred temperatures, the Tibetan midcrust does not host widespread melt.

Geochemistry, Geophysics, Geosystems
Significant work has been done by numerous authors in constraining the thermochronological and thermobarometric evolution of the Himalayas, in many cases as a way of investigating the structural evolution and sequence of faulting in the Himalayan prism. Much of this has focused on low temperature (e.g., 40 Ar/ 39 Ar; apatite fission track, U-Th-He) thermometers and do not constrain the deeper thermal structure of the plateau. Additionally, as we shall discuss, the finer-scale detail of the faulting regime in the Himalayas, while having a major influence on the shallow thermal structure of the mountain range, has little impact on the lower crustal thermal structure of southern Tibet.

The Thermal Structure of Southern Tibet
A difficulty with southern Tibet lies in reconciling the geophysical data, that imply a relatively intact India underthrust beneath the southern plateau, with geochemical data that require widespread mantle melting beneath areas of the southern plateau until ∼10 Ma (Chung et al., 2005;Yakovlev et al., 2019), with inferred temperatures around the Moho that are inconsistent with the underthrusting of cool Indian lithosphere.
One proposed solution has been the removal of large amounts of lithospheric mantle from beneath the plateau, most probably through the break-off of underthrust Indian lithosphere (DeCelles et al., 2002;Guillot et al., 2003;Replumaz et al., 2010). Over the last decade, increasing evidence, ranging from the identification of a relict "Indian" slab in the deeper mantle beneath Tibet (Replumaz et al., 2010(Replumaz et al., , 2014, the progressive migration and termination of magmatic activity across Tibet (Chung et al., 2005;Webb et al., 2017;Yakovlev et al., 2019), changes in sedimentation in both the Himalayan foreland (Mugnier & Huyghe, 2006), in Northern India (Najman et al., 2018), and within the plateau interior (Carrapa et al., 2014;Leary et al., 2016), and a switch from contraction to extension (Stearns et al., 2013), all support a model in which an initial phase of underthrusting by colder Indian material was followed by a prolonged phase of southward slab retreat back under southern Tibet, ending with at least one break-off, which was then followed by renewed underthrusting of Indian material beneath the plateau-a model with major implications for the present-day temperature field beneath the plateau in the lower crust and upper mantle.
These observations therefore raise the question of whether a phase of slab break-off, and the implied removal of Indian material from beneath the plateau, is compatible with the presence of hot upper mantle recently enough, and shallowly enough, to match the geochemical data described above, while also being able, following re-underthrusting of continental lithosphere from India, to match geophysical constraints on the present-day plateau structure. To address this question, we model the thermal evolution of the southern Tibetan plateau, focusing on the last 20 Myr, to investigate the range of options for the kinematic evolution of the plateau that fit all available geophysical and geochemical constraints on the geological and thermal structure of the southern Tibetan plateau. Here, we specifically focus on the thermal structure at depths greater than ∼20 km-an approach that allows us to simplify the shallow fault-dominated kinematics of the collision zone as expressed in the complexity of Himalayan tectonics.

Thermal Modeling
We construct a set of two-dimensional forward models for the thermal evolution of the southern Tibetan plateau, solving for the advection, diffusion, and internal production of heat through time using where T is temperature, t is time, v is the advective velocity field, H r is radiogenic heat production, and ρ, C p , and k are the density, heat capacity, and thermal conductivity, respectively. We use a finite difference approach and separately solve for the diffusive and advective terms at each time step using operator splitting. An iterative Crank-Nicolson approach is used for the diffusive step and an upwind scheme for the advective step.
In the following sections, we show that it is possible to fit all available constraints with a range of combinations of model parameters. We discuss the uniqueness of the parameter choices we make, the various trade-offs among parameters, and what ranges are plausible for those parameters.

Geochemistry, Geophysics, Geosystems
Our models are kinematic, rather than dynamic, and the model geometry and the advective velocity field are prescribed. We follow previous approaches (e.g., Bollinger et al., 2006;Cattin et al., 2001;Craig et al., 2012;Henry et al., 1997;Herman et al., 2010;P. Nábělek & Nábělek, 2014;Wang et al., 2013) in modeling a scenario in which Indian material is underthrust beneath Asia-derived crustal material, with the accretion of Indian upper crustal material into the Himalayan prism, and with the erosive removal of material in the Himalayas. Model kinematics and boundary conditions are summarized in Figure 3a.
For completeness, we initially run our model, with cool Indian material being underthrust beneath Tibet, for (50−t b ) Myr, where t b is the time of material removal. As an approximation to the removal of Indian material through slab break-off (or by any other mechanism), at time t b , we replace all Indian-derived material beneath a given depth, and north of x b km beyond the Himalayan front, with mantle material at the temperature of the mantle adiabat (Figure 3b). We then continue underthrusting Indian material for a further t b Myr, corresponding to the time between "break-off" and the present day, during which the emplaced material is allowed to cool conductively, until we reach the model "present day" (Figure 3c).
Major uncertainties in such models arise from a variety of underlying causes: The thermal profiles used as boundary conditions for incoming "India lithosphere" and overthrusting "Tibetan crust" the distribution of radiogenic heat production, particularly in the Tibetan overthrust; the model kinematics, their partitioning between overthrusting, underplating through the Himalayan megathrust, and the removal of material through the Himalayan front, and their evolution through time; and the timing, location, and mode of removal from the lower sections of the plateau interior. We address each of these in turn in the following sections, describing the range of plausible values, and the effects that they have on model outputs.
For the thermal structure of incoming Indian lithosphere, we use a model in one-dimensional steady state, based on cratonic areas of Indian south of the Himalayan front ( Figure 4). We assume a lithospheric thickness of 200 km, consistent with geophysical estimates based on the variation of shear wave speed with depth, and with pressure-temperature estimates from a suite of mantle xenoliths from Wajrakarur, beneath the Dharwar craton (Roy & Mareschal, 2011). In all cases, we use 0.02 μW m −3 for the lithospheric mantle heat production, although this value has a negligible effect. Constraints on crustal radiogenic heat production are sparse, largely informed by surface exposure of geological units believed to be representative of the deeper crustal column. In consequence, we therefore trial a range of models using crustal heat production columns based on surface measurements of radiogenic heat production, combined with surface heat flow measurements.
In our preferred model of the initial steady-state Indian thermal structure (shown in Figure 4), based on crustal data from the cratonic western Dharwar province (Roy & Mareschal, 2011), we use a 19-km upper crust of tonalite-trondhjemite-granodiorite gneiss, overlying a 19-km lower crust of mafic granulite, with radiogenic heat production of 1.1 and 0.36 μW m −3 , respectively, yielding a surface heat flux of 41.0 mW m −2 . Given the paucity of comprehensive data on crustal structure and heat production, we also trial three further models (Figure 4) for crustal heat production, based on values for the eastern Dharwar (Roy & Mareschal, 2011) and Bundelkhand (Podugu et al., 2017) provinces, all summarized in Table S1. Figure 5 then shows a comparison between our preferred model (western Dharwar) and the other three models. The model for Bundelkhand (Figures 4 and 5d) is significantly colder than any for Dharwar in terms of upper mantle temperatures (Moho temperature of ∼460°C, reaching 600°C at 56-km depth) and is likely too cold to represent the averaged thermal structure of India, as it would result in an extremely deep seismogenic thickness from, and an incredibly strong lithosphere-neither of which is supported by continent-wide data. Models for eastern Dharwar, with significantly higher radiogenic heat production in the upper Indian crust (2.4 or 2.0, rather than 1.1 μW m −3 ), lead to a hotter Himalayan prism (reaching a difference in places of 200°C), but only minor differences beneath southern Tibet, in either the relict Indian lower crust or Tibet-derived material. The difference between our models for western and eastern Dharwar ( Figure 5) shows that the thermal structure of the Himalayas is highly dependent on the distribution of radiogenic heating in the upper crust and as such is likely to be variable along strike and through time for the orogen, but that the deeper thermal structure of the lower crust and uppermost mantle is relatively insensitive to variations in shallow upper crustal radiogenic heating. Although beyond the scope of this study, we note the potential for shallow thermochronology from the Himalayas, combined with more advanced kinematics within the Himalayan prism, to further constraint the shallow thermal structure and 10.1029/2019GC008837 Geochemistry, Geophysics, Geosystems determine the appropriateness of these these two input models (e.g., Herman et al., 2010;Robert et al., 2011). From this point on, we show results only using the western Dharwar model.
In solving Equation 1, heat capacity, density, and thermal conductivity are all treated as temperature-dependent variables. We use experimentally derived parameterizations (see Figure S2) for For the initial temperature structure of the Tibetan overthrust, we use a linear geothermal gradient based on fitting crustal xenolith data from Taipinghu in northern Tibet (Figure 2h, Hacker et al., 2000), although we note that the location of these xenoliths within northern Tibet places them beyond the northward extent of our model. These xenoliths constrain the temperature beneath the northern plateau interior between 30-and 50-km depths. A linear geotherm fit through these data also matches the estimated regional Curie depths from satellite magnetic data across the whole plateau, putting the 550°C isotherm at ∼15 km. Similarly, the seismogenic thickness of the plateau interior, along with geodetic constraints on the locking depth of major intraplateau faults, suggests the seismic-aseismic transition at 10-12 km, also linked to the 350°C isotherm. All of these constraints are based on recent (<5 Ma) proxies and may not apply equally to the plateau back through time. We also note that a simple linear geotherm is unlikely to correctly represent the true thermal structure of north Tibetan crust. However, as will be discussed below, our more-complex thermal models produce a crustal thermal profile in central/northern Tibet equally able to fit the xenoliths data from Taipinghu, along with shallow geophysical constraints (see Figure 3f). Other than their use in constraining a crustal geotherm for northern Tibet and verifying that this does not vary significantly in the northern sections of our model, we do not use the Taipinghu xenoliths otherwise.
Few constraints exist on the distribution of radiogenic heating in the non-Indian-derived component of Tibetan crust. Much of this material was sourced from the various Cenozoic collision fronts and arc systems that accreted onto the southern edge of Asia during the closure of the Tethys, prior to contact with continental India at ∼55-50 Ma (e.g., Aitchison et al., 2007;Dewey et al., 1988;Zhu et al., 2013). As such, we can reasonably expect that this material will have been progressively enriched in incompatible radiogenic elements (including U and Th) as a result of arc magmatism. It is therefore likely that the average radiogenic heat production of Tibetan crustal material is higher than that of depleted continental lower crust, but the degree to which this enrichment has taken place and its distribution within the plateau crust are unconstrained. We hence adopt a model in which radiogenic heat production is uniform throughout the Tibetan crust and trial a series of models where the average heat production span a range of typical values of continental upper crust from 1.5 to 2.5 μW m −3 (Jaupart & Mareschal, 2003;Rudnick et al., 1998; results shown in Figure S3). The model shown in Figure 3 uses a value of 2.0 μW m −3 . As expected, increased internal heat production within the Tibetan overthrust leads to higher temperatures in the upper crust and midcrust of the plateau but has little effect on the deeper thermal structure within underthrust India ( Figure S3).
Convergence rates in our model are based on a discretized version of the India-Asia convergence rates of (Copley et al., 2011) (see Figure 7a). These rates, however, represent the total convergence between peninsula India and stable Eurasia, only part of which is accommodated by the Himalayas. Constraints from modern geodesy (Ader et al., 2012), fluvial morphology (Lavé & Avouac, 2000), and low-temperature thermochronology (Herman et al., 2010) in the Himalayas show that convergence between northern India and central Tibet presently accommodates roughly half the total convergence rate, although this partitioning in unlikely to have been constant through time-at the present, significant shortening is accommodated in the distributed deformation belts north of Tibet (the Tien Shan, Altai, Qilian Shan, etc.), but many of these did not become active until the later Miocene/Pliocene (e.g., Bullen et al., 2001;Lease et al., 2011;Molnar & Stock, 2009). Additionally, the removal of Indian material from beneath the plateau during the early Figure 4. Input temperature profiles for Indian-derived lithosphere using four crustal heat production models. Brown line indicates the Moho depth. Gray line shows the mantle adiabat, for a mantle potential temperature of 1,315°C. Black dots show P-T estimates for xenoliths from Wajrakarur (Roy & Mareschal, 2011). Colored lines show steady-state geotherms for four separate crustal heat production models. Numbers in brackets indicate surface heat flux.

Geochemistry, Geophysics, Geosystems
Miocene likely resulted in an increase in the rates of underthrusting, as forces acting to inhibit convergence decreased, due to higher temperatures and lower viscosities beneath the plateau. In the model shown in Figure 3, we include a 50% increase in the proportion of convergence accommodated along the Himalayas at the time of break-off (to 75% of the total convergence), tapered linearly to the present-day convergence rate (50% of the total). Noting that this partitioning is unconstrained, we test a variety of temporal variations in kinematics, as summarized in Figures 7 and S4 and discussed below.
The convergence rate between northern India and central Tibet is subdivided in our modeling approach into a rate of underthrusting (i.e., the motion of Indian material relative to the Himalayan front), and a rate of overthrusting (the motion of Tibetan material relative to the Himalayan front), with overthrusting countered by the removal of material via erosion across the Himalayas. In all our models, we use the present-day partitioning between these two components of ∼80% underthrusting to ∼20% overthrusting (Bollinger et al., 2004;Herman et al., 2010).
We include the accretion of upper-crustal material from the top of the incoming Indian plate into the Himalayan prism. The amount of material accreted is constrained by thermochronological data from the Himalayas (Bollinger et al., 2004;Herman et al., 2010) and from seismic images that detect an intracrustal layer in southern Tibet inferred to be the contact between underthrust Indian lower crust and Asian-derived material (J. Nábělek et al., 2009;Sheehan et al., 2014). Approximately half (15-20 km) of Indian crustal material is thought to be transferred from the top of the Indian plate into the Himalayan prism. While detailed modeling of the shallow thermal structure of the Himalayas indicates that accretion potentially takes place in a concentrated zone beneath the prism (Herman et al., 2010), leading to complex intraprism kinematics, such complexity has little effect on the deeper thermal structure, and so we instead model this accretion as taking place uniformly over a 250-km-wide zone where the Indian plate dips beneath southern Tibet. The model shown in Figure 3 uses 15 km of accretion from the top of the Indian plate into the Himalayas. As shown in Figure 6, varying the amount between 10 and 20 km has a major effect on the thermal structure of the Himalayan prism, and the temperature structure of the Tibetan midcrust, with increased accretion leading to a significantly hotter midcrust and Himalaya. Varying the amount of Figure 5. Final-stage thermal models using the four different input crustal heat production models shown in Figure 4. All other parameters are as in Figure 3 and discussed in the text. Left column shows thermal structure, and right column shows temperature differences relative to the model shown in Figure 3c. The top row shows the model from Figure 3. The second row shows a model using the crustal heat production values from eastern Dharwar from Roy and Mareschal (2011). The third row shows a model using a modified set of values for eastern Dharwar with a less-depleted lower crust. The fourth row shows a model using the heat production values from the Bundelkhand craton in northern India (Podugu et al., 2017).

10.1029/2019GC008837
Geochemistry, Geophysics, Geosystems accretion within the range of possible values leads to the second largest variation in the absolute values of the temperature structure beneath southern Tibet out of any of our model parameters, after the assumed structure for incoming Indian lithosphere. The effect on the Tibetan midcrust comes partly as a result of increased thickness of the Tibetan overthrust (required to maintain the same total crustal thickness beneath the plateau), and partly from the increased removal of the coldest sections of the incoming Indian crust. Overall, however, while absolute temperature values change, there is little effect on the general shape of the thermal structure.
For the purposes of the simple models presented here, where material is removed and instantaneously replaced by asthenosphere at a temperature given by the mantle adiabat for a potential temperature of 1,315°C, the mechanism of removal (convective delamination, slab break-off, or the progressive cascade of lithospheric material into the deeper mantle) makes minimal difference to the model outputs, provided that the rate at which removal takes place was fast in comparison to the diffusive cooling time of the material involved (in this case ≤5 Myr). More critical is the depth at which removal takes place, the region over which it occurs, whether it removes significant crustal material from the plateau interior, and the time at which this removal occurs.
The timing of the removal of material from the plateau interior (t b ) can be constrained by the timing of widespread magmatic activity across the central plateau from ∼30 to 10 Ma (Figures 1d and 2c). As discussed elsewhere (Leary et al., 2016;Webb et al., 2017), the removal of material was probably diachronous and may have propagated across the plateau. Similarly, the southward extent of the zone of removal (x b ) can be constrained by the spatial distribution of volcanic rocks (Figure 2d), although these locations must be corrected for the southwards motion relative to the Himalayan front since emplacement. Mantle-derived volcanism predominantly cuts off ∼250 km north of the Himalayan front (although a small number of instances occur further south) and ended by ∼12 Ma. Ascent times for erupted magmas are unknown but are unlikely to be >1 Myr (Mutch et al., 2019). We hence determine that mantle melting must have been possible into the later Figure 6. Thermal models where different thicknesses of the Indian upper crust are transferred into the Himalayan prism. All other parameters are as in Figure 3, and discussed in the text. Tibetan crustal thickness is maintained at 75 km, with differing amounts of accretion compensated for by varying the thickness of the Tibetan overthrust. Left-hand column shown thermal structure, right-hand column shows temperature differences relative to the model in Figure 3.

Geochemistry, Geophysics, Geosystems
Miocene within 300 km of the Himalayan front. The model shown in Figure 3 uses values of t b ¼17.5 Ma and x b ¼ 250 km. In Figures S5 and  S6, we show that the effect on the model thermal structure of varying these values in the range 12.5 < t b < 17.5 Ma and 200 < x b < 350 km.
We also note the uncertainty in the structural geometry of the collision zone prior to removal, and the potential for there to have been more than one phase of material removal. In both of these cases, such complexity would most probably predominantly influence the crustal geotherm in Tibet used as the northern boundary condition, has little impact on the thermal structure from the final phase of removal onwards, with the principal alteration being in the temperatures of the shallow crustal material within the plateau.
The P-T values for the the Namring xenoliths provide the major constraints on the depth to which removal takes place, resulting in the presence of hot mantle material to shallow depths. In the model shown in Figure 3, we include the removal of the lower 10 km of the underthrust Indian crust, allowing mantle material at near-asthenospheric temperatures to reach 65 km beneath the plateau, allowing the fit to the xenolith data shown in Figure 3d. However, whether that is the result of the removal of any crustal material or instead a product of a lesser crustal thickness of the plateau in the Miocene is unconstrained by our modeling. Figure 3 shows one thermal model for southern Tibet consistent with all available required constraints. As discussed, in this model, all Indian mantle and the lower 10 km of the Indian crust north of 250 km north of the Himalayan front are removed at 17.5 Ma and replaced with adiabatic mantle (Figure 3b). The top 15 km of the Indian crust is transferred into the Himalayan prism during underthrusting. Figure 3c shows the modeled thermal structure for the present day, for comparison with geophysical observables. Figures 3e and 3f show a series of temperature profiles for three present-day distances from the plateau margin, illustrating the temperature evolution through time, for comparison with geochemical and geophysical proxies. We will use this model as the basis for the following discussion, but emphasize that it is only one of a range of possible models which are capable of fitting much, if not all, of the available temperature information beneath southern Tibet.

Discussion
The overall "shape" of the thermal structure shown in Figure 3 is common to all our models and to previous studies. Figure 3 represents one scenario capable of appropriately fitting available constraints on the absolute temperature values. The variability in absolute temperature values shown in Figures 5 and 6 serves to illustrate the degree to which certain areas of the model are unknown-in particular, the temperatures present in the Tibetan midcrustal section, around 30-to 40-km depth. The only direct constraints in temperatures at this depth available for southern Tibet are those from estimates of the ABQT, and the indirect impact that changes in midcrustal temperatures have on the better-constrained temperatures in the shallow crust, and near the Moho. As such, the midcrust may vary by up to 200°C from the values shown in Figure 3, while still being able to fit much of the data. Such variation could result most plausibly from variations in accretion (Figure 6), input thermal structure ( Figure 5), or the internal heat production to Tibetan material ( Figure S3). Temperatures in the shallow crust and around the Moho are much better constrains from the geophysical and geochemical data available. We assess that, Geochemistry, Geophysics, Geosystems within the framework of our models as described above, these are unlikely to vary by more than 50-100°C, without leading to significant discrepancies in matching the constraints shown in Figure 2.

Comparison to Previous Models of the Thermal Structure of Southern Tibet
Several previous studies have constructed numerical models for the thermal structure of southern Tibet (e.g., Bollinger et al., 2006;Craig et al., 2012;Henry et al., 1997;Herman et al., 2010;Hétenyi et al., 2007;P. Nábělek & Nábělek, 2014;Wang et al., 2013;Whipp Jr. et al., 2007). The geometry of these models are all, to first order, similar, with India underthrusting beneath the Himalayas, and, in most cases where models extend sufficiently far northwards, flattening out beneath southern Tibet. Depending on the focus of the study in question (typically shallow thermochronometry in the Himalayas, Bollinger et al., 2006;Henry et al., 1997;Herman et al., 2010;Whipp Jr. et al., 2007; or deeper thermal structure beneath the plateau, Craig et al., 2012;Hétenyi et al., 2007;Wang et al., 2013), slightly different approaches to the geometry of the main interface beneath the Himalayas and the kinematics of crustal accretion during underthrusting are incorporated, although, as we note above, these have little impact on the deeper thermal structure. Similarly, as discussed above, the inclusion of a shear (or strain) heating term in Herman et al. (2010) and P. Nábělek and Nábělek (2014), while potentially playing a role in defining the shallow thermal structure of the Himalayas, has little effect on the thermal structure on Indian crustal and mantle beneath southern Tibet.
All of these prior models focus on the present-day thermal structure and, in many case, assume it is in steady state. As a result, the predictions of our model with regards to the past thermal history are not included in this comparison.
The most significant differences between past models, and those that we present here, arise from variations in the thermal structure of incoming Indian lithosphere, and the model kinematics, particularly the rate of underthrusting and the amount of material removed from the top of the Indian plate. As discussed above, the input thermal structures used here are built on data from several regions of peninsula India south of the collision zone. That used in the calculations in Figure 3 is based on crustal heat production data and heat flow from western Dharwar, along with mantle xenoliths from beneath the Dharwar craton. This differs from previous input models, in part due to the inclusion of temperature-dependent material properties and in part due to different assumptions about crustal radiogenic heat production.
The model of P. Nábělek and Nábělek (2014) uses an Indian thermal structure with 100-km-thick lithosphere, and is therefore significantly hotter than any of our models or the other models considered. As a result, while the same overall shape emerges, with a tongue of colder Indian material thrust beneath the plateau, underthrust material does not persist and <700°C beyond ∼200 km from the Himalayan front. Such thin, hot lithosphere from incoming India, while potentially valid for material already thrust beneath the plateau, is inconsistent with measurements of the thermal structure of the cratonic areas of the Indian peninsula, or with its rheology, and unable to fit the majority of the constraints we describe in sections 2 and 3. Similar issues apply to the model of Henry et al. (1997), which also uses a 100-km-thick India lithosphere, as well as relatively high (2.5 μW m −3 ) upper crustal heat production, leading to a very hot incoming Indian lithosphere.
The models of Bollinger et al. (2006), Hétenyi et al. (2007), Herman et al. (2010, and Whipp Jr. et al. (2007) all use a restricted model domain, typically only including the upper ∼80 km of the collision zone, and as such are not appropriate for considering the deeper thermal structure, close to the lower boundary condition at the base of their models (60-to 80-km depth). However, for the crustal section, Bollinger et al. (2006), Whipp Jr. et al. (2007, and Herman et al. (2010) use extremely cold geotherms, most closely comparable to our model using a geotherm from Bundelkhand ( Figure 5, lowermost panels). While that of Herman et al. (2010) is determined through the inversion of thermochronological data from the Himalayas, as discussed above, a geotherm this cold for incoming India is highly unlikely to be representative of much of incoming Indian lithosphere, as it would produce an extremely strong incoming plate, seismogenic to depths of ≥60 km, inconsistent with estimates of the seismogenic and elastic thickness in India. The thermal structure from either model would (unless limited by the removal of material beneath the plateau) lead to the persistence of temperatures ≤600°C significantly further beneath the plateau that in all our models, bar the Bundelkhand model from Figure 5, and would also lead to a decrease in midcrustal temperatures, unless countered by even higher radiogenic heating in Tibetan material than used in our models.
Most similar to our models, especially that shown in Figure 3 are those of Hétenyi et al. (2007) and Craig et al. (2012). Hétenyi et al. (2007) use a similar input geothermal for incoming India to the models for the Dharwar craton we use (Figures 4 and 5), but, as with Bollinger et al. (2006) and Herman et al. (2010), they use a restricted model domain unsuitable for studying the deeper thermal structure of southern Tibet, with a constant heat flow boundary condition imposed at 80 km depth. In contrast, the models presented here demonstrate that the temperature and heat flux at a depth of 80 km vary in space and time. The models presented here are broadly similar to that shown in Craig et al. (2012), updated with improved data on crustal radiogenic heat production, model kinematics, and the thermal parameters of crustal rocks, all described above.

Fit to Temperature Proxies Beneath Tibet
Our model is capable of matching the distance distribution of deeper seismicity beneath India and Tibet and of shallow seismicity across the plateau, with an aseismic region between, based on the expected temperature cut-offs discussed above. The emergence of the region of inefficient head-wave propagation beyond ∼500 km north of the Himalayan front fits with the transition to a negative temperature gradient across the Moho at this distance (right-hand side, Figures 3c and 3f), although the nature and origin of the lower under most of this region is not encompassed by our model. Similarly, the development of a consistent anisotropic fabric at a similar distance fits with the progressive northwards heating and weakening of the crustal material, and the development of significant ductile deformation within the midcrust, as temperatures exceed 750°C, although we note that without knowledge of the mineralogy and stress state, it is not possible to make a quantitative comparison between our model results and the anisotropy data.
As the present-day temperature profiles (red lines in Figures 3d-3f) show, our model does a good job of fitting the estimated depth of the Curie isotherm across Tibet (Alsdorf & Nelson, 1999), fitting to within ∼50°C or a 1-km error in depth. Fits for the ABQT depths for southern Tibet from Mechie et al. (2004) (Figures 3e  and 3f) similarly are well explained with our models. Fitting the ABQT depth of Sheehan et al. (2014) in the Himalayan midcrust is more problematic (see Figures 3d and 3e), with the implied high midcrustal temperatures failing to match our models by ∼250°C, or ∼100°C and a 5-km decrease in depth. Fitting the pressure and temperature estimates from the ABQT at this location requires either extreme amounts of radiogenic heat production within the plateau, implausibly rapid erosion and exhumation of material through the Himalayas (and therefore likely incompatible with shallow thermochronology Herman et al., 2010), or potentially the inclusion of shear heating along the interface, which would boost the temperature along the plate contact and along narrow zone above the interface in the Himalayan prism (P. Nábělek & Nábělek, 2014).
We now consider the fit to geochemical data across southern Tibet. In comparing our thermal models to geochemical data (Figures 3d-3f), we recalculate the position of geochemical samples back through time based on the movement of material within the plateau toward India as a result of the erosive removal of material through the Himalayan front, using the integrated erosive component of the advective velocity field between the present and the time of emplacement of each sample. Profiles prior to 0 Ma on Figures 3e and 3f therefore do not represent the distance quoted, but the greater distance representative for that point at each period in time.
Our model fits the constraints derived from north Tibetan xenoliths from Taipinghu (Figure 3f), although the linear geotherm assumed for the northern boundary condition is not entirely appropriate. As a result of the removal of material from beneath the plateau at 17.5 Ma, our model reproduces the distribution and latest occurrence of mantle melting beneath the plateau. As shown by the 15-Ma profile on Figure 3d, our model also matches the P-T-t estimates from the Namring xenoliths. We note that, in our modeling approach, matching both the present-day crustal thickness of the plateau, the P and T values for these xenoliths and a mantle source, requires the removal of the lowermost Indian crust during break-off, but this could equally be accommodated by having a thinner crust beneath the plateau prior to the removal of mantle material, and neither possibility can be discounted. A removal age of 17.5 Ma is best able to fit the Namring xenolith data, but this remains uncertain to within a few million years. The time at which this removal of material occurs has little effect on the final thermal structure, due to subsequent re-underthrusting (see Figure S5), and trades off with the postremoval convergence rate (see Figure 7).

Geochemistry, Geophysics, Geosystems
As discussed earlier, the lack of age constraints for the P-T estimates from xenoliths at Sailipu and Yibuchaka, along with the uncertain nature of pressure estimates based purely on the lack of garnet, mean they place only limited constraints on the deeper thermal structure of the plateau. However, we note that for for Sailipu (Figure 3e), our model would come close to replicating the required P-T conditions shortly prior to their emplacement near the surface at 17 Ma. However, a complete match to the Sailipu or Yibuchaka xenolith P/T ranges (shown in Figure 3e) would require the removal of a greater thickness of the lower crust than in the model shown in Figure 3, in order to get the required temperatures at a sufficiently shallow depth to be above the garnet-spinel transition.
Overall, the model shown in Figure 3 is capable of matching both the majority of the geophysical constraints available on the present-day structure of the plateau and the geochemical constraints on the plateau evolution in the early-mid Miocene. The difficulty of matching both the geophysical data on the modern plateau and the geochemical data for the prior thermal structure of the plateau requires two critical tectonic events: the removal of significant volumes of mantle material to depths close to the Moho either at or continuously until the mid Miocene and the re-underthrusting of colder Indian material beneath southern Tibet since that removal. To fit the modern-day geophysical structure of southern Tibet, this underthrusting is required to emplace a cold, subhorizontal Indian slab beneath southern Tibet to a minimum of ∼500 km north of the Himalayan front. Given that the tip of this slab will also be heated by lateral diffusion at its end, its true extent is likely greater.
For the tip of the underthrust Indian lithosphere to reach such distances following removal of lithosphere beyond ∼250 km north of the Himalayan front, while also having break-off at a sufficient distance and time to match the occurrence of high Mg# volcanism, requires that the average rate of underthrusting (i.e., the motion of India relative to central Tibet) since break-off was faster than it is at present, although the rate of underthrusting trades off against the time and southernmost extent of break-off. Figure 7 plots convergence rates and penetration distances of underthrust Indian material for a series of convergence models, demonstrating that models (i) with removal extending farther south than that shown in Figure 3, (ii) ending later than in Figure 3, or (iii) without acceleration after removal would fail to match the present-day northwards extent of India beneath southern and central Tibet. This is subject to the caveat that our model assumes uniformity along strike in terms of both the break-off distance and the penetration distance of present-day India beneath the plateau, and some variability in these values could potentially reduce the required slow down in convergence rates, while still fitting available data constraints.
Underthrusting since the late Miocene at a rate greater than the current rate at which Indian material is thrust beneath Tibet is conceptually in alignment with the reduction in resistive force likely to result from the removal of stalled Indian mantle beneath central/southern Tibet. Additionally, deformation in the distal regions of the India-Asia collision zone (e.g., Tien Shan, Qilian Shan, Altai, and Baikal Rift) either initiated or increased over the period following break-off (Bullen et al., 2001;Lease et al., 2011;Molnar & Stock, 2009;Petit & Déverchère, 2006). Given that the total convergence rate varied relatively little over this time period (Copley et al., 2010;Molnar & Stock, 2009), this deformation requires an implicit decrease in the rates of deformation along the Himalayan front over the last ∼15 Ma.

Conclusion
Simple models for the temperature evolution of southern Tibet are able to reconcile all available geophysical and geochemical information on the deep thermal structure of the plateau through time. Three key features from these models emerge. First, present-day observations require the presence of a cold Indian plate emplaced beneath southern Tibet by large-scale underthrusting. Second, geochemical evidence requires the removal of Indian-derived mantle (and potentially lower-most crust) from beneath the plateau during the early Miocene. Our models are unable to resolve whether crustal removal is also required, or if crustal thickness of the plateau has varied through time. Lastly, to match both the southwards extent of Miocene mantle-derived volcanism, and the present-day geophysical structure of the plateau, the average rate since the middle Miocene at which India underthrust Tibet must have been faster than it is at present.

Data Availability Statement
All data used in this study have previously been published and are available as described in the text.