Geodynamic Modeling of Lithospheric Removal and Surface Deformation: Application to Intraplate Uplift in Central Mongolia

Intraplate surface deformation is enigmatic and the underlying mechanisms responsible are not fully understood. We use thermo‐mechanical numerical modeling to explore the conditions under which lithospheric removal processes can occur and investigate the timing and amplitude of consequent surface deformation. We allow lithospheric removal to develop dynamically and self‐consistently by applying a phase transition and density jump that is hypothesized to be a consequence of metamorphic eclogitization in a thickened crust, rather than simply imposing an initial dense block. By systematically varying parameters, we test their influence and control on lithospheric removal. We confirm that a weak and dense lower crust, a hot crust‐mantle boundary, and a convergent regime are critical requirements for delamination‐style removal. We apply these results to central Mongolia, which is an ideal natural laboratory for studying intraplate surface uplift because of its high topography and location in the continental interior. We determine that the physical parameters and structure inferred in this region match the conditions for delamination. Thus, model outputs are evaluated against the available observational evidence. We find that the removal of the lithosphere by delamination can explain: the dome‐shaped topographic pattern and elevated surface; the thin lithosphere and its structure; the elevated temperature at the crust‐mantle boundary. Additionally, it is hypothesized that delamination leads to magmatism due to mantle decompression melting, consistent with intraplate volcanism. Ultimately, the results suggest that lithospheric removal by delamination is a physically plausible mechanism and a potential explanation for intraplate uplift.

. A simplified schematic illustrating the concept of lithospheric delamination-induced surface uplift. (a) Continental crust in a convergent regime (b) thickens, and a dense root may form (e.g., by eclogitization; bright yellow). (c) A weak lower crust allows a portion of the lithosphere to decouple, peel away, and be removed, sinking into the asthenosphere below, which is concurrent with rapid surface uplift. Mantle thermal buoyancy must also be considered. (d) A portion of the lithosphere is replaced with a buoyant and hot asthenospheric upwelling that supports surface uplift. Ultimately, cooling leads to renewal of the lithospheric mantle.
This study investigates the potential origin of intraplate deformation and its link to dynamic crust-mantle interactions by modeling a lithospheric removal event and exploring the conditions under which it can occur. The results are evaluated with respect to the regional setting in Mongolia to test whether they can offer a plausible explanation for the observed intraplate uplift.

Regional Background: Central Mongolia
Central Mongolia is an ideal natural laboratory for studying intraplate uplift and surface deformation processes because it is a unique region with high topography located far into the continental interior, several thousand kilometers from tectonic margins ( Figure 2). In contrast, other regions of intraplate deformation may not be completely isolated from the effects of nearby tectonic plate boundaries (e.g., Dyksterhuis & Müller, 2008), presenting complications. Central Mongolia (approximately latitudes of 45°N to 50°N and longitudes of 95°E to 105°E) is dominated by the Hangai Dome and plateau, an uplifted, high elevation, low relief region that is ∼1-1.5 km above the regional elevation (with local peaks up to 4,000 m above sea level), which covers a large area (∼200,000 km 2 ; Cunningham, 2001). In the south, an internally drained basin, known as the Valley of Lakes, lies between the Hangai and Gobi-Altai mountains (e.g., Cunningham 2001   (a) Map of the regional tectonic setting of Mongolia, which is located between the Siberian Craton and the Tibetan Plateau. Central Mongolia is dominated by a high elevation dome (Hangai) that is bounded by major fault zones. The Bulnay and Gobi-Altai/ Bogd fault zones, as well as the South Hangai fault system (SHF), are marked (gray lines; e.g., Calais et al., 2003;Walker et al., 2007). The dotted line marks the approximate location of the Mongol-Okhotsk suture zone (MOS; Van der Voo et al., 2015). GPS vectors (arrows) display northward motion in northern China and eastward motion in central Mongolia . Bold line is section shown in (c). (b) Map of geophysical data collected across Mongolia showing locations of magnetotelluric measurements (MT, red;Käufl et al., 2020) and seismic measurements (blue; Meltzer et al., 2019). Dashed yellow line denotes the −250 mGal Bouguer anomaly . Gray lines are faults. (c) Smoothed elevation (above sea-level) on a ∼1,000 km profile along a longitude of 100.5°E. The elevated Hangai Dome is visible. Fault zones are identified with gray bands. (d) A simplified schematic of the subsurface structure, based on geophysical, geochemical, and geological measurements and models (modified from Comeau et al., 2017). Beneath the Hangai Dome, evidence indicates a thin lithosphere, a small-scale asthenospheric upwelling, shallow melting, and a weak lower crust.

Tectonic Setting and Geology
Central Mongolia is part of the Central Asian Orogenic Belt (CAOB), a long-lived accretionary orogen that covers a large area of Central and Eastern Asia (e.g., Yin, 2010). The CAOB is composed of several Precambrian to Triassic (lithotectonic) units that are distinguished based on their tectonics and lithological successions, and the formation of these different units is interconnected. Furthermore, they are known to have diverse geophysical signatures (e.g., Guy et al., 2014; see also Comeau, Becken, Käufl, et al., 2020). In particular, the region beneath the Hangai Dome consists of an ancient (Precambrian) micro-continental basement block (Badarch et al., 2002). The CAOB occupies a unique position, to the north is the Siberian craton, which is relatively stable (e.g., Calais et al., 2003), and to the south are the North China and Tarim cratons, with the whole region having sustained periods of (northward) convergent motion during Phanerozoic continental growth (e.g., Jiang et al., 2017;Yin, 2010). The closure of paleo-oceans in Mongolia during this time (e.g., the Paleo-Asian Ocean and the Mongol-Okhotsk Ocean; e.g., Buchan et al., 2001;Jiang et al., 2017;Van der Voo et al., 2015) likely introduced fluids and volatiles into the lithosphere beneath Mongolia.
Central Mongolia is bounded by large, seismically active, strike-slip faults (e.g., Calais et al., 2003;Walker et al., 2007) (Figure 2). The distribution of these faults across the landscape is striking, with long ruptures and large offsets visible (Rizza et al., 2015). Notably, both the left-lateral Bulnay (north) and Gobi-Altai/ Bogd (south) fault systems have experienced several intracontinental earthquakes larger than magnitude 8 within the last century (e.g., Rizza et al., 2011;. South of the Hangai Dome, a long shear zone is identified (the South Hangai Fault system, SHF; e.g., Walker et al., 2007), which shows a considerable amount of small-scale seismicity (Welkey et al., 2018). This fault system does not result from stresses introduced by the locally elevated topography (Walker et al., 2007), rather it is believed to be a reactivated shear fault system that is associated with a pre-existing lithospheric weak zone attributed to an ancient (Cambrian) collisional suture zone (e.g., Buchan et al., 2001), which is crustal-scale or possibly lithospheric-scale.
Evidence from GPS measurements indicate that the Gobi-Altai/Bogd fault system accommodates residual northward compressional motion from Northern China of <5 mm/year, and as a consequence the transpressional Gobi-Altai mountains were formed (e.g., Calais et al., 2003;Cunningham, 2001), a large part within the past 10 Myr (Vassallo et al., 2007). In contrast, to the west in adjacent parts of central Asia, measurements indicate up to 15 mm/year of northward-directed motion-highlighting the dominant (horizontal) tectonic forces Yin, 2010). In contrast, the Hangai Dome appears to be little modified by the regional strain (e.g., Cunningham et al., 2001;Walker et al., 2007).
Widespread Mesozoic and Cenozoic volcanism exists across Mongolia. Well-documented, low-volume, basaltic, intraplate volcanism of Cenozoic age is found dispersed throughout the Hangai region (Barry et al., 2003;Hunt et al., 2012). Mesozoic (including Cretaceous) volcanism is broadly distributed, and is found south of the Hangai region (e.g., throughout the Gobi-Altai and Gobi), and in some regions further east and towards the north (e.g., Sheldrick et al., 2018;Vorontsov et al., 2007;Yarmolyuk et al., 2014). Although often attributed to deep-rooted mantle plumes (e.g., Windley & Allen, 1993), the volcanism has no clear spatial-temporal pattern that would indicate a deep-rooted and fixed mantle-plume (Ancuta et al., 2018;Barry et al., 2003). Interestingly, there are no significant changes in the chemistry of the lava samples across the Hangai region (Cenozoic, majority <15 Myr ago), throughout time or space, which suggests that the erupted magmas may have been generated from a single source within the mantle (Ancuta et al., 2018;Barry et al., 2003;Hunt et al., 2012). In fact, geochemical analysis of erupted material indicates that partial melting occurred over an extremely prolonged period and originated at a depth of ≥70 km, with an asthenospheric signature (Barry et al., 2003;Hunt el al., 2012;Ionov, 2002 Early studies inferred that the Hangai Dome was uplifted in the mid-late Cenozoic and was contemporaneous with the initiation of local (Cenozoic) volcanic activity (e.g., Cunningham, 2001;Sahagian et al., 2016). However, recent evidence from apatite (U-Th)/He thermochronology (McDannell et al., 2018) across the Hangai mountains points to uplift and exhumation that created topography ∼100-120 Myr ago, revealing that the Hangai Dome is an older feature. Geological field observations and geomorphological analysis are consistent with this result (e.g., McDannell et al., 2018). Although it may be that only a part of the current height of the Hangai Dome was created at that time, with modest uplift occurring in the Cenozoic (McDannell et al., 2018), a conjecture supported by basalt vesicle paleo-altimetry (Sahagian et al., 2016).

Geophysical Measurements
Seismic studies indicate that the lithosphere beneath the Hangai region is anomalously thin (60-80 km), compared with its surroundings (∼150-225 km), and detect anomalously low shear-wave velocities at depths shallower than 100 km (Chen et al., 2015;Petit et al., 2008;Priestley et al., 2006). These results may suggest that a previously thicker lithosphere was present beneath Mongolia. In contrast, the crust beneath the Hangai region is found to be thick (50-55 km), and is thin towards its edges (∼40 km; Petit et al., 2008;Welkey et al., 2018). Gravity data have determined that the lithosphere beneath the Hangai region is characterized by a very low negative Bouguer anomaly (<−250 mGal), compared to the surrounding region, and models revealed a localized low-density structure at a depth of 80-125 km .
A large magnetotelluric study across central Mongolia detected a deep (70-150 km), localized low-resistivity zone beneath the Hangai Dome that was postulated to represent a shallow, upwelling asthenosphere (Comeau et al., 2018;Käufl et al., 2020). The magnetotelluric data also detected a lower crustal (30-50 km depth) low-resistivity zone that was hypothesized to be a thermally perturbed region that underwent metamorphic dehydration and devolatization reactions, and thus produced fluids (Comeau et al., 2018;Käufl et al., 2020). These fluids can accumulate in fluid-rich domains and become trapped below the brittle-ductile transition zone-as demonstrated by Comeau, Becken, Connolly, et al. (2020) based on a conceptual hydrodynamic model from Connolly and Podladchikov (2004)-and are compatible with low seismic shear-wave velocities (as compared to the surroundings) in the lower crust determined by Mordvinova et al. (2007).
Critically, this implies a weak lower crust in the Hangai region, because the presence of fluids not only reduces the electrical resistivity but substantially changes the crustal rheology and considerably reduces the viscosity (e.g., Liu & Hasterock, 2016;Rosenberg & Handy, 2005; see also Comeau et al., 2015;. This is in agreement with work by Vergnolle et al. (2003) that determined a significantly reduced viscosity (of several orders of magnitude) in the lower crust, compared with the upper crust, by analyzing post-seismic slip measurements along the major fault zones of central Mongolia. Comeau, Becken, Connolly, et al. (2020) determined that a low viscosity in the lower crust is compatible with the conceptual model of crustal fluids and is consistent with reasonable estimates for crustal properties, including the activation energy, the temperature, and the depth to the brittle-ductile transition.

Physical Equations and Numerical Implementation
In this study, numerical models were generated that simulated thermo-mechanical convection processes within the Earth's lithosphere and upper mantle. This was accomplished by employing the finite element code ASPECT Kronbichler et al., 2012) that is freely available and widely used to investigate geodynamic processes at all scales. Computations were done using ASPECT version 2.1.0 (Bangerth et al., 2019). For simplicity, two-dimensional models were generated. Future work is required to investigate the complexities of three-dimensional models.
Considering an incompressible fluid with infinite Prandtl number, the conservation equations can be formed as described below. The temperature (T) equation is written as where ρ is the density, c p is the heat at a specific pressure, t is the time, v is the velocity vector, k is the thermal conductivity, and H is the internal heating rate. The momentum equation is expressed as where p is the dynamic pressure, η is the viscosity, ε is the strain rate tensor, and g is the gravitational acceleration. The continuity equation is written as and the composition (C) conservation as where C i is a set of advected quantities (compositional fields). The compositional fields are used to track the material throughout the model simulation.
ASPECT employs the finite element method for spatial discretization and a fully implicit time-stepping method. The resulting system of linear equations is solved with a geometric multigrid method applying the generalized minimal residual (GMRES) method as a solver . The free surface is implemented in an Arbitrary Lagrangian-Eulerian (ALE) framework (Rose et al., 2017).
The code was adapted so that it can handle both viscosity variations and phase transitions simultaneously. The incompressible Boussinesq equations are considered, in which the density ρ depends on the temperature T and the change in density Δρ, due to a phase transition, in the form where ρ 0 is the reference density (3,250 kg/m 3 ), T 0 is the reference temperature (273 K), α is the thermal expansion coefficient (3 × 10 −5 1⁄K), and Γ represent the phase function. The phase function for a depth z and temperature T is written as where z t is the transition depth (35 km), T t is the transition temperature (1110 K), d is the width of the phase transition (10 km), and  is the Clapeyron slope of the phase transition (1 MPa/K) (values based on van Thienen, 2003; see also Hacker et al., 2003).
The visco-plastic rheology follows the equation for dislocation creep, which represents viscous flow below a yield stress, and is written for a viscosity η as where A is the viscosity prefactor term (material constant), f is a scaling factor (used to control the viscosity, see below), n is the stress exponent, ε ii is the square root of the second invariant of the strain rate tensor, E is the activation energy, V is the activation volume, p is the pressure, and R is the gas constant. See Section 3.4 for the material properties used in the models.
In order to consider crustal deformation, the viscosity is limited by the use of the yield stress σ y , which is a function of pressure p, cohesion C, and angle of friction Φ, and is the same as the Mohr Coulomb surface, expressed as If the viscous stress  2 ii  generated during deformation exceeds the yield stress, the viscosity is limited to the yield surface, that is This is then used to compute the effective viscosity η eff , in the form Aspect has been benchmarked against other mantle convection codes (e.g., Tosi et al., 2015) and we ran additional test cases featuring a phase transition in isoviscous and thermoviscous convection which we tested against the mantle convection code plaatjes (cf. Stein et al., 2004;. Furthermore, Thieulot (2011) presents benchmark results for Rayleigh-Taylor instabilities.

Model Geometry and Boundary Conditions
In order to generate realistic models, the choice of modeling input parameters was guided by available observational data. The modeling domain had a horizontal extent of 500 km and a vertical extent of 450 km ( Figure 3). The model approximates a north-south transect across the Hangai Dome and central Mongolia along longitude 100.5°E, corresponding to the region between the Bulnay and Gobi-Altai/Bogd faults. We limit modeling to the region around the Hangai Dome, in part due to computational restrictions which limit the domain size and for simplification, and leave future work to explore the complexities of features beyond this region. As an initial condition, the crust was defined to be 40 km thick, with the lower crust from 25-40 km depth, in line with values of average continental crust (e.g., Christensen & Mooney, 1995  lithosphere was initially 100 km thick. The model mesh was rectilinear and consisted of 240 elements horizontally and 240 elements vertically. This ensured a high resolution (∼2 km) and enabled fine structures to be discerned.
Within the crust, a central region (175 km across) was defined to be mechanically weaker, achieved by reducing the initial viscosity (see below), compared with the outer edge regions (each 150 km across). This ensured that the crustal thickening occured within a confined zone. In addition, an adjacent pre-existing weakness (25 km wide) was introduced in the crust, on the right side of the central region (model distance of 325-350 km). This was done in order to focus the crustal deformation. It can be interpreted as a local shear zone, a fault zone, or a broad boundary zone between two regions of differing crustal lithologies, which simulates the conditions and setting in Mongolia (i.e., it could correspond to the South Hangai fault system and suture zone). These subregions were modeled by applying a scaling factor to the viscosity (Equation 7; details in Section 3.4). Otherwise, the viscosity variations throughout the model (Figure 3c; Figure S1) were controlled by the depth, temperature, and strain rate, and they evolved with time.
An influx velocity was prescribed to simulate convergent processes provided by external plate forces occurring outside the model domain. The influx was allowed at the right side only and within the lithosphere (depths of 0-140 km; the left side had no influx). For the reference case, a constant rate of 10 mm/year was used, comparable to estimated rates in southern Mongolia and to shortening rates detected elsewhere (a range of values were tested in Section 4.3). For greater depths, a variable outflux velocity was allowed on both sides to account for the conservation of mass. In all cases, the vertical velocity at the model sides was 0 mm/year. The mechanical boundary condition applied at the top of the model consisted of a free surface that allowed for the generation of topography. Note that, for simplicity, erosion was not included. In Mongolia, long term stability of the low-relief landscape is expected due to the high aridity of the regional climate since the Jurassic, with proposed erosion rates in agreement with average global erosion rates (e.g., 5-20 m/ Myr; McDannell et al., 2018;Portenga & Bierman, 2011). The bottom of the modeling domain had a free slip boundary (closed).
Estimated subsurface temperatures for central Mongolia, derived from geochemical analysis (e.g., Barruol et al., 2008;Ionov, 2002), and average continental values were used to guide the choice of initial temperatures throughout the model. The initial temperature at a depth of 40 km, at the crust-mantle boundary, was set to 1113 K (840°C). At a depth of 70 km, within the uppermost mantle, the temperature was set to 1473 K. These values are within the range of typical continental conditions and are similar to values inferred elsewhere (e.g., Jull & Kelemen, 2001). In addition, the temperature was set to 273 K at the top of the model and 1610 K at the bottom. Thermal boundary conditions included isothermal top and bottom boundaries and insulating side boundaries.

Modeling Eclogitization With a Phase Transition
The lower reaches of the thickened continental crust can undergo prograde metamorphism and form eclogite through a phase transition, which is denser than the original material (e.g., Austrheim et al., 1997; cf. Figure S2). In this scenario a dense crustal root is formed. The density depends on the temperature and pressure conditions and on the quantity of material that has been transformed to eclogite. It is known that eclogitization is promoted by the presence of fluid, specifically water, even in small amounts (<2%; fluid also reduces rock mechanical strength significantly), and thus the hydration state of the lower crust may determine the location and quantity of eclogitization (e.g., Austrheim et al., 1997;Bjørnerud et al., 2002;Leech, 2001). Indeed, a source of hydrous fluids, such as the paleo-subduction of a nearby plate, is considered to be a contributing factor of whether or not a region may undergo lithospheric delamination (e.g., Krystopowicz & Currie, 2013;Leech, 2001).
Metamorphic eclogitization is an important process to consider because it plays a role in triggering lithospheric delamination, for the reason that the eclogitized crustal root-which can be more or less dense than the underlying mantle material-is also weak enough to allow the decoupling of the lithosphere (e.g., Jull & Kelemen, 2001;Krystopowicz & Currie, 2013;Le Pourhiet et al., 2006;Li et al., 2016). Furthermore, it has been demonstrated that an initial small density increase (e.g., ∼7%) in the crustal root will subsequently cause further crustal thickening to focus at that location, and thus lead to a region that is more susceptible to gravitational instability and removal (Krystopowicz & Currie, 2013).
In order to include metamorphic eclogitization in the models, a phase transition is introduced when the appropriate temperature and pressure (i.e., depth) conditions are satisfied during the model evolution. The initial phase transition is chosen to be located within the lower crust, at a depth of 35 km and over a transition width of 10 km (see Equation 6). As the simulation evolves, the location of the phase transition may vary. The phase transition is characterized by a density change or density jump-a valid assumption because of indications that the metamorphic transition to eclogite is a rapid and sudden process (Leech, 2001, and references therein)-with values based on field and laboratory measurements (e.g., Austrheim et al., 1997;Doin & Henry, 2001;Leech, 2001). In the reference model the density jump was 400 kg/m³ (a range of values were tested in Section 4.2, some of which produced similar behavior). Implementing a phase transition allows a dense crustal root to arise naturally, rather than simply imposing an initial dense block. Table 1 lists the parameters and material properties used in the modeling, including the rheological and thermal parameters. Values are based on viscous flow laws derived from laboratory experiments: the upper crust approximates wet quartzite (Gleason & Tullis, 1995), the mantle approximates wet olivine (Karato & Wu, 1993), and the lower crust approximates dry diabase (Mackwell et al., 1998), with some deviations to account for the setting in Mongolia, which is inferred to be weaker and warmer. Similar to other numerical studies (e.g., Heron et al., 2016;Krystopowicz & Currie, 2013), we introduce a scaling factor f that is applied to the viscosity prefactor term A (i.e., A' = A/f; see Equation 7). This allows the initial viscosity of the material to be varied (compared to the assumed viscous flow laws) in certain subregions in order to test the influence of this parameter (see Section 4.4). In the reference model, the crust in the central region uses f = 1, which corresponds to a reduction in viscosity (up to 2 orders of magnitude in the lower crust), compared to the edge regions.

Outcome of the Reference Model
Numerical models were generated to show the dynamical evolution of a system using the parameters described above (i.e., the reference parameters). The evolution of the reference model (Figure 4) a The factor f is applied to the viscosity prefactor term A. The values shown are for the outer edges; in the reference model the crust in the central region uses f = 1.

Table 1
The Parameters and Material Properties Used in the Modeling Figure 4. The evolution of the reference model from the numerical simulation, including surface topography, shown at model times of (a-f) 1.6, 3.0, 4.0, 4.7, 5.0, and 5.5 Myr, respectively. Color denotes temperature and arrows represent the velocity field. Dark colors highlight the initially defined subregions (central lower crust; crust-mantle boundary; lithosphere-asthenosphere boundary), plotted as a tracer cloud for each material phase. Delamination-style lithospheric removal is observed and topography is generated with a broad, dome-shaped pattern. Lithospheric thinning occurs across the central region of the model, in a window <150 km wide, whereas the surrounding regions and outer edges have nearly unchanged lithospheric thickness (∼140 km), creating a doming lithosphere-asthenosphere boundary.
begins to sink (causing return mantle flow), subsequently triggers detachment of the lower crust (which is weak and dense, due to the presence of eclogite in the crustal root), and as the detachment point migrates laterally a slab of material peels off (delaminating, and, at much later times, breaking off completely and sinking) and is replaced with an asthenospheric upwelling, leading to dramatic lithospheric thinning in an interval of ∼5 Myr across the entire central region.
In the reference model, the crust and lithosphere within the central region are observed to thicken within the first ∼2 Myr (Figure 4a). Initially, this is due to the convergence regime and the motion of material from the right side. The thickening takes place locally near the pre-existing weak zone, which acts to focus the deformation. A small peak is formed in the topography above this region. As the material near the pre-existing weak zone encounters the phase transition, it begins to move downwards and leftwards (at a model time of ∼3-3.5 Myr; Figure 4b); the density jump causes a buckling of the sinking current (near the phase transition depth) so that the crust within the central region is pushed (downwards) by the motion of the buckling material.
The dense root of the crust then begins to peel off (due to negative buoyancy; model time of ∼4-4.5 Myr; Figure 4c), as the detachment point migrates laterally. This is subsequently followed by rapid surface uplift, which produces a small dome of topography above the central region, a topographic peak above the pre-existing weak zone, and a topographic peak at the right boundary due to continued convergence. The asthenosphere has, locally, risen to a minimum depth of ∼80 km. Note that at this point in time, velocities associated with convective processes exceed the velocity of the incoming material due to the convergent regime.
The weak crustal root decouples from the overlying material near the pre-existing weak zone, the initiation point, and the mantle lithosphere then begins to sink downwards into the underlying asthenosphere (model time of ∼4.5-4.7 Myr; Figure 4d), in a slab ∼100 km long. The asthenosphere continues to flow into a small window (∼50 km wide) behind the separated slab and reaches a minimum depth of ∼60 km, thereby replacing the mantle lithosphere with much hotter asthenospheric material. This likely leads to widespread decompression melting and the generation of partial melts at this location (e.g., Wang & Currie, 2015). The surface uplift continues and the topographic dome grows, supported by the mantle upwelling.
The window widens and shallows and the mantle downflow follows the detached, sinking slab in a clockwise pattern on the right side. The consequent upflow additionally causes an anticlockwise motion on the left side that acts to push the mantle lithosphere downwards (model time of ∼4.7-5 Myr). Similar behavior has been described elsewhere, for example by Ueda et al. (2012) for tectonic evolution in a post-collisional setting. This is followed by a second, smaller, event as another slab of material begins to sink (model time of ∼5-5.5 Myr; Figure 4e).
By a model time of ∼5.5 Myr ( Figure 4f) the mantle lithosphere has been nearly completely removed in the central region, in a 100-150 km wide window. Thus, at this time the lithosphere is locally the thickness of the crust (∼40 km thick). However, the lithosphere is more than 43-51 km thick at a model time of 7-10 Myr ( Figure S3), which corresponds to up to ∼10 km of lithospheric mantle. In contrast, the surrounding regions (e.g., the outer edges) have practically unchanged lithospheric thickness (∼140 km), creating a doming lithosphere-asthenosphere boundary. With enough time the sinking slabs break off and are completely removed; the small slab (left) first and the large slab (right) much later, after reaching a depth of ∼300 km (by a model time of ∼8 Myr and 20-28 Myr, respectively; not shown).
The pattern of topography produced is asymmetrical, with a broad dome of topography on the left side and a sharp narrow peak on the right side, with the surroundings having, generally, a lower and smoother topography. Such an asymmetric pattern is in contrast to the highly symmetric pattern predicted from modeling lithospheric drips caused by a Rayleigh-Taylor instability (e.g., Göğüş et al., 2017;Göğüş & Pysklywec, 2008). Topography develops over the course of the simulation, starting at a model time of ∼2 Myr. There are different ways to measure the amplitude of surface deformation; here we consider the prominence, the difference between local maximum and minimum topography (the zero line is not relevant in this case). For the dome, the maximum amplitude of the topography is 600-800 m measured from the top of the dome to the surroundings on the left and 1,600-1,700 m on the right (at a model time of 5-5.5 Myr). For the peak, the amplitude is 2-3 km, although at later times it decreases ( Figure S4) In order to test the influence of features at the model boundary on model evolution a special case was set up that included a second weak zone near the right boundary. This was inspired by the setting in Mongolia, that is, it could correspond to a fault system or tectonic block boundary (e.g., the Gobi-Altai/Bogd fault system). The special case showed very similar behavior to the reference model ( Figure S5). However, important differences were the start time of the delamination event, delayed by 1-2 Myr, and a topographic peak generated near the second weak zone. In general, this special case validates that the behavior in the central region is not dependent upon unique model boundary features.

Varying the Density Change
Metamorphic eclogitization of the lower crust coincides with a large density increase (e.g., 10%-17%; Leech, 2001). This effect is included in the models by implementing a phase transition. The average density of the crustal root depends on the proportion that has been transformed to eclogite, which may not be all the material within the eclogite regime (i.e., incomplete metamorphism), although measurements suggesting 40%-90% eclogitization are typical (Bjørnerud et al., 2002). Evidence shows that metamorphic eclogitization may develop in localized, heterogenous zones, on scales below the resolution of the numerical models (e.g., Austrheim et al., 1997;Jin et al., 2001;Leech, 2001).
It is assumed that an average lower crustal density change can be representative; for example, an overall ∼10% density increase can be achieved if approximately half of the total material is converted to eclogite (using typical values for mafic granulite, ∼3,000 kg/m 3 , and mafic eclogite, ∼3,600 kg/m 3 ; e.g., Christensen & Mooney, 1995; see also Austrheim et al., 1997;Doin & Henry, 2001;Leech, 2001), in line with other numerical studies (e.g., Krystopowicz & Currie, 2013). Note that the eclogitized crustal root may be less dense or more dense than the underlying mantle material (e.g., Krystopowicz & Currie, 2013). The reference model uses a density change (Δρ) of 400 kg/m 3 (e.g., corresponding to eclogitization of at least 67% of the material (0.67*[3,600-3,000 kg/m 3 ]). The effect of values over a range of 0-700 kg/m 3 were tested (i.e., from no eclogitization to full eclogitization).
All simulations showed removal of lithospheric material, but the style of removal differed dramatically.
No delamination-style lithospheric removal is observed if no density change nor eclogitization is included in the models, that is, a density change of 0 kg/m 3 (Figure 5a). In this case, only crustal thickening from COMEAU ET AL.

10.1029/2020JB021304
12 of 26 continental convergence is observed and the thickened lithosphere undergoes a Rayleigh-Taylor instability with a stationary detachment location, causing material to sink downwards. Note that the removed material has a distinct behavior and a different velocity pattern than the reference model, due to the different physical mechanisms causing the motion (e.g., compare the velocity arrows in Figures 4 and 5). As a consequence of material removal, moderate surface uplift is exhibited above the unstable dripping regions.
For a density change as low as 100 kg/m 3 delamination-style behavior is observed. However, for a density change of 200 kg/m 3 ( Figure 5b) and above (corresponding to eclogitization of at least 33% of the material), larger amplitude and longer wavelength surface uplift is observed. Models with a density change of 300-500 kg/m 3 showed similar evolution with little change, including a somewhat constant maximum amplitude of the topography and little variation in the initiation times and event durations. This was also true for higher values of 600-700 kg/m 3 , although they may be too large to represent physically possible values for eclogitization (Figure 5c). Although the model evolution appears similar for a range of values, the physical mechanisms responsible are different. In the case of a low density change, the downward motion is attributed to a ridge-push phenomenon, whereas in the case of a high density change the circumstances of the downward motion are similar to a slab-pull scenario.

Varying the Convergence Rate
To simulate convergent processes occurring outside the model domain, a velocity influx was prescribed. The reference model used a convergence rate (v in ) of 10 mm/year and values over a range of 0-70 mm/year were tested. No lithospheric removal of any kind is observed if no continental convergence is included in the models, that is, a convergence rate of 0 mm/year (Figure 6a). This will be true at all times, because, in this case, there are no convective processes that occur and no deformation. The thicknesses of the crust and lithosphere remain practically unchanged, as does the surface topography. For a convergence rate of 5 mm/ year delamination-style lithospheric removal is observed, with smooth and broad topography generated (Figure 6b). Values of 10-25 mm/year showed similar model evolution to the reference case, and the maximum amplitude of the topography was relatively unchanged.
COMEAU ET AL.

10.1029/2020JB021304
13 of 26 For higher convergence rates the topography was rougher (i.e., more spikes) and occurred over a broader area. This may be a consequence of folding, with the horizontal flow in the surface layer, which is decoupled from the convection below, dominating over the vertical flow. In addition, the higher values (i.e., faster convergence) had an earlier initiation of the first removal event (up to 1 Myr earlier) and a shorter duration of the event (up to 2 Myr shorter). Although, the duration of the secondary removal event was unaffected. A rate of 30 mm/year resulted in (large amplitude) undulating topography with no distinct dome shape, and no large peak above the weak zone (Figure 6c). It was observed that if the rate is too high (e.g., 70 mm/year), the models exhibit a type of restrained delamination. This is possibly because the dense crustal root can be removed piecemeal before there is enough material present for delamination to occur.

Varying the Viscosity of the Crust in the Central Region
The crust in the central region was defined to be mechanically weaker than the outer edges. To do this a scaling factor f was applied to the viscosity (Equation 7). The reference model used a factor of f = 1 for the weak region, which corresponded to a viscosity decrease of up to 2 orders of magnitude compared with the case of no weak region (∼10 19 -10 20 Pas compared to 10 21 -10 22 Pas in the lower crust; Figure S6). The factor was varied and tested from a value of f = 10 7 , which was similar to no weak region (i.e., the outer edges), to a value of f = 10 −8 , which had a viscosity up to 1.5 orders of magnitude lower than the reference case.
Delamination-style lithospheric removal was not observed if a weak crust was not included (Figure 7b), and crustal thickening was also limited. Delamination was also not observed for a slightly weakened crust (Figure 7), for example, for values of f ≥ 10 4 , which corresponded to a viscosity decrease of up to 1 order of magnitude (∼10 20 -10 21 Pas in the lower crust). In these cases, Rayleigh-Taylor instabilities developed and caused material to be removed. However, the process occurred over a longer time period than in the reference case (>9 Myr), and material did not detach completely and break off (within the modeled time). Furthermore, the lower crust was largely unaffected. In addition, the topography did not show a broad and pronounced dome-shaped pattern. For values of f ≤ 10 3 , delamination was observed to proceed as in the reference case. However, for values of f ≤ 10 −4 , which represent a very weak crust, the topographic increase was diminished for both the dome and the peak. Therefore, the results indicate that a moderately weak crust, relative to its surroundings, is required to achieve delamination-style removal with an elevated and dome-shaped topography.

Varying the Viscosity of the Pre-Existing Weak Zone (Initiation Point)
A pre-existing weakness, introduced on the right side of the central region, acts as an initiation point for detachment because it focuses deformation. It can be interpreted as a local shear zone, a fault, or a boundary between different crustal lithologies. To control its strength a scaling factor f was applied to the viscosity. The reference model used a factor of f = 1 for the weak zone, which corresponded to a viscosity decrease of 1.5-3.5 orders of magnitude in the lower crust compared with the case of no weak zone (∼2 × 10 18 Pas compared to 10 20 -10 22 Pas; Figure S7). The factor was varied and tested over a range of values, from f = 1 to f = 10 18 , which was similar to no pre-existing weakness (i.e., the outer edges).
Delamination was not observed if a pre-existing weak zone was not included in the model (Figure 8b). Delamination was also not observed if a slightly weakened zone was included (Figure 8), for example, for values of f ≥ 10 12 , which corresponded to a viscosity decrease of 0.5-1 orders of magnitude in the lower crust. In these cases, the crust exhibited a buckling behavior, and subsequently generated clear Rayleigh-Taylor instabilities that caused large drips of material to be removed. The initiation of this process was delayed and the process occurred over a longer time period than in the reference case (>9 Myr). Surface uplift was found to occur, across the entire modeling region, with rough topography and no discernable dome-shaped pattern (and no large peak above the weak zone). For values of f = 10 9 to 10 12 , removal of the lithosphere was observed in only one event (beginning at a model time of [3][4]. For values of f ≤ 10 8 , delamination was observed to proceed as in the reference case, with multiple events. However, for values of f ≤ 10 5 a notable increase was observed in the maximum amplitude of the topography and a dome-shaped pattern was generated

Varying the Lower Crustal Temperature
A high temperature in the lower-most crust has been demonstrated to be important for triggering lithospheric removal (e.g., Jull & Kelemen, 2001;Morency & Doin, 2004). In the reference model the temperature at a depth of 40 km (T 40 ) was 1113 K (840°C). To test the effect of the temperature near the crust-mantle boundary, values over a range of 723-1423 K were tested. For temperatures of 873-1423 K delamination was observed and model evolution was somewhat similar to the reference model. The exception was the topographic pattern: the higher the temperatures the smaller the amplitude of the dome (Figure 9). This may arise because the temperature difference between the buoyant upwelling and the ambient material becomes small.
COMEAU ET AL.

10.1029/2020JB021304
15 of 26 Distinct behavior was observed for lower temperatures of 723-773 K. Note that temperatures of ≤773 K (≤500°C) have been proposed for Archean cratons, including the Siberian craton (Artemieva & Mooney, 2001). In these cases, delamination-style removal was not observed (Figure 9a). Rayleigh-Taylor instabilities caused material to be removed over a longer time period than in the reference case (>9 Myr). Furthermore, a part of the lower crust was considerably unaffected and never fully removed. Therefore, the sinking of the material is observed to be hindered and delayed by the lower temperatures.
COMEAU ET AL.

10.1029/2020JB021304
16 of 26 no weak zone (model times of 5.8 and 9.0 Myr). In both cases, a Rayleigh-Taylor instability occurs; delamination is not observed. The resulting topography is rough, with no discernable dome-shaped pattern.

Discussion and Application to Central Mongolia
In the following section the simulations presented in this study will be evaluated with respect to central Mongolia. This is an intraplate region with elevated, dome-shaped topography located in the continental interior far from tectonic boundaries that also has a thin lithosphere compared to its surroundings (details in Section 2). Therefore, an episode of lithospheric removal by delamination is hypothesized as an explanation for its origin. First, the conditions revealed to be important for delamination processes (summary in Figure 10  Alternatively, other mechanisms that lead to lithospheric removal and a mantle upwelling are possible. Convective removal of drips of material from a stationary location is possible (e.g., Göğüş et al., 2017), as observed in some of the models in this study, depending on the chosen parameters. Note that this process takes much longer to remove material than delamination, as confirmed by Beall et al. (2017). Another possibility is an upwelling that originates at the mantle transition zone related to a stagnant slab (e.g., Chen et al., 2015). In this scenario, an upwelling is caused by piles or fragments of cold subducted oceanic plates ponded at the mantle transition zone sinking into the lower mantle, possibly triggered by thermal perturbation of the upper mantle, for example by changes in the subduction of the Pacific plate (e.g., Dash et al., 2015). A further potential mechanism is edge-driven convection, in which the flow of mantle material is influenced by a lithospheric step (e.g., between Mongolia and the Siberian craton) that can trigger lithospheric erosion and removal. Although the size of the possible upwelling, the amount of lithospheric material that could be removed, and the potential amplitude of the surface uplift from this effect is not fully known.
Structurally, a locally thinned lithosphere is the key signature that would be expected after an event caused lithospheric removal (i.e., thinning), such as a delamination event (Figure 1), illustrated by the modeling results (e.g., Figure 4). Other studies have shown that a thinned lithosphere will remain thin for a long time-possibly on the scale of geological periods-because the asthenosphere, which should cool to form a new lithospheric mantle, may maintain a high temperature due to ongoing convection processes and cratonic insulation (e.g., Gorczyk et al., 2015;Ueda et al., 2012).

Convergent Tectonic Regime
The results of numerical studies, including those presented here, have shown the importance of a convergent tectonic regime to initiate crustal thickening that leads to the triggering of a delamination event.
Although, some studies required a limited convergence rate to avoid suppressing the detachment process (dependent on the mantle strength; Krystopowicz & Currie, 2013), and others required convergence only to start the thickening process (Ueda et al., 2012).
Mongolia has sustained long periods of convergence throughout the Phanerozoic, including from (northward-directed) accretion and paleo-ocean closures (e.g., Jiang et al., 2017;Van der Voo et al., 2015) as well as from the effects of the India-Asia plate collision during the Cenozoic (e.g., Yin, 2010). Present-day GPS measurements in northern China (south of Mongolia) indicate 5 mm/year of northward motion, whereas the Siberian craton (north of Mongolia) is considered to be stable (Figure 2; Calais et al., 2003).

Weak Lower Crust
The simulations illustrate that the existence of weak lower crust is a critical requirement, necessary to decouple a portion of the lithosphere and to trigger removal by delamination in a timely manner. This confirms previous results that determined the importance of this condition (e.g., Krystopowicz & Currie, 2013;Le Pourhiet et al., 2006;Meissner & Mooney, 1998). However, in practice, a weak lower crust is not always straightforward to identify.
In central Mongolia two pieces of evidence point to the existence of a significantly weak lower crust. First, electrical resistivity models from Käufl et al. (2020) and Comeau et al. (2018) clearly detect a broad and pervasive low-resistivity zone in the lower crust beneath the Hangai region. Comeau, Becken, Connolly et al. (2020) determined that this feature can be attributed to hydraulically interconnected domains of fluids, consistent with a conceptual hydrodynamic model from Connolly and Podladchikov (2004). The presence of fluids is known to substantially change the rheology of the crust and to significantly reduce the viscosity, thereby reducing the mechanical strength (e.g., Rosenberg and Handy, 2005; and the electrical resistivity, see; Comeau et al., 2016). Second, post-seismic slip analysis along major faults in central Mongolia by Vergnolle et al. (2003) determined that a significantly reduced viscosity (of several orders of magnitude), compared with the surroundings, was necessary in the lower crust. Furthermore, the presence of eclogitized lower crustal rocks (see below) can also reduce rock mechanical strength (e.g., Leech, 2001).

Dense Eclogite in the Lower Crust
The results agree with previous studies that determined that a dense crustal root, which can be generated as a consequence of metamorphic eclogitization, acts to focus crustal thickening and strongly facilitates the decoupling of the mantle lithosphere and the initiation of a delamination event (e.g., Gray & Pysklywec, 2012;Krystopowicz & Currie, 2013;Li et al., 2016;Meissner & Mooney, 1998). A hydrated lithosphere may promote metamorphic eclogitization, especially in a convergent tectonic regime (e.g., Gray & Pysklywec, 2012;Jin et al., 2001). Past subduction-related events, for example, may cause lithospheric hydration, and could act as pre-conditioning events, even when occurring long before delamination (e.g., ∼100 Myr; Ueda et al., 2012).
Mongolia has a complex tectonic history that includes ocean closures in the Phanerozoic (e.g., Jiang et al., 2017;Van der Voo et al., 2015;Buchan et al., 2001). These include the (latest Jurassic-earliest Cretaceous) closure of the Mongol-Okhotsk ocean, possibly by double-sided subduction, which left a prominent fault and suture zone across Mongolia, including adjacent to the Hangai Dome (e.g., Van der Voo et al., 2015). These events likely introduced hydrous fluids into the lithosphere, and thus they may have caused metamorphic eclogitization. In fact, eclogite is found in samples from nearby suture zones (Štípská et al., 2010), as well as from the northern Hangai region (Perepelov et al., 2020). Occurrences of (Cretaceous) adakitic lavas, often interpreted to be formed by partial melting of eclogite, have been discovered in the eastern Gobi-Altai region , and, in a continental setting, may be coincident with lower-crustal delamination (Xu et al., 2002). Therefore, it is inferred, based on geochemical data, that metasomatism of a hydrated lithosphere occurred beneath Mongolia (e.g., .

Pre-Existing Weak Zone (Initiation Point)
The modeling results demonstrated that a pre-existing weakness in the crust acts to focus deformation, and consequently becomes an initiation point for detachment of the lithosphere. Other studies have demonstrated that inherited crustal structures influence tectonic evolution (e.g., Heron et al., 2016). In intracontinental settings, a pre-existing weakness may be caused by fault zones or post-collisional suture zones. Such zones can be broad regions of deformation and typically undergo intense shearing as well as lithospheric metasomatism, and consequently they possess significantly different rheological and chemical properties from the surrounding regions (e.g., Gorczyk et al., 2013). Therefore, these zones will respond differently if exposed to far-field tectonic stresses, such as readily deforming and focusing strain.
South of the Hangai Dome, an ancient (Cambrian) collisional suture zone is found (e.g., Buchan et al., 2001), which is crustal-scale or possibly lithospheric-scale, that is today associated with a re-activated fault system, the South Hangai fault system (Walker et al., 2007). The structure appears as a pronounced feature in electrical resistivity models and is interpreted to be a weak zone between rheologically distinct tectonic blocks (Comeau et al., 2018(Comeau et al., , 2021Käufl et al., 2020). Comeau et al. (2018) established that the South Hangai fault system and suture zone terminates the weak crust inferred beneath the Hangai region. Similarly, south of the Gobi-Altai/Bogd fault zone the lower crust was found to be distinct, more resistive and inferred to be stronger, and the lithosphere was found to be thicker .

Surface Responses: Surface Deformation
Surface responses to crust-mantle dynamics should be observable. Surface uplift is an instantaneous response to a lithospheric removal event and the subsequent buoyant asthenospheric upwelling (e.g., Kay & Kay, 1993;Magni & Király, 2019). The models display an asymmetric pattern of surface deformation, with topography generated on the order of 1 km. For the reference model (and some other models with moderately different parameters, see Figure 10), the final topography consists of a broad, smooth, dome-shaped pattern, which is located above the asthenospheric upwelling. In addition, a sharp, narrow peak is produced near the pre-existing weak zone. Topographic deformation starts at a model time of ∼2 Myr and develops over the course of the simulation, with the majority of the surface uplift occurring in an interval of ∼5 Myr, which matches the timing of the initial instability and removal of material. The duration (and rate) of surface uplift is similar to findings of some numerical modeling studies (e.g., Göğüş et al., 2017), but is shorter (and faster) than others (e.g., Krystopowicz & Currie, 2013).
The topographic pattern generated in the models is qualitatively similar to that observed in central Mongolia (compare Figures 2 and 4), with an elevated, broad, low-relief dome that is ∼200 km wide. There may be some deviation between the amplitude of the modeled topography and that observed in central Mongolia (which is ∼1-1.5 km). However, it has been hypothesized that the present-day topography may have been generated during more than one period (e.g., the surface may have been uplifted by up to 1 km over the past ∼30 Myr; McDannell et al., 2018;Sahagian et al., 2016). Furthermore, in agreement with the models, exhumation rates predicted for the Hangai Dome (McDannell et al., 2018) indicate that ∼0.5-1.5 km of surface uplift would be possible within ∼5 Myr.

Surface Responses: Increase of Volcanism
It is predicted that lithospheric removal by delamination is accompanied by magmatism, because the subsequent asthenospheric upwelling can lead to decompression melting and the generation of partial melts (e.g., Kay & Kay, 1993;Wang & Currie, 2015). At the surface, volcanic activity would be spatially (and temporally) associated with the event. However, volcanic activity may continue over a broad time period, rather than an isolated and distinct pulse, as melts are generated and extracted from the asthenosphere (cf. arc settings Comeau et al., 2015;Pritchard et al., 2018). In addition, the continental crust overlying the region of removed lithosphere undergoes basal heating, which can significantly increase the temperatures in the lower crust (e.g., Göğüş & Pysklywec, 2008;Houseman et al., 1981). This effect is observed in the numerical models, with the temperatures at depths of 40-70 km increasing by up to several hundred degrees over the course of the model evolution. Furthermore, the models show that the temperatures within the upwelling asthenosphere stay elevated for long time periods.
Widespread intraplate, mafic volcanism has occurred throughout Mongolia in the Mesozoic and Cenozoic (see Barry et al., 2003;Sheldrick et al., 2018). Geochemical analysis has revealed that its origin is sustained partial melting from a sub-lithospheric mantle source (Barry et al., 2003;Hunt et al., 2012), possibly within an asthenospheric upwelling. The volcanism lacks a pattern of age progression (e.g., Ancuta et al., 2018), but if the asthenospheric upwelling and decompression melting were extensive and occurred rapidly, then no spatial-temporal pattern would be expected (note that a Rayleigh-Taylor instability would produce a radially symmetric pattern). In addition, from geochemical studies, the temperature of the lower crust is inferred to be elevated (e.g., Ionov, 2002). Thermobarometric data predict that temperatures beneath the Hangai Dome at depths of ∼60-80 km range from ∼1100°C to ∼1300°C, and xenolith data suggest decompression (e.g., Harris et al., 2010). These temperatures are in broad agreement with the numerical models, which show temperatures of ∼1100-1250°C (i.e., >1350 K) within the asthenospheric upwelling (depths of ∼40-100 km).

Timescale for Lithospheric Removal and Renewal
Lithospheric removal and the consequences-such as surface deformation, asthenospheric upwelling, and an increase in volcanism-should occur in a timely manner, and are ultimately followed by cooling and lithospheric renewal. The simulations showed that material removal starts at a model time of ∼2 Myr and develops over the course of the simulation, with the removal of the majority of the lithospheric mantle over a small region occurring in an interval of ∼5 Myr. The duration of lithospheric removal shown in this study is generally in agreement with the range reported in other numerical studies, which vary widely from 1-15 Myr and are variously defined (e.g., Gorczyk et al., 2015;Krystopowicz & Currie, 2013;Le Pourhiet et al., 2006).
A further constraint on the timing is the fact that in central Mongolia the present-day lithospheric mantle is on the order of 20 km thick (assumed to be a small fraction of its original thickness). A portion of lithospheric mantle would be expected to be present if the lithosphere was not completely removed but only moderately thinned, or if the detachment occurred below the crust-mantle boundary, a theoretical possibility (Meissner & Mooney, 1998). Alternatively, it is more likely that the presence of a lithospheric mantle suggests renewal (and cooling) has commenced. This is interesting because it implies (significant) time has passed since the removal event, implying an age in the early Cenozoic or Mesozoic.
The renewal of a mechanically rigid lithospheric mantle by cooling of the high-temperature convecting asthenosphere occurs slowly, although there is no general consensus on the timescale. Some studies suggest cooling with a characteristic time of ∼60 Myr (Nelson, 1992; see also Turcotte, & Schubert, 2014), although it is expected to be several times higher for non-static cooling scenarios (e.g., Ueda et al., 2012). Numerical models from Ueda et al. (2012) demonstrated that renewal of 60-65 km of lithosphere can take place within 130-200 Myr of cooling, but they caution that considering mantle flow will affect the cooling rates. Other studies have shown that the asthenosphere may maintain a high temperature for longer than expected due to cratonic insulation, as well as ongoing mantle convection processes (e.g., Gorczyk et al., 2015), and therefore there may be a lengthy delay before cooling begins (possibly tens of millions of years; e.g., Ueda et al., 2012).

Constraints on the Age of a Removal Event
Recent geochemical analysis of data from Mongolia may help shed light on the possible timing of a lithospheric removal event. Geochemical analysis, including investigating high-field strength elements and isotopic signatures, from Sheldrick, Barry, Dash, et al. (2020) and Sheldrick et al. (2018) suggests the following: (1) magmatism prior to ∼107-130 Myr ago was derived from metasomatized subcontinental lithospheric mantle; (2) between ∼130 Myr ago and ∼107 Myr ago, evidence (including a gradual reduction in the involvement of garnet) indicates a significant shallowing of melt depths; (3) after ∼107 Myr ago melt was predominantly derived from an asthenospheric source, from a depth of ≥70 km with no to little variation over time; (4) after ∼50 Myr ago geochemical variations are observed and are best explained by the incorporation of metasomatically enriched material (e.g., lithospheric mantle or ancient slab material). Therefore, geochemical evidence supports the idea of a delamination event in Mongolia and provides critical timing constraints, namely that much of the lithospheric mantle material was weakened, removed, and replaced with an asthenospheric upwelling prior to ∼107 Myr ago.
A convergent regime would have existed ∼107 Myr ago from prior ocean closures; for example, the Mongol-Okhotsk ocean closure was completed 30-70 Myr previously in the latest Jurassic-earliest Cretaceous (e.g., Van der Voo et al., 2015). Previous subduction events would have likely hydrated and metasomatized the lithosphere, and may have pre-conditioned the region for a later delamination event. Crucially, within the same timeframe, thermochronology data indicate uplift and exhumation that created topography in the Hangai mountains (∼100-120 Myr ago; McDannell et al., 2018). What's more, in Mongolia, long term stability and preservation of the landscape (i.e., relic topography) is expected, on the scale of geological periods (McDannell et al., 2018).
Thus, lithospheric mantle removal by a delamination event, such as that modeled in the present study, that occurred in the Cretaceous would agree with the timing of other events in this area, most importantly the time constraints for the uplift and exhumation of the Hangai mountains determined by McDannell et al. (2018) and for the transition to an asthenospheric melt source established by Sheldrick, Barry, Dash, et al. (2020). In addition, the duration of the delamination removal event modeled in the present study does not exceed the temporal constraint of <20-30 Myr predicted from the geochemical data (Sheldrick, Barry, Dash, et al., 2020).
The numerical simulations do not preclude a removal event in the Cenozoic. However, a removal event that occurred in the Cenozoic, although potentially able to explain the Cenozoic volcanism (and any Cenozoic uplift), may not be corroborated by the available geochemical evidence. Therefore, an alternate explanation may be required to explain the Cenozoic intraplate volcanism in Mongolia. One possible mechanism, proposed to explain intraplate volcanism in other regions, is the long-lived and intermittent triggering of volcanism from bursts of paleo-slab flux (Mather et al., 2020;cf. Yang & Faccenda, 2020). This explanation is compatible with geochemical signatures from surface lavas in Mongolia that point to ancient slab material/ fluids originating in the upper mantle (Sheldrick, Barry, Dash, et al., 2020;Sheldrick et al., 2018), which may be attributed to oceanic slabs ponded at the mantle transition zone (e.g., from the Paleo-Asian, Mongol-Okhotsk, and Paleo-Pacific plates; Sheldrick, Hahn et al., 2020).

Conclusions
In this study we investigate aspects of lithospheric removal, including the conditions under which it can occur and its effect on surface deformation, such as the amount and timing of topographic uplift. To do this we generated realistic thermo-mechanical numerical simulations, using the geodynamic code ASPECT, with the choice of modeling input parameters guided by geological, geochemical, and geophysical data. We find that lithospheric removal by delamination develops dynamically and self-consistently if a dense and weak lower crust is present. We generate a dense crustal root in the simulations by implementing a phase transition and a density jump to reproduce the effect of metamorphic eclogitization. In addition, a moderate to high temperature at the crust-mantle boundary and a moderate convergence rate are also required. A pre-existing weak zone, interpreted to be a local shear zone, a fault, or a post-collisional suture zone, is included in the models and is found to promote and facilitate the delamination by focusing deformation and creating an initiation point for detachment.
We apply the results of these simulations to central Mongolia, which is a prime example of intraplate surface uplift because of its location far into the continental interior and its enigmatic high topography. We find that there is evidence that the conditions and requirements for delamination exist in central Mongolia-including a weak lower crust. Furthermore, the main signature of a delamination event, an anomalously thin or removed lithosphere, has been detected beneath central Mongolia.
The outputs of the simulations generally match the observational evidence in Mongolia. Notably, a broad dome-shaped pattern of topography and an elevated surface are produced. Most significantly, variations in the thickness of the lithosphere (i.e., doming of the lithosphere-asthenosphere boundary) and the size of an asthenospheric upwelling are consistent with that inferred from electromagnetic and seismic imaging, and supported by geochemical results. A mantle upwelling leads to an elevated temperature at the crust-mantle boundary and (shallow) mantle decompression melting, both of which are compatible with available petrological evidence, and links a delamination event to an increase in magmatic and volcanic activity. It is hypothesized that previous ocean closure and subduction events in Mongolia likely hydrated and metasomatized the lithosphere, and may have pre-conditioned the region for a delamination event. The results of this study suggest that removal of the lithospheric mantle by a delamination event is a physically plausible explanation for the observed surface uplift, as well as the volcanism, in central Mongolia.

Data Availability Statement
The geodynamic modeling software, ASPECT v2.1.0, is available from Bangerth et al. (2019), https://doi. org/10.5281/zenodo.2653531, and the input parameters for numerical models and tests can be found in Table 1 and Figure 10.

Acknowledgments
Additional explanatory details can be obtained in the Supporting Information, in addition to a video of the simulation. Computational resources were provided by the Institute for Geophysics at the University of Münster. This research was financially supported by DFG (Deutsche Forschungsgemeinschaft; German Research Foundation) grant BE5149/6-1 awarded through the DACH program. MB is funded through a DFG Heisenberg Grant (5149/7-1, 5149/9-1). We thank Alexey Kuvshinov for discussions relating to this study. We thank the reviewers for providing very detailed, thoughtful, and constructive comments on the manuscript. We thank the Computational Infrastructure for Geodynamics (geodynamics.org) which is funded by the National Science Foundation under award EAR-0949446 and EAR-1550901 for supporting the development of ASPECT. Open access funding enabled and organized by Projekt DEAL.