Reproducing past subduction and mantle flow using high‐resolution global convection models

Plate subduction drives both the internal convection and the surface geology of the solid Earth. Despite the rapid increase of computational power, it remains challenging for geodynamic models to reproduce the history of Earth‐like subduction and associated mantle flow. Here, based on an adaptive approach of sequential data assimilation, we present a high‐resolution global model since the mid‐Mesozoic. This model incorporates the thermal structure and surface kinematics of tectonic plates based on a recent plate reconstruction to reproduce the observed subduction configuration and Earth‐like convection. Introduction of temperature‐ and composition‐dependent rheology allows for incorporation of many natural complexities, such as initiation of subduction zones, reversal of subduction polarity, and detailed plate‐boundary dynamics. The resultant present‐day slab geometry well matches Benioff zones and seismic tomography at depths < 1500 km, making it possible to hindcast past subduction dynamics and mantle flow. For example, the model produces a flat Farallon slab beneath North America during the Late Cretaceous to Early Cenozoic, a feature that has been geodynamically challenging to reproduce. This high‐resolution model can also capture details of the 4‐D evolution of slabs and the ambient mantle, such as temporally and spatially varying mantle flow associated with evolving slab geometry and buoyancy flux, as well as the formation of shallow slab tears due to subduction of young seafloors and the resulting complex mantle deformation. Such a geodynamic framework serves to further constrain uncertain plate reconstruction in the geological past, and to better understand the origin of enigmatic mantle seismic features.


Introduction
Plate tectonics is a unique characteristic of the solid Earth. It provides a framework to the study of the coupled crust-mantle system. Since the establishment of the theory of plate tectonics in the 1960s, much progress has been made in better understanding the Earth's tectonic history. A plate reconstruction back to 180 Ma ago proposed in the 1980s (Engebretson et al., 1985) provided a quantitative analysis for the past motion of tectonic plates. This type of reconstruction has been evolving with increasing amounts of detail, as observational constraints on palaeomagnetics and reference frames have accumulated . Ultimately these reconstructions provide an important linkage between surface plates and the underlying mantle. Inside the Earth, seismic tomographic inversions have, with increasing resolving power, outlined fast anomalies beneath the subduction zones (e.g. Grand, 2002;Ritsema et al., 2011;Obayashi et al., 2013), most of which could represent subducted slabs. These slabs, on the one hand, are thought to drive large-scale mantle convection, and, on the other hand, could pull their surface plates into the mantle, both of which may help to drive the tectonic plates (Conrad and Lithgow-Bertelloni, 2002;Stadler et al., 2010;Becker and Faccenna, 2011).
However, the subduction history and resulting mantle evolution inferred from plate reconstruction and that from seismic tomography do not always agree (van der Meer et al., 2010;Wu J et al., 2016). This reflects the uncertainties in these two different approaches: the former suffers from inadequate information on subducted seafloors; the latter from limited seismic resolution and non-unique interpretation of mantle structures. In addition, true polar wander may depend on modulation of the internal mass distribution of Earth by subduction and convection (Steinberger and Torsvik, 2010), complicating interpretation of palaeomagnetic data and the reconstruction of plate motions. Finally, our understanding of the geometry and dynamics of subducting slabs remains limited, as does the resulting tectonics processes, such as trench migration (e.g. Schellart, 2008), orogenic processes (De-Celles et al., 2009), and the temporal and spatial distribution of volcanisms (e.g. Trumbull et al., 2006). Geodynamic modeling with data assimilation represents a promising way to reconcile the different inferences on past subduction and to further understand the relevant tectonic processes (Liu LJ and Stegman, 2011;Hu JS et al., 2016;Zhou Q et al., 2018). Global models, without artificial vertical walls as those in regional mod-els, have obvious advantages in resolving large-scale, wholemantle dynamics. However, limited computational power has been a rate-limiting factor. Until now, published global time-dependent models (McNamara and Zhong SJ, 2005;Foley and Becker, 2009;Crameri et al., 2012;Coltice et al., 2012;Bello et al., 2015;Bower et al., 2015) have fallen largely into two categories. The first category adopts a free-slip upper surface with self-merging plate motion (e.g. Crameri et al., 2012;Coltice et al., 2012). Such models are often used to investigate the physical properties that control the plate behavior (Mallard et al., 2016;Coltice et al., 2012). The other category of models assimilates plate motions on the surface that guide subduction over time (e.g. McNamara and Zhong SJ, 2005;Bello et al., 2015;Bower et al., 2015). These models, although not fully dynamic, are geographically oriented and better for studying Earth processes, such as the evolution of LLSVPs (Mc-Namara and Zhong SJ, 2005), the motion of hotspots (Hassan et al., 2016), and topographic evolution (Flament et al., 2014). However, none of these models (including both categories) could simultaneously generate all subduction characteristics, including single-sided subduction, temporally varying slab dip angle, complex slab deformation and tearing, robust subduction-induced mantle flow, and observed present-day mantle structures.
In this paper we present a high-resolution global mantle convection model since the mid-Mesozoic using an adaptive approach of sequential data assimilation. In this model we incorporate a relatively realistic rheology, including temperature-and depth-dependent viscosity and composition-based weak oceanic crust, ensuring the generation of Earth-like subduction. The model also assimilates seafloor age and plate motion history from a recent plate reconstruction . Relative to published works of its kind, our model could better represent the evolution and morphology of slabs, especially at upper mantle depths, allowing us to investigate, in greater detail, local-scale subduction dynamics, the temporal evolution of mantle deformation, and the origin and implication of upper mantle seismic structures.

Method
We simulate global subduction since 150 Ma using the sequential data-assimilation technique (Liu LJ and Stegman, 2011;Hu JS et al., 2016). Compared to our previous work, this new model is more advanced in that it now assimilates a thermal-chemical lithosphere structure and allows more realistic trench configuration and subduction dynamics, including formation of new subduction zones and reversal of subduction polarity. The simulation was carried out using the 3-D spherical finite element code CitcomS (Tan E et al., 2006;Zhong SJ et al., 2008). We assume that the mantle is incompressible and satisfies the Boussinesq approximation. The model covers the whole mantle domain, which is discretized into 12 caps, and each cap has 257×257×97 nodes. The model has a uniform lateral resolution, which is ~23 km. The vertical resolution varies with depth, and is finer as it approaches the surface (~12 km) and the core-mantle boundary (~20 km). The average resolution of this model is higher than that in earlier similar studies. For example, the total number of grids is about 6 times higher than in two recent global models Bello et al., 2015). The model properly resolves the thermal-mechanical structure of the lithosphere, such as the weak oceanic crust, which is critical in generating single-sided subduction (Crameri et al., 2012).

Assimilating Plate Motion and Thermal Lithosphere
The model takes plate motion from the recent plate reconstruction of Müller et al. (2016) as the surface boundary condition ( Figure 1a). These data are processed using the open source paleo-geographic software GPlates (www.gplates.org/; Gurnis et al., 2012), and output at every Ma. At each model time step that is much smaller than 1 Ma, the surface velocity condition is interpolated between two adjacent input conditions. The assimilation of plate motion ensures the slab subducts at the proper geographic location, which is important for matching the observed presentday slab geometry (Liu LJ and Stegman, 2011) and the associated mantle flow (Hu JS et al., 2017).
We also assimilate and update the thermal structure of the lithosphere at each time step according to the reconstructed seafloor ages (Figure 2a). In an oceanic lithosphere, we use a modified halfspace cooling model to define its thermal structure: where T l is lithosphere temperature, T m mantle temperature, z depth, t seafloor age, κ thermal diffusivity, and A modification factor as a function of t. As in most global models, the limited numerical resolution restricts its capability of reproducing Earth convection vigor while maintaining a thin thermal boundary layer and subducting slabs (e.g., Davies et al., 2007). Therefore, we use a Rayleigh number of 3.12×10 8 (with the length scale being Earth's radius) and a non-dimensional mantle temperature of 0.58 that corresponds to 700 °C, about half of that of the Earth. In order to maintain the buoyancy flux of the subducting oceanic plate and thus the correct convection vigor, we introduce an age-dependent modification factor A, that decreases the temperature gradient of the lithosphere, so that the depth integral of lithospheric thermal buoyancy equals that of the exact half-space cooling profile with a mantle temperature of 1300 °C. This way, we make sure that the buoyancy of the subducting plates remains similar to that in Earth, which is important for the study of mantle flow and associated surface expressions. Within the continental lithosphere we assume an initial conductive thermal profile whose temperature increases from 0 on the surface to 0.5 at a depth of ~160 km ( Figure 2a). This way, interiors of the continental lithosphere eventually grow thicker than most oceanic plates, maintaining a quasisteady state thermal structure (e.g., Figure 5) similar to that of the Earth (McKenzie et al., 2005).
To guarantee that a smooth thermal profile translates from within the oceanic plates into the down-going slab, and to make sure that other plate boundary properties such as composition and viscosity could evolve self-consistently, we avoid updating the temperature of the subducting and overriding plates within 600 km horizontal distance from the trench. As discussed later, we also avoid updating the corresponding composition in the vicinity of the trench, such that tracer-carried properties (density and viscosity) will behavior naturally over time. We name this an adaptive data assimilation approach.

Assimilating Chemical Composition
We assume chemical layering both in the oceanic and the continental lithosphere ( Figure 3). The buoyancy and viscosity associated with different compositions significantly affect the dynamics of the system (e.g. Crameri et al., 2012). In our model, the compositional density anomaly is defined using the ratio method: where B is the buoyancy ratio (Table 1), an input parameter that varies with different compositions, Δρ ch the compositional density anomaly, ρ 0 the reference density, α 0 the reference thermal expansivity, and ΔT the temperature contrast across the mantle. The compositional effect on viscosity is achieved through a multiplication factor that is a geometric mean for all the compositions within the element: where C is a prefactor, C i is the viscosity multiplier for composition i (Table 1), and r i is the fraction of composition i over the total elemental composition.
Within the oceanic plates we define two chemical layers distinct from the ambient mantle: a weak crustal layer and a buoyant layer (Figure 1c, f; Figure 3). The former represents the uppermost portion of the oceanic plate and has a thickness of ~25 km. This layer is neutrally buoyant but has a low viscosity at 10 19 Pa·s, thus acting as a lubricating layer to decouple the subducting plate from the overriding plate ( Figure 1e; Crameri et al., 2012). However, since we impose the surface velocities, we turn off this viscosity effect beyond 200 km from the subduction zone on the side of the subducting plate. We also assimilate a buoyant layer within the oceanic plates at a depth of 40-60 km. This layer has a buoyancy ratio of -0.83. The total buoyancy of this layer is similar to that of an oceanic crustal layer with a thickness of 7 km and a density of 3.0 g/cm 3 . Effectively, this layer carries the buoyancy of the oceanic crust, but avoids getting accreted to the upper plate during subduction since it is within the relatively strong slab.
We approximate the continental plate as two compositional layers, the crust and the mantle lithosphere ( Figure 1f; Figure 3). The crust extends from 0 to 38 km depth and has a buoyancy ratio of -4.4, while the mantle lithosphere extends from 38 km to 160 km depth and has a buoyancy ratio of -0.42. The buoyancy ratio of the crust is consistent with a mean density of 2.8 g/cm 3 , while that of the mantle lithosphere largely cancels out that of the thermal  ity of 10 19 Pa·s, mimicking those of the arcs in real Earth. It helps to decouple the two plates on both sides of the subduction zone during subduction. In addition, to allow the composition of the mantle wedge to evolve self-consistently, we do not assimilate the composition of the overriding plate < 600 km from the trench.
At the initial time step we also define a uniform thermal-chemical layer above the core-mantle boundary (CMB), which represents the source of the LLSVPs beneath Africa and Pacific. This layer initially has a thickness of 250 km, similar to that in Zhang N and Li ZX (2018), but maybe slightly thicker than estimates from tomography (e.g. Burke et al., 2008); it has a buoyancy ratio of 0.67, equivalent to a density anomaly of 2.4%, similar to that in Hassan et al. (2016); and it has a viscosity multiplier of 0.1. The temperature is 500 °C higher than the ambient mantle, which is not as high as in earlier studies such as Hassan et al. (2016) and Zhang N and Li ZX (2018), due to the relatively low Rayleigh number adopted in our model. It eventually evolves into two thermal-chemical piles at the present, similar to findings in previous studies (e.g. Mc-Namara and Zhong SJ, 2005). However, we emphasize that the structure of LLSVPs is not the focus of this study. Future studies will elaborate the properties of this basal layer and will better constrain its structure and evolution.

Viscosity Structure
The viscosity of this model (Figure 1b, e) is depth-, temperatureand composition-dependent: , where η is the non-dimensional viscosity, η 0 is the background viscosity, C is viscosity multiplier due to composition as defined in equation (3), E η is the activation energy, T η is the temperature offset, T m is the non-dimensional mantle temperature (Table 1), and T is the non-dimensional temperature. The activation energy is depth-dependent, ranging from 17 kJ·mol -1 in the lower mantle to 42 kJ·mol -1 in the upper mantle, whose values are lower than laboratory estimates (e.g. Karato and Wu P, 1993). In such a global model, however, due to the lack of extra weakening mechanisms such as plastic yielding and dislocation creep, as well as faulting and hydration, we need to use a lower activation energy to reduce overall strength of slabs during subduction.
The radial variation of viscosity for the initial condition is shown beneath both continental and oceanic plates ( Figure 2b). In these viscosity profiles, we use a strong lithosphere whose viscosity reaches 10 23 Pa·s, and a relatively weak asthenosphere whose viscosity is 10 20 Pa·s. The viscosity increases from the asthenosphere to the transition zone, and then to the middle lower mantle with a maximum of 10 23 Pa·s. Then it decreases to 2×10 20 Pa·s as it approaches the CMB. This radial viscosity profile is similar to those constrained by geophysical inversions (e.g. Steinberger and Calderwood, 2006) and by reproducing present slab geometry (e.g. Hu JS et al., 2016). The increase of viscosity at ~1000 km depth is in accordance with recent inferences based on geoid (Rudolph et al., 2015), seismic tomography , and petrological experiments (Marquardt and Miyagi, 2015). The initial thermal structure of the lithosphere has a step at 120 km depth ( Figure 2a), but the step is quickly smoothed out after a few million years due to thermal diffusion. Note that we reduce the viscosity of the slab hinge ( Figure 3b) by imposing a maximum viscosity cut-off at 2.5×10 22 Pa·s, to mimic the effect of plastic yielding (Billen and Hirth, 2007).

The Effect of Rheological Features on Subduction
The model has significant lateral viscosity variations (Figure 1b, e), which reach 4 orders of magnitude, ranging from 10 19 to 10 23 Pa·s. This lateral viscosity variation is critical in generating Earth-like plate subduction. For example, the absence of the weak oceanic crust causes a thickened slab that tends to get congested at the subduction zone (Figure 4a, b), while the incorporation of the weak oceanic crust can effectively reduce the coupling between the subducting and overriding plates, resulting in an asymmetric, slim slab entering the trench (Figure 4c, d).
Compared with the weak oceanic crust, the arc-like weak zone and weak slab hinge have less significant impact on the slab evolution, but both of them help in generating Earth-like subduction. For example, without the arc-like weak zone the subduction is congested for the first several million years, as can be seen from the thickened slab at 145 Ma ( Figure S1a, b). However, this congested slab can subsequently get through, due to the imposed plate velocity boundary at the surface. Once the subduction is fully established, the weak oceanic crust would act as a lubricating layer, preventing slab congestion ( Figure S1a, b). On the other hand, a strong slab hinge tends to reduce slab bending, resulting in a flatter slab ( Figure S1c, d) and preventing the slab from retreating. This rigid system significantly slows down the convergence rate and causes numerical difficulties. For the following sections, we focus on the standard model that implements all three of these rheological features.

Results
As stated earlier, published global models using data assimilation do not always produce realistic upper mantle processes. In this section, we present several examples of how the complex subduction dynamics are realized in our model, including formation of a new subduction system and reversal of subduction polarity. Then we will show the predicted present-day slab structure and further compare model predictions with actual seismic observations. Importantly, our model reproduces the past evolution of subduction and mantle flow, which are useful for better connecting tectonic/geologic phenomena with relevant deep mantle processes. Examples include the formation of abnormal subduction such as the Farallon flat slab and the evolving 3-D mantle flow.

Simulating Complex Subduction Dynamics
The Earth's subduction history is very complex. For example, a number of trenches have formed during the past 150 Ma, along with several events of subduction polarity reversal . In addition, some trenches migrate very fast, such as the South Sandwich Trench in the east of the Scotia Sea (Barker, 2001). Sometimes trenches interact with each other, resulting in a complex situation such as that observed surrounding the Philippine Sea Plate (Wu J et al., 2016). Modeling the appearance, disappearance, and changing morphology of trenches has been challenging for geodynamic models, even when surface kinematics are imposed, such as that done in regional models (Liu LJ and Stegman, 2011;Hu JS et al., 2016). Modeling Earth-like subduction with complex trench kinematics requires both high numerical resolution and sophisticated rheology; neither is practical in global simulations. To overcome such difficulties, some earlier models adopted a "slab assimilation" approach in which the morphology, thermal structure, and sinking velocity of slabs are imposed down to upper mantle depths . Obviously, these kinematic prescriptions limit the models' capability to capture the dynamics associated with slab evolution.
Here, we simulate subduction in a more dynamic fashion, by us-  ing a weak oceanic crust, a weak slab hinge, and self-evolving plate boundary properties to reproduce Earth-like asymmetric subduction under complex trench kinematics (Figure 3). Take the formation of a new subducting slab as an example: when a new trench appears, the composition of the emerging mantle wedge becomes less viscous and more buoyant, mimicking the property of a volcanic arc, as mentioned in Method (Figure 3a). This arc-like structure increases lateral density variation within the subduction zone, which helps to initiate slab subduction. In addition, the low viscosity of this structure facilitates convection and return flow inside the mantle wedge. Meanwhile the weak oceanic crust effectively decouples the slab from the overriding plate, and the weakened slab hinge further facilitates slab bending (Figure 3b). Consequently a new slab develops in respond to the convergent surface boundary condition, without formation of congested slab at shallow depth (Bello et al., 2015). As subduction proceeds, the weak mantle wedge initially accommodates convergence through internal deformation (Figure 4c, d), followed by establishment of a well-behaving slab lubricated by the weak oceanic crust from the overriding plate (Figure 4c, d). Asymmetric subduction will continue as long as the two plates keep converging. It is also important to realize that, during this simulation, the adaptive data assimilation scheme avoids overwriting properties near the convergent boundaries and thus allows the key subduction features to evolve self-consistently over time. Figure 5 illustrates how this process works during the formation of an intra-oceanic subduction zone, where the Farallon plate subducted beneath the proto-Caribbean plate . At 85 Ma, the subduction had just started (Figure 5a), so there was not yet a subducting slab (Figure 5b). Initially the mantle wedge was weakened due to the emerging arc-like structure and there was a weak spot on the surface of the plate right on the trench representing the weak oceanic crust (Figure 5b). 5 Ma later a clear asymmetric thermal structure beneath the trench developed, from the upwelling within the mantle wedge and downwelling of the slab along the trench (Figure 5d). At 70 Ma the subducting slab was fully developed and had reached ~350 km depth (Figure 5f).
In this model, we could reproduce mantle dynamics associated not only with subduction initiation but also with trench polarity reversal. For example, part of the southwestern Pacific trench shifted its polarity at 85 Ma ( Figure 6 Müller et al., 2016). Before 85 Ma, the Pacific plate subducted obliquely beneath the Australian Plate, resulting in a segmented westward-dipping slab (Figure 6a, b). After the polarity reversal occurred, the slab started to subduct eastward (Figure 6c, d). At 70 Ma an eastward-dipping slab was formed (Figure 6e, f). Significantly, because we allow the temperature and composition of the overriding plate near the trench to deform dynamically, the resulting evolution of these features in the mantle wedge look very natural (Figs. 6f, 5h).
Our adaptive data assimilation scheme could also reproduce complex subduction systems where multiple trench segments intersect with each other or where fast trench migration occurs. One such example is the double subduction system surrounding the Philippine Sea Plate (Figure 7). The Izu-Bonin-Mariana trench intersected with the Ryuku-Nankai trench and migrated northeast-ward rapidly between ~50 Ma and ~30 Ma . In a cross section that transects the Philippine Sea and South China, we found that the double slab system beneath the two trenches was well resolved (Figure 7b). Even the prominent slab beneath China that resulted from the subduction of the former Izanagi Plate was reproduced. Therefore we demonstrate that the fine model resolution and sophisticated compositional functions of our adaptive data assimilation approach, likely for the first time, is able to simulate complex subduction dynamics in a global scale.

Reproducing the Present-Day Mantle Structure
The present-day mantle seismic structure ( Figure 8) is one of the most important constraints on the validity of all realistic subduction models. On the global scale, the predicted present thermal structure (Figure 8) shows that some of the oldest (150 Ma) slabs have already piled up above the CMB, consistent with earlier studies of the characteristic time scale of mantle convection (Bunge et al., 1998). In addition, the boundary between the subducted slab and the thermal-chemical piles ( Figure 8) largely resembles the observed geometry of the LLSVPs beneath Africa and Pacific revealed by tomography models (e.g. Grand, 2002;Ritsema et al., 2011), implying that the predicted lowermost mantle structure is also consistent with observation.
At regional to local scales, the predicted slab structures are evaluated against observations that include Benioff zones (Hayes et al., 2012) and seismic tomography (e.g. Ritsema et al., 2011;Obayashi et al., 2013). We first consider a recent P-wave model GAP_P4 (Figure 9; Obayashi et al., 2013) and a S-wave model S40RTS ( Figure S2; Ritsema et al., 2011) in map view from 250 to 1500 km depth. In the upper mantle (above 660 km), the predicted slabs almost always coincide with seismically fast regions (Figure 9a, b, c; Figure S2a, b, c). In the lower mantle, the predicted slabs are also generally consistent with fast anomalies in these tomography models (Figure 9d, e, f; Figure S2d, e, f). For example, the predicted Nazca slab fits the tomography down to 1000 km beneath South America for both the P-wave and S-wave models. The Sumatra slab, too, fits the tomography down to 1000 km and its disappearance at 1500 km is also consistent with tomography models.
To better examine the fit between the predicted and observed slab structure, we further plot several cross sections along different locations (Figure 10). The first is an east-west cross section through the North American plate, under which the ancient Farallon slab subducts (Figure 10a). The predicted present Farallon slab is located under the eastern part of North America, with an eastdipping orientation. This slab is now disconnected from the trench, due to the subduction of the Pacific-Farallon mid-ocean ridge since ~30 Ma and the cessation of convergence afterwards . Both the location and geometry of the Farallon slab are consistent with those outlined by tomography, including GAP_P4 , S40RTS (Ritsema et al., 2011) and LLNL_G3Dv3 (Simmons et al., 2012).
In another cross section that runs across eastern Asia (Figure 10b), the predicted slabs are also largely consistent with fast anomalies in tomography. For example, the subducted Izanagi slab piles up above the CMB, resembling the fast anomalies in the lower mantle, especially for the P-wave tomography GAP_P4 and IINL_G3Dv3 (Figure 10b). The dip angle of the Pacific slab at the Izu trench is slightly shallower than the interpolated Benioff zone (Hayes et al., 2012). This is mainly because the slab is too close to the trench-trench-trench triple junction on the north, which is dynamically challenging to reproduce. The last cross section shows the slab structure across the Sumatra-Java trench (Figure 10c). The predicted slab extends to ~1000 km depth, which is consistent with tomography models. The interpolated Benioff zone (Hayes et al., 2012) reveals a normally dipping slab that also fits our prediction quite well. To avoid redundancy, we show only three cross sections here. Nevertheless, we emphasize that the predicted slabs in our model also match well the fast anomalies beneath other major subduction zones, especially in the upper mantle, such as those in South America, Tonga, and the Philippines (Figure 9 and Figure S2). However, the cross sections in both eastern Asia ( Figure 10b) and Sumatra (Figure 10c) show some mismatches in the transition zone or at about 1000 km depth. Compared with tomography, the predicted slab has less volume at these depths, and the slab has a dipping geometry in contrast to the flat-lying fast anomalies shown in tomography. Traditionally these fast anomalies have long been interpreted as stagnant slabs (Fukao et al., 1992(Fukao et al., , 2009Huang JL and Zhao DP, 2006;Wu J et al., 2016). Although a stronger mechanical barrier at 660-km or 1000-km depth might result in more slabs accumulating inside the transition zone or above 1000 km depth, this does not seem to be the case because slabs at other parts of the world in the current model match observations quite well. Alternatively, these transition zone fast anomalies may represent structures not associated with oceanic slabs, and one such example is delaminated lithospheric roots as suggested below the southern Atlantic (Hu JS et al., 2018). A more detailed understanding of these features will require a significant amount of future work, beyond the scope of the current paper.

Temporal Evolution of Subduction and Mantle Flow
Relative to earlier global studies, our model has the advantage of better reproducing upper-mantle subduction and mantle flow, both of which are crucial for understanding the Earth's tectonic evolution. Here we use the formation of the Farallon flat-slab and the history of mantle flow as examples.
The particular tectonic history of the western United States, including the Laramide orogeny (e.g. Saleeby, 2003;DeCelles, 2004), inboard migration of crustal deformation and arc magmatism (e.g. Snyder et al., 1976;Coney and Reynolds, 1977) and the formation of the Western Interior Seaway (e.g. Cross and Pilger, 1978;De-Celles, 2004), suggests that the Farallon slab was flatly subducting during the Late Cretaceous to early Cenozoic (Bird, 1988 (O'Driscoll et al., 2009;Liu SB and Currie, 2016) or assumed a known geometry of flat slab (Bird, 1988;English et al., 2003). Nor have recent global models with data assimilation properly reproduced the evolution of this subduction process Bello et al., 2015).
Here we show that a Farallon flat slab can be reproduced in our model (Figure 11). At ~80 Ma, the Farallon slab developed a shallower-than-normal dip angle (Figure 11a). Since then, the slab dip  angle progressively decreased and reached a maximum flatness at ~50 Ma (Figure 11b), when the flat slab extended more than 1000 km inland from the trench. Subsequently, the flat slab began to detach from the upper plate. At the present day this slab is located mostly beneath the eastern US ( Figure 11c); its predicted configuration matches tomography quite well (Figure 9, 10a). According to our analysis the slab anchoring in the lower mantle, the fast west-ward motion of the overriding plate, and the progressive younging of the Farallon plate since the Cretaceous all contribute to the flattening of the slab. However, the predicted slab in the mid-lower mantle is slightly too steep compared with tomography, which seems to suggest that other mechanisms, such as hydrodynamic suction of the thick Colorado Plateau (O'Driscoll et al., 2009;Liu SB and Currie, 2016) or the subduction of Shatsky and Hess conjugates (Liu LJ et al., 2010), may also have played a role.
Mantle flow represents another key aspect of Earth evolution, since it could drive plate motion (e.g. Becker and Faccenna, 2011), transport chemical compositions (e.g. Hassan et al., 2016), and produce seismic anisotropy (e.g. Hu JS et al., 2017). Our high-resolution global subduction model has produced a complex, evolving mantle flow field that may have many important implications. At the plate scale, the upper-mantle flow pattern changes both with space and with time ( Figure 12). Spatially, the flow beneath a subducting plate is mostly influenced by the surface plate motion through Couette flow, such as that below Pacific, Australia, India, and Africa. In contrast, mantle flow below an overriding plate sees a stronger influence from the nearby sinking slab that forms the Poiseuille flow above, as is seen in the Americas and Eurasia. These observations are consistent with our recent conclusions based on a regional subduction model (Hu JS et al., 2017). Similarly, the temporal variation of upper-mantle flow is also affected both by changes in surface motion and by the evolution of slabs. For example, the flow history beneath the Pacific mostly follows the plate reorganization at the surface (Figure 12a vs. Figure 12b), as suggested by Seton et al. (2015), while that below North America and South America is more subject to the geo-  (Figure 12a vs. Figure 12b). It is also clear that predicted mantle flow history should have important implications for the resulting mantle and lithosphere deformation and the interpretation of seismic anisotropy. At the regional scale, our high-resolution model could also capture the complex mantle deformation. An example is the northwestern South America, where the model predicts a slab tear. This generates fast mantle flow connecting the Pacific upper mantle with the Atlantic mantle ( Figure 13). Earlier studies relate the change of SKS orientations along the arc to a postulated slab tear and the resulting mantle flow (e.g. Porritt et al., 2014;Idárraga-García et al., 2016). Our model provides a confirmation of this proposed slab tear and change in mantle flow. Mechanically, this tear formed due to the overpressure beneath the mechanically weak young slab, similar to what we find in the Juan de Fuca slab . Our predicted slab geometry is also consistent with the recent tomographic image that reveals a semi-vertical wall of fast anomalies below the region (Levander et al., 2014), providing an additional validation to the associated mantle flow.

Discussion
Modeling the Earth's past subduction in a fully dynamic way, i.e., without imposed lithosphere structures or surface kinematics, is still impractical at the moment. Currently, data-oriented geody-namic modeling remains a promising means to the study of Earth's past tectonics (e.g. Liu LJ and Stegman, 2011;Bower et al., 2015;Bello et al., 2015). Although generally successful in solving many questions, recent global studies of this kind pose some limitations in simulating subduction dynamics. For example, Bower et al. (2015) prescribed thermal structure and sinking velocity of slabs in the upper mantle, which restricts the self-evolving of the slabs. Bello et al. (2015) used a low Rayleigh number and scaled down the surface velocities, which results in a sluggish mantle convection and strong coupling between subducting and overriding plates. Along this direction, the high-resolution global convection model using an adaptive data assimilation scheme presented in this paper represents one of the most advanced accomplishments. On the one hand, it can simulate self-evolving asymmetric slabs and sophisticated historical geodynamic processes such as formation of new subduction zones and reversal of subduction polarity, as well as complex subduction systems with multiple slabs near each other (Figures 5, 6, 7). On the other hand, the resultant present-day mantle structure matches well with tomography and Benioff zones ( Figures. 9, 10) and the associated mantle flow history provide new insights for understanding mantle and lithosphere deformation ( Figure 12).
We emphasize that, compared to recent studies of its kind, our model could better represent the slab structure and dynamics, especially those along subduction zones and at shallow mantle  , S40RTS (Ritsema et al., 2011) and LLNL_G3Dv3 (Simmons et al., 2012). (a) Cross section that goes through the U.S.. (b) a cross section through East Asia. (c) Cross section through the Sumatra-Java trench. The upper panels show the temperature field from our model, while those lower panels show the velocity anomaly from the tomography models. All panels are overlain with black lines that represent the nondimensional temperature contours with a dimensional interval of 200 K. The thick purple lines are the interpolated Benioff zones (Hayes et al., 2012). There are no Benioff zones beneath the eastern U.S.. Note that, due to the selected projection, the distortion in lateral distance increases with increasing depth.

200
Earth and Planetary Physics doi: 10.26464/epp2018019 depths. These model properties are important, because the upper mantle represents the most dynamically active portion of the solid Earth. For example, instantaneous geodynamic models usually rely on the slab geometry inferred from interpolated Benioff zones (e.g., Stadler et al., 2010;Eakin et al., 2014), but our recent study in South America suggests that such a slab configuration is problematic because these flat slabs are internally torn and are not well represented by the sparse intra-slab earthquakes (Hu JS and Liu LJ, 2016). This global model, expanded from the regional model in Hu JS et al. (2016), provides a framework for systematically evaluating the geometry and dynamics of subducting slabs on Earth.
While consistency at different spatial scales between our predicted and observed mantle structures validates the past evolution of subduction, there are some local discrepancies, which we interpret as reflecting, at least partly, incorrect plate motions in the region. For example, the slabs beneath the Indian plate do not fit P-wave tomography at 1500 km depth (Figure 9f). The fit for Swave tomography is better, although still not perfect ( Figure S2f). In our model these slabs resulted from the intra-oceanic subduction within the Neo-Tethys, part of a two-stage India-Eurasia collision system (van de Voo et al., 1999). This intra-oceanic subduction likely terminated at ~50 Ma , in agreement with our model. In fact, the slabs in this region fit the deep- er fast anomalies better, which seems to suggest that the intraoceanic subduction should have stopped slightly earlier or the sinking rate of the slab should have been slightly faster. However, since most other lower mantle slabs match tomography fairly well, we suggest the misfit below the Indian plate implies that plate motion history for this region is incorrect. Another obvious mismatch between predicted and observed mantle structure is the widespread transition zone fast anomalies below East and Southeast Asia, most of which are not reproduced by the known subduction history of the region. As discussed earlier, we suggest these fast anomalies may represent delaminated lithospheric roots (Hu JS et al., 2018), although detailed studies on their physical or chemical properties are still required to verify this hypothesis.
Since the model properly captures the evolution of subduction, it also has potential to better reproduce the driving mechanisms for geological events related to past subduction. Examples include the formation of the Farallon flat slab ( Figure 11) and the past upper-mantle flow ( Figure 12). Consequently, these allow us to further understand the deformation of the mantle on a global scale and the evolution of continental lithosphere. Immediate future research is required to connect the predicted mantle flow with seismic anisotropy, as that recently done in South America (Hu JS et al., 2017). Similarly, further implementation of realistic overriding plate structure and rheology as those done in regional models (Leng W and Gurnis, 2011) could help further improve our understanding of continental lithosphere dynamics.   vigor associated with the thermal-chemical piles tends to be underestimated. This may explain why this model does not produce as many isolated plumes as in Hassan et al. (2016) and Zhang N and Li ZX (2018). To better represent the thermal boundary layer at CMB, a higher temperature is probably required, as suggested by the high CMB heat flux (5-15 TW) estimated by core temperature, geodynamo energetics, plume flux, and phase change in magnesium silicate perovskite (Lay et al., 2008). However, a higher Rayleigh number would increase numerical difficulties. Due to limited computational power, global models with a higher Rayleigh number have to sacrifice detailed subduction zone structures and dynamics such as in Bower et al. (2015). Therefore, we await either more powerful supercomputers or better data assimilation schemes to solve this problem. Another limitation is related to the Boussinesq approximation applied in our model. Due to the omission of shear heating and internal heating, the lower mantle tends to be cooled gradually by the incoming slabs. A model with extended Boussinesq approximation and sufficient basal and internal heating (Hassan et al., 2016), can likely solve this problem. However, we emphasize these limitations would not affect conclusions regarding shallow mantle structures, which are the focus of this study.

Conclusion
We present a high-resolution geodynamic model that simulates global subduction since 150 Ma using adaptive data assimilation. Figure 14 summarizes the predicted present-day mantle structure and flow. The model's assimilated surface kinematics and lithosphere thermal profile, and its composition-dependent rheologic-al features, such as the weak oceanic crust, arc-like weak zones and weak slab hinges, facilitate realization of asymmetric subduction and dynamically evolving down-going slabs. This improved model capability can better deal with many other natural complexities, such as the initiation of subduction zones, reversal of subduction polarity, and detailed plate-boundary dynamics. More importantly, the model produces a global mantle slab structure that well fits the Benioff zones and seismic tomography (especially above 1500 km depth). Consequently, the model captures the temporal and spatial variation of slab dip angle, such as that associated with the formation and cessation of the Farallon flat slab. It also predicts the temporally and spatially varying mantle flow associated with evolving slab geometry and buoyancy flux. Due to the nature of data assimilation, i.e. geographically oriented, this model can be used to further study many geological events related to subduction processes.