Assimilating lithosphere and slab history in 4-D Earth models

We develop methods to incorporate paleogeographical constraints into numerical models of mantle convection. Through the solution of the convection equations, the models honor geophysical and geological data near the surface while predicting mantle ﬂow and structure at depth and associated surface deformation. The methods consist of four constraints determined a priori from a plate history model: (1) plate velocities, (2) thermal structure of the lithosphere, (3) thermal structure of slabs in the upper mantle, and (4) velocity of slabs in the upper mantle. These constraints are implemented as temporally-and spatially-dependent conditions that are blended with the solution of the convection equations at each time step. We construct Earth-like regional models with oceanic and continental lithosphere, trench migration, oblique subduction, and asymmetric subduction to test the robustness of the methods by computing the temperature, velocity, and buoyancy ﬂux of the lithosphere and slab. Full sphere convection models demonstrate how the methods can determine the ﬂow associated with speciﬁc tectonic environments (e.g., back-arc basins, intraoceanic subduction zones) to address geological questions and compare with independent data, both at present-day and in the geological past (e.g., seismology, residual topography, stratigraphy). Using global models with paleogeographical constraints we demonstrate (1) subduction initiation at the Izu-Bonin-Mariana convergent margin and ﬂat slab subduction beneath North America, (2) enhanced correlation of model slabs and fast anomalies in seismic tomography beneath North and South America, and (3) comparable amplitude of dynamic and


Introduction
Plate tectonics is the fundamental Earth sciences paradigm that provides a framework to interpret surface features and the geological record.Processes associated with plate motions and subduction of oceanic lithosphere form island arcs and volcanic belts, accrete material to continental margins, deform plate interiors, and drive vertical motions of ocean basins and continents.Subducting oceanic lithosphere that extends from present-day subduction zones into the upper mantle is revealed by earthquake locations and high seismic velocity anomalies by seismic tomography (e.g., Grand, 2002).Similarly, the circum-Pacific belt of high velocity seismic anomalies in the lower mantle can be correlated with slabs subducted during the Cenozoic and Mesozoic (Richards and Engebretson, 1992).A density model for subducted slabs can explain the degree 4-9 components of the observed long-wavelength geoid (Hager, 1984) and the negative buoyancy of slabs can drive present-day plate motions (Conrad and Lithgow-Bertelloni, 2002).
These first-order inferences suggest an intimate connection between the history of subduction and present-day mantle structure.Seismic tomography reveals lower-mantle slabs that can be used to determine the longitude of paleo-subduction zones and therefore potentially constrain absolute plate motions (van der Meer et al., 2010).Intra-Panthalassa subduction zones are predicted by combining the present-day position and timing of formation and accretion of extinct intra-ocean volcanic arcs with a plate reconstruction to compare with seismic tomography (van der Meer et al., 2012).Beneath the Americas, seismically fast anomalies around 800 km depth relate to subduction of the Nazca, Cocos, and Juan de Fuca plates (Ren et al., 2007).The broad sheet-like high velocity anomaly in the lower mantle beneath eastern North America is typically thought to have originated from the Cretaceous subduction of the Farallon plate (Liu et al., 2008;Grand, 2002;van der Hilst et al., 1997), although an alternative hypothesis suggests it is slab that originated from intra-oceanic subduction zones (Sigloch and Mihalynuk, 2013).
The large low-shear velocity provinces (LLSVPs) are long-wavelength structures (∼ 1000 km) at the core-mantle boundary that are positioned at present-day beneath Africa and the Pacific Ocean.The role of paleosubduction in determining the location and morphology of the LLSVPs remains debated.Some models favor substantial mobility of the LLSVPs in response to plate motions (Davies et al., 2012;Zhang et al., 2010;McNamara and Zhong, 2005), while others suggest that LLSVPs are insensitive to Wilson cycles (e.g., Torsvik et al., 2010) and may even organize plate tectonics (e.g., Dziewonski et al., 2010).Broad upwellings can be produced in dynamic models with imposed slab buoyancy flux for the past 300 Myr that share some similarities to the inferred distribution of plumes at the edges of LLSVPs (Steinberger and Torsvik, 2012).Therefore, understanding the dynamic interplay between surface tectonics and lower mantle structure enables us to interpret seismic images as a time-integrated record of an evolving thermochemical mantle.Furthermore, we can elucidate potential connections between surface geology and deep structure such as large igneous provinces derived from deep-seated mantle plumes.Convection models with data constraints are used to estimate global longterm sea-level change since the Cretaceous by predicting isostatic and dynamic topography to determine the volume of ocean basins (Müller et al., 2008), eustatic sea level (Spasojevic and Gurnis, 2012) and differential, vertical motion of continents (e.g., Moucha et al., 2008;Spasojevic and Gurnis, 2012).Present-day mantle thermal heterogeneity derived from seismic tomography is often used as an initial condition for inverse (backward advection) models (e.g., Spasojevic and Gurnis, 2012) or a present-day constraint for adjoint models (e.g., Bunge et al., 2003;Liu et al., 2008;Spasojevic et al., 2009).Stratigraphic data constrains the vertical motions of North America and provides an estimate of the viscosity ratio across the 660 km discontinuity by associating subduction of the Farallon slab with the widespread flooding of the western interior of North America (Spasojevic et al., 2009).
Convection models can be either entirely physics-based or semi-empirical.Entirely physics-based models honor the physics of fluid flow through conservation laws and are best suited to investigate the fundamental physics of convection in an Earth-like body.Using realistic constitutive relations, they can be used to explore evolving subduction in two dimensions and give rise to asymmetric subduction with slabs that are of the same thickness as the incoming oceanic plate (Burkett and Billen, 2009).However, given the high resolutions required, of the order of 1 km, physics-based convection models require enormous computational resources and so cannot yet be applied to large-scale, three-dimensional and long duration models that would be required to address the spatial and temporal characteristics of subduction preserved in the geological record.By contrast, semi-empirical models impose "known conditions" on the system to ensure that models are consistent with the history of subduction: for example, kinematic boundary conditions are commonly applied at the surface to model plate motions (e.g., Bunge et al., 1998).These models are not physically-self consistent because behavior that is enforced by applied conditions (e.g., plate motions) would not otherwise necessarily evolve naturally from the physics or parameters that the models include.Such models are still computationally expensive albeit feasible with existing technology.
The progressive data assimilation method that we develop here is a semiempirical approach that can be used to investigate the flow associated with specific tectonic environments, such as back-arc basins and intra-oceanic subduction zones.This enables us to address geological questions at both regional and global scales and compare with independent data at present-day and in the geological past (e.g., topography, geoid, gravity, seismic images, stratigraphy, rock uplift, etc.).Oceanic lithosphere, continents, slabs, and LLSVPs are volumetrically significant components of the mantle buoyancy field and are closely linked to plate tectonic history.Therefore, it is necessary to ensure that the temporal and spatial distribution of these buoyancy sources is consistent with plate history to create 4-D Earth models to compare with other data.Concurrent developments in paleogeographic software such as GPlates (Gurnis et al., 2012) are enabling construction of high temporal and spatial resolution plate history models (Seton et al., 2012) that include deforming regions (Flament et al., 2014).Variants of the assimilation method presented herein have already been incorporated into models to understand the paleogeography of Australia (Matthews et al., 2011), the influence of plate reconstructions on deep Earth structure (Bower et al., 2013), and the topographic asymmetry of the South Atlantic (Flament et al., 2014).

Method
We devise our progressive data assimilation method to produce global convection models with prescribed subduction zones that are consistent with an a priori plate history model.The method is comprised of four a priori data constraints that are applied to the convection model at each time step: (1) plate velocities (including velocities in deforming regions), (2) ther-mal structure of the lithosphere, (3) thermal structure of slabs in the upper mantle, and (4) velocity of slabs in the upper mantle.We denote nondimensional variables with primes and non-dimensionalize lengths with respect to Earth radius (R 0 , Table 1).Temperature is non-dimensionalized using T = ∆T (T + T 0 ), where ∆T is temperature drop and T 0 = T 0 /∆T is surface temperature.

Plate history model
A prescribed time-dependent plate history model is required to guide the evolution of the lithosphere and mantle in numerical models with progressive data assimilation.The plate model describes the evolution of plate boundaries, plate velocities, oceanic lithosphere ages, continent boundaries (i.e., non-oceanic regions), velocities in deforming regions, and internal boundaries such as cratons.Typically the model has a temporal resolution of 1 Myr.We generate a simple synthetic regional plate history model to demonstrate the assimilation method and then use a global plate history model modified and extended back in time from Seton et al. (2012) to demonstrate the method in spherical convection models.We create a sequence of a priori data files using the plate history model that are input to the mantle convection code at each time step.

Lithosphere assimilation
"Lithosphere assimilation" is a method to constrain the horizontal buoyancy flux of the lithosphere in a numerical model according to a plate history model.The method is comprised of two constraints; (1) plate velocities, and (2) thermal structure of the lithosphere.Plate velocities are often applied as a kinematic boundary condition to organize the convective flow in the mantle beneath the lithosphere (e.g., Han and Gurnis, 1999).Furthermore, plate velocities broadly steer the location of paleosubduction in models by promoting downwellings (slabs) to form at the margin of converging plates (e.g., Bunge et al., 1998).We also use the kinematic boundary condition to encode the deformation of continental lithosphere in addition to rigid plate velocities (e.g., Flament et al., 2014).
The thermal evolution of the lithosphere is constrained by modifying the thermal profile of the upper thermal boundary layer according to the plate history model.This also suppresses large-scale thermal instabilities away from downwelling regions and modulates the surface heat flux.For these reasons the method is unsuitable to investigate small-scale convection at the base of the lithosphere.We create a sequence of a priori data files from the plate history model of the thermal age, A, for the top surface nodes in the computational mesh for each time.The thermal age depends on both position (e.g., longitude and latitude) and time (t) (or age).For oceanic regions, the thermal age is usually determined from a model for the reconstructed seafloor age and we can additionally apply a maximum age constraint to model flattening of the seafloor.For non-oceanic regions, the age can be constant or vary according to the geological age (e.g., Archean, Proterozoic, Phanerozoic).
The thermal ages of the lithosphere are assimilated.At each time step in the convection model, we read in two of the a priori data files containing seafloor ages before and after the present model age.The thermal age for the present model time step is computed by linearly interpolating the age between the data files.We create an idealized thermal boundary layer for the lithosphere using thermal ages and the half-space cooling model: where T l is lithosphere temperature, T m mantle temperature, z depth, and τ thermal diffusion timescale (Table 1).Other cooling models for the lithosphere, such as the plate model, can be similarly implemented.At the end of each time step, we compute an updated temperature (T u ) by assimilating the idealized lithosphere thermal structure (T l ) with the advected temperature (T a ): where f l is the lithosphere assimilation factor that smoothly merges the plate model and advected temperature fields together: where z l is lithosphere assimilation depth and β lithosphere assimilation parameter (typically 0.5).Note that z l is not equivalent to the maximum depth of the thermal lithosphere and can either be constant (typically 0.01) or a function of lithosphere age (Eq.4-126 of Turcotte and Schubert, 2002): Furthermore, to investigate certain geological phenomena such as thermal subsidence following stretching and lithosphere-plume interaction we can construct "no assimilation regions" in which lithosphere assimilation is not applied.To achieve this we set z l to a small non-zero value for nodes in the computational mesh that are contained within a predefined no assimilation region (e.g., Flament et al., 2014).

Slab assimilation
Downwellings in models with kinematic boundary conditions do not always originate (nor remain fixed) at convergent margins and they are often symmetric drips, inconsistent with seismic observations of asymmetric subduction.Overcoming such artifacts would require very fine resolutions and more realistic rheologies that are currently not feasible in time-dependent global models.Therefore, we constrain the evolution of slabs in the upper mantle using "slab assimilation".This method enables numerical models to include asymmetric slabs that are consistent with data (e.g., dip angle) and prevent the not Earth-like shallow advective thickening of the over-riding plate that can occur when kinematic boundary conditions are applied.It additionally ensures that the thermal buoyancy flux of slabs in the upper mantle is consistent with the thermal buoyancy flux of the lithosphere prior to subduction at convergent margins in the plate history model.

General case
For slab assimilation we produce a priori data files that encode the slab model, which includes the thermal structure and velocity of slabs in the upper mantle.Subduction zone locations are exported from the digitized plate boundary dataset and stored as data files.The data files may also contain unique header data for each subduction zone that can be used to construct the slab model, such as slab dip.For each subduction zone we determine the position of the slab at each depth in the computational mesh by assigning the slab a radius of curvature in the uppermost mantle (R c ) and a constant dip in the uppermost mantle (θ) (Fig. 1A, Table 1).This generates a slab center line to guide the construction of the slab model.
The thermal age of subducting seafloor at convergent margins (A) is determined from the plate history model.We extrapolate this age along the slab center line to provide an estimate of the thermal age of the slab at each depth.We then create a thermal slab model by constructing a thermal boundary layer either side of the center line at each depth (Fig. 1B).The temperature drop across each boundary layer is T m /2 to ensure that the thermal buoyancy of the slab is equal to the thermal buoyancy of lithosphere of the same age.An additional sine term is necessary to conserve down-dip buoyancy since the thermal profiles are constructed at fixed depth rather than normal to slab dip.
where T s is slab temperature, θ slab dip, x horizontal coordinate from the center line, and other parameters as in Eq. 1. Regions away from slabs are assigned the mantle temperature (T = T m ) and the temperature at each depth is smoothed using a Gaussian filter (Fig. 2B).We additionally require the slab thermal model at each depth to be consistent with the thermal structure of the lithosphere, thus the "combined" slab temperature (T c ) is constructed similar to the lithosphere temperature profile (Eq. 1, Fig. 2C): We construct a thermal slab stencil (Ψ) to assimilate the idealized slab thermal structure (T c ) and the advected temperature field (T a ): where λ is a length scale that defines the lateral (x) and depth (z) extent of the stencil and µ is stencil smoothing.The stencil ranges from 0 for ambient mantle to 1 near the slab (Fig. 2A).We set λ x = 300 km and µ x = 50 km to define the lateral transition from the stenciled to non-stenciled region.Trial experiments with a smaller λ x revealed that the upper thermal boundary layer neighboring the subduction zone was unrealistically drawn down with the prescribed slab.Both λ z and µ z are determined independently for each subduction zone at each time because they are functions of the slab depth z s (t) that grows during subduction initiation (Fig. 5).We set λ z (t) = z s and additionally enforce a minimum (λ min = 75 km) and maximum (λ max = 350 km) stencil depth.Stencil smoothing (µ z ) increases linearly from 0.01 (small positive number) to µ max as the slab depth (z s ) increases from zero (surface) to λ max (Fig. 5).Note that the stencil depth (λ z ) formally defines the maximum depth at which the thermal slab model and advected temperature are combined in equal parts (0.5 contour, Fig. 2A,C).The stencil is smoothed with a Gaussian filter with a width of 110 km.
The combined slab temperature (T c ) and stencil (Ψ) are written to a priori data files that are read by the convection code at each time step.A linear interpolation routine calculates the temperature and stencil for the present time step similar to lithosphere assimilation.Slab assimilation follows lithosphere assimilation and the updated temperature is: where T f is final temperature with both lithosphere and slab assimilation, and T u is temperature with lithosphere assimilation (Eq.2).The final temperature (T f ) feeds into the right-hand side of the Stokes flow equation for the next iteration of the solver.

Shallow-dipping slabs
Slabs often do not subduct into the upper mantle with a single dip.For example, flat slab subduction is a leading hypothesis to explain the North American Laramide orogeny (e.g., Coney and Reynolds, 1977) and presently occurs at several subduction zones including the Andean margin (Gutscher et al., 2000).Non-Newtonian rheology and low viscosity weak zones can produce a flat slab in dynamic models (Van Hunen et al., 2004;Manea and Gurnis, 2007) but these are technically challenging to implement in global models.It is therefore useful to develop alternative methods to incorporate flat slabs in dynamic models to enable us to test and explore the implications for flat slab scenarios that may be constrained by independent data, such as the migration of volcanism on the over-riding plate.
We incorporate flat slabs in our slab assimilation method by modifying the geometry of the slab center line along which the temperature profile and stencil are constructed (Fig. 3).The digitized plate boundary dataset encodes flat slab location and lifetime, depth nearest the trench, and depth of the leading edge.This enables us to define a flat slab interface with shallow dip and we determine the (thermal) age along the interface using paleo-age grids from the plate history model.The thermal profiles for the steeply dipping segments of slab at the trench and leading edge are constructed as before (Fig. 1, Eq. 5).However, we modify the method for the slab at the leading edge to ensure that the thermal profile is only constructed for depths greater than the leading edge depth (Fig. 3).Along the shallow interface with the steeply dipping segments we use Eq. 5 with different definitions to construct the thermal profile: here, x is vertical distance to the center line and the sine term is obsolete (set θ = π/2).The slab stencil (Ψ = 1) is extended to include the shallow dipping interface using the location of the flat slab (Fig. 3).

Internal velocity conditions
We can further prescribe the velocity of the slab in the upper mantle by applying a velocity condition to select internal nodes in the computational mesh in the same manner as plate velocities are applied to surface nodes.We compute the along-strike and orthogonal-to-strike velocity components of the subducting plate for each subduction zone in the digitized plate boundary dataset using the plate kinematic model.The orthogonal-to-strike component is partitioned into vertical (depth) and horizontal components using slab dip.This enables us to construct a velocity vector for the slab in the upper mantle using the vertical, horizontal, and along-strike components.In cross-section, the slab velocity is applied to an internal node in the convection model that is closest to the slab center line (and within a certain distance) at one or several depths (e.g., Fig. 2A,C).Typically each depth is greater than 250 km (to avoid complications with the prescribed surface velocities) and less than the thermal stencil depth.In 3-D this produces a line of nodes with assigned velocities for a given depth level in the computational mesh.We can also coarsen the mesh that is used to export the velocities by factors of two in each dimension to downsample the number of prescribed velocity nodes to ensure continuity.Similar to lithosphere and slab assimilation, we construct a series of a priori data files that contain the three-components of slab velocity and the velocity stencil.The velocity stencil is equal to one for internal nodes that have a velocity applied ("on") and zero for free nodes ("off").Imposed velocity nodes turn on and off as subduction zone geometry in the plate history model evolves.At each time step in the numerical model we determine the velocity and stencil from the a priori file that contains data at the closest time (age) to the current model time (age).Prescribing the velocity of the slab may not always be appropriate, such as for investigations of slab break-off.

Solution method
At each time step in the convection code, data from the a priori plate history model at the current (convection) model age (Ma) are assimilated according to the following steps: 1. Solve the Stokes flow equation using plate velocities (and optionally, internal velocity conditions for the slabs) and with temperature (and optionally, composition) inserted in the right-hand side.2. Advect temperature (T a ) (and optionally, composition) using the flow velocity derived from the Stokes flow equation.3. Linearly interpolate the data-derived temperature fields for both the lithosphere (Eq. 1) and slabs (Eq.6) to the current model age using data at neighboring integer ages.4. Update temperature by assimilating the lithosphere (Eq.2). 5. Update temperature by assimilating the slabs (Eq.9).6. Update plate velocities to the current model age using linear interpolation.7. Optional: update internal velocity conditions for the slabs to the current model age (Section 2.3.3).8. Next time step (return to 1).

Convection models
We modify the finite-element code CitcomS (Zhong et al., 2000(Zhong et al., , 2008) ) to solve the equations for the conservation of mass, momentum and energy for incompressible flow (Boussinesq approximation) with our data assimilation method.CitcomS non-dimensionalizes length with Earth radius (Table 1) which means the Rayleigh number is approximately an order of magnitude larger than a definition with mantle depth.We apply kinematic and isothermal (T = 0) boundary conditions at the top surface and free slip and isothermal (T = 1) boundary conditions at the core-mantle boundary.
Our regional model domain is centered on the equator and spans 28.6 • (0.5 radian) in latitude (65 nodes), 57.3 • (1 radian) in longitude (129 nodes), and the depth of the mantle (2870 km, 65 nodes).Side boundaries are stress-free.The full sphere is constructed of 12 caps, each with 128 x 128 x 64 elements, giving a total of 12.6 million elements.Near the top surface the lateral resolution is about 45 km and the vertical resolution is about 18 km.For the regional cases (K, 1-8) and global case G1 the viscosity structure (η profile = 1, Table 2) is temperature-and pressure-dependent: where η 0 is 1 for 0-410 km depth, η T Z for 410-660 km depth, η LM for depths greater than 660 km, and η A is non-dimensional activation energy.For global cases G2 and G3 the viscosity structure (η profile = 2) is temperature-and composition-dependent (e.g., Flament et al., 2014): where η 0 is 1 for 0-660 km depth and η LM is 100 in the lower mantle.η C = 100 is the compositional pre-factor used for chemically distinct continents.
The viscosity is additionally truncated to be between 0.1 and 100.

Plate models and initial conditions
We create a regional synthetic plate history model from 100 Ma to 50 Ma with a temporal resolution of 1 Myr to show how the assimilation method can prescribe surface tectonic evolution in a numerical model.We define a ridge (thermal age=0 Ma) at 0 • longitude with a half-spreading rate of 5 cm/yr in the east-west direction.Prior to 100 Ma, a passive margin is specified east of the ridge with the same geometry as the subduction zone that evolves at later times.Subduction initiates at 100 Ma (the initial condition for the geodynamic model) and the trench is located at 46 • longitude (0.8 radians) for negative latitude and arcs northwestward for positive latitude (radius of curvature is approximately 1600 km, Fig. 4A).The trench is stationary from 100 Ma to 90 Ma and rolls back westward at 2 cm/yr from 90 Ma: displacement of the over-riding plate accommodates this motion (Fig. 4B,C).
We use a global plate history model from 230 Ma to present day that is extended from Seton et al. (2012) and includes a revised reconstruction of the Arctic (Shephard et al., 2013), Southeast Asia (Zahirovic et al., 2013), Gondwana breakup (Gibbons et al., 2013), and full-fit reconstructions for Australia-Antarctica (Whittaker et al., 2013).Furthermore, the model incorporates deforming continental regions such as the basin and range province of western North America, and the Laramide (∼100-50 Ma) and Andean (∼40-0 Ma) flat slabs.Encoding the model using GPlates ensures that plates (represented as polygons) are continuously closing and obey the rules of plate tectonics (Gurnis et al., 2012).GPlates outputs data files that contain subduction zone locations and properties (which can be time-dependent, such as slab dip) and plate velocities.Complementary age grids are built specifically for each reconstruction and provide the age of the lithosphere.
The global plate history model is used to construct the initial condition for the global cases (G1-3).G1 begins at 110 Ma and slabs are inserted to 660 km depth, and G2 and G3 start at 230 Ma and slabs initially extend to 1200 km depth.However, if subduction has recently initiated then slab depth is calculated using an upper mantle sinking rate instead.Slab dip is determined for each subduction zone from a database of slab properties that is compiled from seismic and geologic information (if available), and otherwise defaults to 45 • .A lower thermal boundary layer is included in G2 and G3 (mantle temperature, T m = 0.5) and is omitted for all other cases (T m = 1).Internal heating is not included.For all cases the initial thermal structure of the lithosphere is calculated from the age grids (Eq.1).The continental lithosphere is assigned a thermal age according to the geological age (e.g., Archean, Proterozoic, Phanerozoic) for G2 and G3 and is constant (200 Ma) for other cases.Furthermore, G2 and G3 include tracers to model continents (see Flament et al., 2014) and a deep Earth chemical layer of initial thickness 113 km and buoyancy number of 0.5.
Using GPlates, surface velocity at nodes within one degree of a plate boundary are smoothed linearly using the velocity at the node, the average velocity at the closest boundary point, and the distance between the node and the boundary.For the regional model the velocity of a slab in the upper mantle is prescribed at both 252 km and 336 km depth once the thermal stencil has grown to its maximum depth.Furthermore, we taper the surface velocity to zero at the edge of the domain and do not apply internal velocity boundary conditions within three nodes of the edge.For the global model G1 the slab velocity is prescribed at a single depth less than the slab depth when the slab is between 252 km and 336 km, and at 336 km depth once the thermal stencil has grown to its maximum depth.In addition, we coarsen the internal velocity boundary condition mesh by a factor of 2 to ensure continuity (Section 2.3).Slab velocities are not applied for cases G2 and G3.

Results
We demonstrate the assimilation method by varying parameters that govern the vigor and style of convection: Rayleigh number (Ra), radial viscos-ity structure (η T Z , η LM ), and activation energy (η A ) (Table 2).A typical Rayleigh number using an upper mantle reference viscosity (and Earth radius) is between 10 8 and 10 9 .The viscosity increase in the lower mantle likely ranges between 10 (Hager, 1984;Paulson et al., 2007) and 100 (Spasojevic et al., 2010;Forte and Mitrovica, 1996).Temperature-dependent viscosity variations in global models are usually less than 3-4 orders of magnitude to ensure numerical convergence.
Case 1 is our reference that we describe in detail and compare with the other cases.The initial temperature field (at 100 Ma) reveals oceanic lithosphere with increasing thickness with distance from the ridge and continental lithosphere with uniform thickness (Fig. 5A).Following subduction initiation (also at 100 Ma), the thermal stencil depth increases with the duration of subduction from λ min = 75 km to a maximum depth of λ max = 350 km as the slab depth increases (Fig. 5D,G).During this time, the stencil smoothing also increases to µ max = 75 km (Fig. 5D,G).The thermal stencil reaches its maximum depth by 89 Ma and the slab descent rate is applied to nodes at 252 km and 336 km depth that are closest to the slab center line for subsequent time steps (Fig. 5H,K,N).Slab rollback begins at 90 Ma (Fig. 5H) which produces a dipping slab in the lower mantle (Fig. 6A).Lithosphere assimilation prevents downwellings from occurring away from the subduction zone and also prevents the formation of a ridge by suppressing the upwelling at the edge of the domain beneath the continent (Fig. 6D).The slab develops smoothly and continuously from the prescribed stencil region into the lower mantle and is consistent with the first-order interpretation of slab geometry from seismology (Fig. 5M, 6A,D).Note that if subduction ceased the slab thermal and velocity stencil would instantaneously be removed and the structure of the slab would subsequently evolve solely by the solution of the convection equations.
We analyze the temperature, velocity, and temperature flux (proportional to buoyancy flux) of the seafloor prior to subduction (located 100 km outboard of the trench, Fig. 6D, black line) and post-subduction when the slab is at 336 km depth (the maximum depth that the slab velocity is applied, Fig. 6D, purple line) using an equatorial cross-section (Fig. 6A).We normalize the integrated slab profile (temperature and flux) by the integrated lithosphere profile at each time step to compare the relative magnitude of the temperature (T int ) and flux (Q int ); a value around one is expected.
where T sp (q sp ) and T lp (q lp ) are temperature (flux) profiles across the slab and lithosphere at a given time step, respectively, and p is distance along the profile.For time steps after the thermal slab stencil has reached its maximum depth extent following subduction initiation we compute the time-averaged means: T int and Q int .The standard deviation reports the fluctuation (Table 2).We also analyze the velocity field to determine if the magnitude of the applied slab velocity (= |v sub |) is greater ("M") or less ("L") than the average velocity outside the prescribed nodes (Table 2) .Case K has the same parameters as Case 1 except only kinematic boundary conditions are applied at the surface (no lithosphere or slab assimilation).In Case K the downwelling has a larger lateral length scale in the upper mantle (500 km versus 100 km in Case 1) and the depth of the slab is shallower, for example 625 km versus 800 km in Case 1 at 66 Ma (compare Fig. 6D  and 7A).The return flow at x ≈ 5800 km thermally erodes the continental lithosphere to a greater extent in models without lithosphere assimilation.Furthermore, for these parameters the magnitude of flow velocities around the slab are reduced when the slab descent rate is not applied as an internal velocity boundary condition (Fig. 7A).We determine T int = 2 and Q int = 0.3 at 67 Ma for a slab profile at 450 km depth and lithosphere profile 200 km outboard of the trench.In Case K the anomaly is too broad to use a profile at 336 km depth and 100 km outboard of the trench as for the other cases.
We increase the lower mantle viscosity in Case 2 to η LM = 100 which causes the slab to flatten at the 660 km interface (Fig. 7B).The thermal structure of the lithosphere (oceanic and continental) is almost identical to Case 1.In Case 3, relative to Case 1, we decrease the activation energy (η A = 2.303) which reduces the temperature-induced viscosity contrast by two orders of magnitude (Fig. 7C).Therefore, the lithosphere is more mobile and a second instability develops inboard of the trench at x ≈ 5100 km.Lithosphere assimilation with an assimilation depth of z l = 64 km is unable to suppress this instability and it forms a downwelling.Flow velocities are larger in Case 3 than Case 1.In addition, the slab descends deeper and has a steeper dip in the lower mantle, although these features might, in part, be due to the second instability that modifies flow.We reduce the transition zone viscosity in Case 4 (η T Z = 1) and the slab evolution is comparable to Case 3 (Fig. 7D).Around 58 Ma the slab detaches at approximately the thermal stencil depth (λ max ), although this is likely a consequence of negative pressure in the mantle wedge that is artificially accentuated by proximity to the boundary of the domain.The same phenomenon occurs in Case 3, albeit at later time (51 Ma).
The cases with a reduced Rayleigh number, Ra = 10 8 (relative to Ra = 10 9 for Cases 1-4) exhibit quite similar behavior, regardless of the pressureor temperature-dependence of viscosity (Fig. 7E-H).In all cases the slab flattens at the upper-lower mantle interface and does not penetrate into the lower mantle.In Case 7 (η A = 2.303, Fig. 7G), the slab has a slightly larger dip and does not extend as far laterally in comparison to the other Ra = 10 8 cases.
Case G1 demonstrates the application of progressive data assimilation to a global model.This model begins at 110 Ma and has Ra = 5.5 × 10 8 and the same viscosity structure as Case 3 and 7. We analyze two regions from this model: Laramide flat slab subduction beneath North America (Fig. 8) and the initiation of the Izu-Bonin-Mariana (intra-oceanic) convergent margin (Fig. 9).Laramide flat slab subduction begins around 100 Ma and the leading edge of the flat slab moves eastward and underrides North America (Figs.8A,B).Slabs that subducted prior to the formation of the flat slab are evident in the temperature field (e.g, cold material at 15 • , 700 km depth, Fig. 8C).The flat slab reaches its maximum lateral extent around 70 Ma (Fig. 8D).Izu-Bonin-Mariana subduction initiates around 52 Ma (Figs. 9A,C) and the slab descends to approximately 450 km depth by 39 Ma (Fig. 9G).Inboard of the trench are two spreading centers: the Philippine Sea mid-ocean ridge ("PSR", also referred to as the spreading center in the central region of the west Philippine Sea plate) and the proto-Izu-Bonin-Mariana (back-arc) ridge ("BAR").These ridges are revealed in the temperature cross-sections, for example "PS" at x ≈ 500 km and "BAR" at x ≈ 1500 km in Fig. 9A.
Relative to G1, cases G2 and G3 have a slightly larger Rayleigh number (8.6 × 10 8 ), a lower thermal boundary layer (T m = 0.5), chemically distinct continents (see Flament et al., 2014), and a different viscosity profile (Eq.11).G2 uses progressive data assimilation whereas G3 only has the surface kinematic boundary condition applied.The power spectra of temperature for G2 reveals relatively more power in degrees 4 and 6 in the lower mantle and less power in degrees 6-15 between 410 and 1000 km depth compared to G3 (Fig. 10D,F).G2 also has a reduced amplitude of thermal heterogeneity as revealed by the maximum power with depth (Fig. 10E,G) and slab buoyancy flux (Fig. 10C).The power spectra of S-wave seismic heterogeneity for S40RTS is shown for comparison (Fig. 10A).We also compare the dynamic topography at present day for G2 and G3 (Fig. 12) and cross-sections of temperature and tomography beneath North and South America (Fig. 11; see discussion).

Discussion
The objective of the new methods is to develop (mostly global) 4-D Earth models that are compatible with plate history by incorporating constraints on the thermal structure of the lithosphere and slabs.The regional cases demonstrate the key features of the assimilation method and the global cases illustrate applications to address geodynamical questions while demonstrating the shortcomings of the traditional use of only kinematic boundary conditions.

Regional cases
The synthetic regional models demonstrate that the assimilation method produces slab buoyancy flux that is consistent with plate history for a range of model parameters that dominantly control the vigor and style of convection: the Rayleigh number, and temperature-and pressure-dependence of viscosity.In all cases with assimilation, the slab in the upper mantle is narrow and asymmetric and retains its prescribed dip of 45 degrees independent of model parameters.Note that the slab is prescribed to 350 km depth but evolves solely according to the flow equations at depths greater than 425 km.In essence, the evolution of the uppermost part of the slab is independent of the local Rayleigh number.By contrast, Case K (Fig. 7A) has kinematic boundary conditions only and produces a slab with a large lateral length scale due to accumulated buoyancy that is caused, in part, by shallow advective thickening.The morphology and temperature distribution of this slab is notably different from the slab generated with data assimilation, which has a profound consequence for predicting surface observables such as dynamic topography and gravity.Steady trench roll-back produces an asymmetric slab in Case K, although models with kinematic boundary conditions often exhibit symmetric downwellings.
In most cases, lithosphere assimilation introduces an upper thermal boundary layer that is consistent with plate history and inhibits downwellings away from the subduction zone.The passive upwelling at the edge of the domain (x ≈ 5800 km, Fig. 7B-H) is suppressed at shallow depth by lithosphere assimilation which imposes the thermal structure of continental lithosphere (thermal age = 200 Ma).A secondary downwelling occurs in Case 3 (Fig 7C) which suggests that a larger or spatially-dependent (Eq.4) lithosphere assimilation depth (z l ) should be considered to stabilize the boundary layer.Therefore, z l is effectively a function of the local Rayleigh number and larger values are necessary for models without a stiff lithosphere.Large z l prevents flow in the convection model from modifying the thermal lithosphere, whereas small z l permits the thermal boundary layer to evolve according to model parameters which will result in a thick (low Ra) or thin (high Ra) thermal boundary layer that does not produce the dichotomy of continental and oceanic thermal structure.The models indicate that constant z l = 64 km is generally an appropriate choice for the range of convection parameters that we explore.
Assimilated slabs evolve continuously from the prescribed stencil into the transition zone and lower mantle (e.g., Figs.6D).Most of the higher Rayleigh number cases (Ra = 10 9 , Cases 1,3,4) produce slabs that descend into the lower mantle to ∼ 1000 km depth.Case 2 (Ra = 10 9 ) and the lower Rayleigh number cases (Cases 5-8) produce slabs that flatten at the upperlower mantle interface.Slab morphology is determined, in part, by trench retreat encoded in the plate history model.However, at later times the flow is almost certainly artificially influenced by the edge of the domain which causes the slab to flatten from enhanced negative pressure in the mantle wedge.Therefore we do not use these regional models to quantify the evolution of slab morphology, but simply note that the assimilation method can produce continuous slabs that exhibit a range of dip angles in the lower mantle.
The integral of temperature and temperature flux (proportional to buoyancy flux) across the lithosphere outboard of the trench and across the slab at 336 km depth are consistent with the prescribed plate history model (e.g., Fig. 6G,H).The time-averaged normalized integrated temperature (T int ) (Table 2) compares the temperature of the slab to the temperature of the subducting lithosphere and is slightly less than one (0.95 -0.97) for all models with assimilation except Case 3 which is 1.01: a value of one is expected for a rigid slab that subducts at the same speed as the plates are converging toward one another.This discrepancy may arise from warm return flow in the mantle wedge that is advected into the stencil and assimilated with the idealized slab thermal structure, which slightly reduces the temperature integral across the slab.
There is a variability of the time-averaged normalized integrated temperature flux (Q int ) between models (0.83-1.15).This is anticipated because Q int is computed from the temperature and velocity profiles of the lithosphere and slab.While the slab temperature profile is highly constrained by the thermal stencil, the slab velocity profile is only constrained by the velocity prescribed along two lines of nodes at 252 km and 336 km depth.Similarly, the lithosphere velocity profile depends on coupling to plate motions that are applied at the surface.Therefore, the velocity profiles are more sensitive to the convection model parameters (Ra, η, etc.) than the temperature profiles.For most cases, the internal velocity boundary condition increases the descent rate of the slab in the region of the prescribed nodal velocities in comparison to outside this region (Fig. 6E).In Case K, with kinematic boundary conditions only, the buoyancy flux is very small (Q int =0.3) which confirms that progressive data assimilation ensures that slab and lithosphere buoyancy fluxes are comparable.
We do not account for the dynamic coupling between temperature and velocity when we construct the a priori data files from the plate history model.This enables us to develop a simple and flexible method in which we determine the temperature and velocity constraints independent of each other.It is then straightforward to implement these constraints in a timedependent dynamic model by combining the idealized data with output from the previous time step.Consequently, the models are not entirely physically self-consistent.Appropriate model parameters should be selected such that the broad scale characteristics of flow are compatible with the plate history.Furthermore, we can avoid the need for internal velocity boundary conditions, which tend to be computationally expensive, by selecting model parameters that produce a slab sinking rate in the upper mantle comparable to average plate velocity.

Applications
A motivation for our data assimilation method is to enable geodynamic models to honor data constraints.The high seismic velocity anomaly beneath North America (Fig. 11C, at ∼ 80 • W) correlates well with a slab from G2 (with assimilation, pink contours) and poorly with a downwelling from G3 (surface velocities only, purple contours).Relative to the seismically-imaged slab, the downwelling in G3 is far too broad and displaced 20 • farther west.The downwelling in G3 is also too large beneath South America (Fig. 11D) and offset 15 • east compared to tomography.By contrast, data assimilation in G2 ensures that slab buoyancy flux is closely tied to the thermal structure of subducting lithosphere and thus produces slabs that compare more favorably in terms of size and location to high seismic velocity anomalies as revealed in tomography.Furthermore, the South American cross-section illustrates that much of the overriding plate is entrained in the mantle downwelling in the kinematic-only boundary condition case, but not in the slab assimilation case.The large spatial offsets between predicted cold anomalies from the model with only surface velocities prescribed and high seismic velocity anomalies (Fig. 11) illustrates that spectral analysis (Fig. 10) alone is insufficient to assess the success of a mantle flow model.
The large spatial offsets in the position of the slabs in models with and without assimilation is manifest in large shifts in the position of dynamic topography lows.The changing position of these lows has become an important geological constraint on geodynamic models (Flament et al., 2013).Furthermore, slabs are connected to the surface in models with slab assimilation which can be important to explain dynamic topography lows.We compare the dynamic topography of case G2 (with assimilation, Fig. 12A) and G3 (without assimilation, Fig. 12B) with the residual topography of Panasyuk and Hager (2000) (Fig. 12C).The dynamic topography predicted by the mantle flow models is generally of longer wavelength than the residual topography.This is particularly true of the positive dynamic topography because in the forward models the active mantle plumes have not yet reached the upper mantle (Flament et al., 2014).
The amplitude of the dynamic topography predicted for G2 (between -1670 m and 1470 m) is consistent with that of the residual topography (between -1700 m and 1120 m).In contrast, the amplitude of the dynamic topography predicted for G3 (between -6500 m and 2800 m) is much larger than the amplitude of residual topography (note that a scaling factor of 2 is used in Fig. 12B to visually compare the results of G2 and G3).The distribution of dynamic topography lows further illustrates the benefits of using slab assimilation for these calculations: G2 generally produces higher frequency dynamic topography lows than G3.For instance, G2 captures the two distinct lows occurring in residual topography to the northwest and southeast of India, and the negative dynamic topography extending from Arabia into northern Africa.The extent of the residual topography lows in eastern North America and between Australia and Antarctica (Fig. 12C) are also better captured by G2 than G3.
Fig. 10C shows the slab buoyancy flux for G2 (with assimilation, black lines) and G3 (surface velocities only, red lines).From 230 to 190 Ma the flux is low for G3 which is compatible with the regional models that show reduced flux for the cases with a kinematic boundary condition only relative to data assimilation (Table 2).During this time, downwellings in G3 are broad and do not penetrate the transition zone, akin to the thermal anomaly in Fig. 7A (see also Table 2, Case 10).Between 190 and 150 Ma the flux increases for G3 as broad downwellings impinge on the 660 km viscosity jump between the upper and lower mantle.Since the start of the model (230 Ma) the upper thermal boundary layer has thickened unabated because lithosphere assimilation is not applied.The downwellings source their buoyancy from this thick thermal boundary layer which explains why the amplitude of buoyancy flux is 4-7 times larger for G3 relative to G2 for much of the model run time (Fig. 10C, note different y-axes).Around 120 Ma the downwellings have accumulated sufficient buoyancy to descend through the lower mantle.This causes the mass flux to steadily decrease from 120 Ma to present day as the upper thermal boundary is drained (Fig. 10C).
In summary, models with surface velocities only may be characterized by thickening of the upper thermal boundary layer followed by sinking of cold material in the mantle.Therefore, downwellings in these models are not slab-like because buoyancy is not smoothly and continuously introduced into the mantle throughout the duration of the model (∼ 200 Myr).By contrast, models with assimilation (e.g., G2) prevent large variation of slab buoyancy flux (unless dictated by data) and over-thickening of the upper thermal boundary layer.This explains why the amplitude of thermal heterogeneity (Fig. 10E,G) and dynamic topography (Fig. 12), and size of thermal anomalies (Fig. 11), is less for G2 than G3.

Conclusions
Progressive data assimilation methods are powerful techniques to ensure the spatial and temporal distribution of buoyancy in a numerical model is consistent with plate history.This enables us to test geological hypotheses by comparing model outputs (e.g., dynamic topography, mantle structure) with observational constraints (e.g., residual topography, seismology).The methods apply plate motions and velocities in deforming regions at the surface and produce a thermal lithosphere consistent with seafloor age (and optionally, age of continental regions) and asymmetric narrow slabs at subduction zones.The thermal structure of slabs is constrained by the (thermal) age of subducting lithosphere to ensure buoyancy conservation.Furthermore, the slab buoyancy flux is constrained by additionally prescribing the velocity of slabs in the upper mantle, although this may not be necessary if convection models inherently produce slab sinking rates comparable to average plate velocities.Synthetic regional models show the methods are robust for a range of viscosity laws and Rayleigh numbers.Global models demonstrate application of the methods to understand mantle structure and residual topography.(A) Schematic of slab geometry (cross-section) with upper mantle descent velocity applied to the node at P 2 .For visual clarity the slab depth (z s ) is the same as the maximum slab stencil depth (λ max ).The thermal slab is constructed at each horizontal layer in the computational mesh about a center line.(B) Thermal profile at fixed depth.See Table 1 for parameters.Figure 2: Cross-section of the slab thermal and velocity stencil, slab temperature, and combined temperature using parameters derived from the synthetic regional plate history model (Section 3.2, Table 1).(A) Slab thermal stencil (Ψ) and velocity stencil.Grey circles are nodal points for a computational mesh with 50 km resolution and the slab velocity in the upper mantle is applied to the node at P 2 (large purple circle).2).(Panasyuk and Hager, 2000).Oceanic regions are water-loaded and continental regions are air-loaded in (A), (B) and (C).
Figure 1: (A) Schematic of slab geometry (cross-section) with upper mantle descent velocity applied to the node at P 2 .For visual clarity the slab depth (z s ) is the same as the maximum slab stencil depth (λ max ).The thermal slab is constructed at each horizontal layer in the computational mesh about a center line.(B) Thermal profile at fixed depth.See Table1for parameters.

Figure 3 :Figure 4 :Figure 5 :
Figure2: Cross-section of the slab thermal and velocity stencil, slab temperature, and combined temperature using parameters derived from the synthetic regional plate history model (Section 3.2, Table1).(A) Slab thermal stencil (Ψ) and velocity stencil.Grey circles are nodal points for a computational mesh with 50 km resolution and the slab velocity in the upper mantle is applied to the node at P 2 (large purple circle).(B) Slab temperature.(C) Combined temperature with slab thermal stencil contours and applied slab velocity.

Figure 6 :
Figure6: Case 1 at 66 Ma. (A) Equatorial cross-section of temperature with velocity, (B) Plate velocity at surface, (C) Temperature of uppermost mantle with theoretical isotherms (T = 0.25, 0.5, 0.75) computed from the age distribution and the half-space cooling model shown as black lines.(D, E, F) Temperature, magnitude of velocity (in plane), and viscosity, respectively, with contours of the thermal stencil (Ψ = 0.1, 0.5, 0.9) shown as green lines and isotherms (T = 0.25, 0.5, 0.75) shown as black lines.The thick black and purple lines show the location of profiles across the lithosphere and slab, respectively.(G, H, I) Profiles of temperature, flux, and velocity, respectively, across the lithosphere and slab.We integrate the profiles at each time step to compute T int and Q int (Table2).

Figure 12 :
Figure12: Present-day dynamic topography (meters) computed for sources deeper than 350 km (as inFlament et al., 2014) with free-slip boundary conditions for (A) G2 and (B) G3.For G3 the dynamic topography has been divided by a factor of two to enable comparison with G2. (C) Residual topography(Panasyuk and Hager, 2000).Oceanic regions are water-loaded and continental regions are air-loaded in (A), (B) and (C).

Table 1 :
Input parameters for a priori assimilation data files and convection models.† The second value applies to cases G2 and G3 only.

Table 2 :
Assimilation model parameters and results.† Influence of internal velocity boundary conditions on slab descent."M" means applied velocity is generally greater than total slab descent rate as determined by model parameters.The opposite is true for "L".‡ Surface velocities only.Upper mantle / lower mantle.