Effects of upper mantle heterogeneities on the lithospheric stress field and dynamic topography

The orientation and tectonic regime of the observed crustal/lithospheric stress field contribute to our knowledge of different deformation processes occurring within the Earth’s crust and lithosphere. In this study, we analyze the influence of the thermal and density structure of the upper mantle on the lithospheric stress field and topography. We use a 3-D lithosphere–asthenosphere numerical model with power-law rheology, coupled to a spectral mantle flow code at 300 km depth. Our results are validated against the World Stress Map 2016 (WSM2016) and the observationbased residual topography. We derive the upper mantle thermal structure from either a heat flow model combined with a seafloor age model (TM1) or a global S-wave velocity model (TM2). We show that lateral density heterogeneities in the upper 300 km have a limited influence on the modeled horizontal stress field as opposed to the resulting dynamic topography that appears more sensitive to such heterogeneities. The modeled stress field directions, using only the mantle heterogeneities below 300 km, are not perturbed much when the effects of lithosphere and crust above 300 km are added. In contrast, modeled stress magnitudes and dynamic topography are to a greater extent controlled by the upper mantle density structure. After correction for the chemical depletion of continents, the TM2 model leads to a much better fit with the observed residual topography giving a good correlation of 0.51 in continents, but this correction leads to no significant improvement of the fit between the WSM2016 and the resulting lithosphere stresses. In continental regions with abundant heat flow data, TM1 results in relatively small angular misfits. For example, in western Europe the misfit between the modeled and observation-based stress is 18.3. Our findings emphasize that the relative contributions coming from shallow and deep mantle dynamic forces are quite different for the lithospheric stress field and dynamic topography.


Introduction
The stresses building up in the rigid outermost layer of the Earth are the result of both shallow and deep geological processes. The dynamics of the lithosphere is determined by a combination of plastic, elastic and viscous flow properties of the lithospheric material (Burov, 2011;Tesauro et al., 2012), while the evolution of the sub-lithospheric mantle is predominantly driven by viscous flow (Davies, 1977;Forte and Mitrovica, 2001;Steinberger and Calderwood, 2006). It has been shown that shallow processes influence both the magnitude and orientation of the lithospheric stresses. Among such processes the most important are slab pull, ridge push, trench friction and continental collision (deformation) (Reynolds et al., 2002) as well as cratonic root resistance (Naliboff et al., 2012). Gravitational effects due to lateral density heterogeneities in the lithosphere and tractions from the mantle flow at the base of the moving plates also play an important role. Superposition of different tectonic forces creates dissimilar orientations and regimes of the lithospheric stress field in different regions, as shown by the World Stress Map Published by Copernicus Publications on behalf of the European Geosciences Union.

650
A. Osei Tutu et al.: Mantle heterogeneities, lithospheric stress and dynamic topography project (Bird and Li, 1996;Heidbach and Höhne, 2007;Heidbach et al., 2010Heidbach et al., , 2016. Furthermore, on a global scale the intra-plate stress orientation follows a specific pattern at a longer wavelength due to a large force contribution from the convecting mantle (Steinberger et al., 2001;Lithgow-Bertelloni and Guynn, 2004). This first-order stress pattern (long wavelength) is dynamically supported, as the controlling forces correlate well with the forces driving the plate motion in most continental areas such as North and South Americas and Europe (Solomon et al., 1980;Richardson, 1992;Zoback, 1992). Ghosh and Holt (2012) and Steinberger et al. (2001) used different approaches to show that the contribution of the crust (shallow density structures) to the overall lithospheric stress pattern is rather small compared to that of the mantle buoyancy forces, amounting to ∼ 10 %, except for regions characterized by high altitudes, especially the Tibetan Plateau, where the contribution is larger. In these previous modeling studies, the effect of the crust was determined separately by computing the gravitational potential energy from a crust model, which was subsequently applied as a correction (Steinberger et al., 2001;Ghosh et al., 2013;Ghosh and Holt, 2012). The contribution of the crust with a shallow lithospheric density contrast generates the second-order pattern (mid-to-short wavelength) in the stress field, mostly coming from topography and crustal isostasy (Zoback, 1992;Zoback and Mooney, 2003;Bird et al., 2006).
Likewise, the long-wavelength signal of the topography is related to the vertical component of the stress field tensor originating from the thermal convection of the mantle rocks (Pekeris, 1935;Steinberger et al., 2001). This generates a high dynamic topography in regions of upwelling over the African and Pacific large low shear velocity provinces (LLSVP) and low topography above downwelling in the regions of subduction (Hager and O'Connell, 1981;Hager et al., 1985). In contrast, at a mid-to-short wavelength, topographic features are influenced by processes such as plumelithosphere interaction (Lithgow-Bertelloni and Silver, 1998;Thoraval et al., 2006;Dannberg and Sobolev, 2015) and small-scale convection in the upper mantle (Marquart and Schmeling, 1989;King and Ritsema, 2000;Hoggard et al., 2016). However, the largest fraction of topography is caused by isostasy due to variations in crustal thickness and density, as well as density variations in the subcrustal lithosphere. Comparatively, the fit typically obtained between the modeled dynamic and observation-based residual topographies is lower (Heine et al., 2008;Flament et al., 2012;Steinberger and Calderwood, 2006;Steinberger, 2016;Hoggard et al., 2016), than for the other mantle-flow related observables. For example, the modeled and the observed geoid models give a relatively higher correlation (Čadek and Fleitout, 2003;Hager et al., 1985;Richards and Hager, 1984) due to the large contribution of the lower mantle but the modeled geoid is sensitive to the choice of mantle viscosity (Thoraval and Richards, 1997). One of the reasons for the poor correlation between modeled and residual topographies is our insufficient knowledge of the petrological properties of the upper mantle (Cammarano et al., 2011), for example in relation to the chemical depletion of cratons in continental regions. Hence, in this study, in addition to evaluating the influence of thermal-density heterogeneities on lithosphere stress field and topography, we test the impact of corrections for the continental depletion on the topography and stress field.
Also, constraining the modeled lithospheric stress with observations is challenging due to previously poor spatial coverage by World Stress Map data (Zoback, 1992;Lithgow-Bertelloni and Guynn, 2004;Heidbach et al., 2010). An alternative method documented in the literature is to compare the strain rate estimated from the modeled deviatoric stresses (Ghosh et al., 2008) with the Global Strain Rate Map (Kreemer et al., 2003). However, the lithospheric stress in plate interiors (i.e., far from the plate boundaries) is not well constrained with the Global Strain Rate Map. Hence, a gradually increasing coverage of the observed global stress field data serves as a motivation for studies attempting a global comparison of the observed and modeled stress field patterns, including our present study.
To date, two distinct approaches have been adopted to study the origin of the lithospheric stress, and each has given a relatively good fit to the observed stress field. On the one hand, Bird et al. (2008) have estimated the lithospheric stress from a model that disregards the mantle flow contribution and used the fit between modeled and observed plate velocities as a sole criterion. On the other hand, Ghosh et al. (2013), Ghosh and Holt (2012), Lithgow-Bertelloni and Guynn (2004), Steinberger et al. (2001) and Wang et al. (2015) have aimed at assessing the influence of mantle flow on the lithospheric stress field and have shown that the bulk mantle flow explains a large part (about 80-90 %) of the stress field accumulated in the lithosphere (Steinberger et al., 2001), in both magnitude and the most compressive horizontal direction. The aim of the present study is to evaluate the contribution of the upper mantle density and viscosity heterogeneities above the transition zone to the observed spatial stress regimes of the lithosphere (Heidbach et al., 2016), while testing different approaches and data sets used to describe the thermal and rheological structure of the upper mantle and the crust. We use a 3-D global lithosphere-asthenosphere finite element model (Popov and Sobolev, 2008;Sobolev et al., 2009) with visco-elasto-plastic rheology coupled to a spectral model of mantle flow (Hager and O'Connell, 1981) at 300 km depth. Deriving all force contributions from a single calculation resolves any inconsistency that might arise from treating individual force contributions to the stress field separately, as has been done in earlier studies (Bird et al., 2008;Steinberger et al., 2001;Lithgow-Bertelloni and Guynn, 2004;Ghosh et al., 2008;Naliboff et al., 2012;Ghosh et al., 2013;Wang et al., 2015). As part of this work, we further estimate dynamic topography and correlate our results with two different residual topogra-A. Osei Tutu et al.: Mantle heterogeneities, lithospheric stress and dynamic topography 651 phy models. One is based on seismic surveys of the ocean floor used to correct for shallow contributions to topography and free-air gravity anomalies on continents (Hoggard et al., 2016). The second model is taken from Steinberger (2016) and is based on actual topography corrected for crustal thickness and density from CRUST1.0 (Laske et al., 2013). Both models are also corrected for subsidence of seafloor with age.

Model description
Our global numerical model of the Earth interior consists of the particle-in-cell finite element model SLIM3D (Popov and Sobolev, 2008) within the top 300 km, which solves coupled momentum and energy equations with a semi-Lagrangian Eulerian grid with a free top boundary condition and a Winkler dynamic bottom boundary condition; this is coupled to a spectral mantle flow code (Hager and O'Connell, 1981) to account for the deep mantle contributions. There is no material exchange across the coupling interface, but the continuity of tractions and velocities is ensured through the Newton-Raphson iteration method. Figure 1c shows a sectional schematic representation of the coupled numerical model with depth-dependent layered mantle viscosity structure (Fig. 1b) and the seismic velocity-to-density scaling profile of Steinberger and Calderwood (2006) (Fig. 1a), which are only considered below the depth of 300 km. The top thermo-mechanical component (SLIM3D) has been used in a wide range of 2-D and 3-D regional numerical studies of crustal and lithospheric deformations (Popov and Sobolev, 2008;Brune et al., 2012Brune et al., , 2014Brune et al., , 2016Quinteros and Sobolev, 2013) with different spatial and temporal resolutions, but the coupled code is used here and in Osei Tutu et al. (2018) for the first time. In this 3-D global study, we distinguish three material layers (phases) within the top component (SLIM3D), the crustal layer, the lithosphere and the sub-lithospheric mantle layers, in order to account for the stress and temperature-dependent rheology in the presence of major continental keels and the uppermost part of the subducted lithospheric plates. The visco-elastoplastic rheology is described in detail by Popov and Sobolev (2008), with specific modeling parameters given in Osei Tutu et al. (2018) and in the Appendix of this paper.
This study complements our previous study (Osei Tutu et al., 2018) about the influence of plastic yielding at plate boundaries on plate velocities in a no-net-rotation reference frame and on lithospheric net rotation. A forward model is run for half a million years with a time step of 50 kyr, and at each time step tractions in the lower mantle due to density heterogeneities are computed using the spectral mantle code and then passed across the coupling dynamic boundary to the top component SLIM3D. Within the upper domain (SLIM3D), the flow velocities are then computed and passed back across the coupling boundary as an upper boundary condition to the spectral mantle code, with the method convergence estimated by comparing the velocity and traction norms of two successive iterations. Within the upper mantle, our crustal rheology is taken from Wilks (1990) and below the crust we have considered dry and wet olivine parameters in the lithosphere and sub-lithospheric mantle layers, respectively, modified after the axial compression experiments of Hirth and Kohlstedt (2004) (shown in the Appendix, Table A1. Adopted from Osei Tutu et al. (2018) for studying the influence of both the driving and resisting forces that generate global plate velocities and lithospheric plate net rotation).

Thermal and density structures of the upper and lower mantle
We assign densities of the uppermost layers according to the crustal model CRUST1.0 (Laske et al., 2013). Underneath, we separately consider the layers below and above the interface between the two codes placed at a depth of 300 km to differentiate between the deep and shallow signals. Here the topographic signal induced by the layers below 300 km is assumed to be due to convection in the viscous mantle, although cold rigid subducting slabs (Zhong and Davies, 1999;Faccenna et al., 2007) and possibly also the deepest cratonic roots (Conrad and Lithgow-Bertelloni, 2006) extend deeper than 300 km. We use a 3-D density structure inferred from the hybrid seismic tomography model of Becker and Boschi (2002) and apply a velocity-to-density conversion profile ( Fig. 1a) for the lower mantle buoyancy. In the upper mantle we test two different models for the representation of the upper mantle thermal and density structures, namely TM1 ( Fig. 1d) and TM2 (Fig. 1e). TM1 is based on the 3-D thermal structure TC1 model (Artemieva, 2006) across continents. This is combined with a 3-D thermal structure inferred from the seafloor age  for the mantle under oceanic regions. We use a half-space cooling model to infer the temperature T ocean as a function of age and depth according to the following: where k = 8 × 10 −7 m 2 s −1 is the thermal diffusivity, τ is the age of the oceanic lithosphere, T s is the reference surface temperature, T m is the reference mantle temperature and z is the depth beneath the Earth's surface. In regions of continental shelf, where there is neither age grid nor heat flow data, we interpolated the resulting thermal structures surrounding these regions while in Iceland where both age grid and heat flow data exist, the TC1 model was assigned. We explicitly include slabs in TM1 as described in Osei Tutu et al. (2018). The second model of the upper mantle thermal structure (TM2) is inferred from the seismic tomography model SL2013sv (Schaeffer and Lebedev, 2013). We have chosen this model because of its detailed representation of the up-  (Steinberger and Calderwood, 2006) and (c) a schematic diagram of the numerical method that couples the 3-D-lithosphereasthenosphere code SLIM3D (Popov and Sobolev, 2008) to a lower mantle spectral flow code (Hager and O'Connell, 1981) at a depth of 300 km. Panels (d) and (e) show the thermal structure at a depth of 80 km from two 3-D thermal models adopted in this study. (d) TM1, a heat flow-based thermal structure inferred from the TC1 model of Artemieva (2006) in the continents and the seafloor age model of Müller et al. (2008) in the oceanic areas. (e) TM2, the thermal structure of the upper mantle inferred from the S-wave tomography model SL2013sv of Schaeffer and Lebedev (2013). The "ringing" visible in panel (d) is an artifact introduced by smoothing sharp boundaries with a spherical harmonic expansion.
per mantle heterogeneities, which has been shown by Steinberger (2016) to allow for a better prediction of the dynamic topography than in previous models. This makes it a good candidate for comparison with the model results obtained using TM1, and for a regional investigation of the upper mantle contribution to the lithospheric stresses and topography.
Here we convert seismic velocity anomalies δV s into thermal anomalies T within the upper mantle according to the following relation: where the subscript P stands for a partial derivative at constant pressure (i.e., depth) based on Steinberger (2007). Firstly, we do not correct for the effect of the chemical depletion in cratons in order to evaluate its influence on the modeled lithospheric stress field and topography. In addition and for comparison purposes, we introduce two other thermal models based on two different seismic tomography models SAW24B16 (Mégnin and Romanowicz, 2000) and S20RTS (Ritsema et al., 2011) to evaluate their performance relative to our reference seismic tomography model SL2013sv (Schaeffer and Lebedev, 2013) ( Fig. 1e). We show in Fig. S1 in the Supplement different depth slices of TM1 and TM2. In our model setup, we define the reference crustal, lithospheric and asthenospheric densities (Table A1) and account for lateral density variations linked to thermal anomalies ( Fig. 1) using the following relation: where ρ ref denotes the reference density at a reference temperature of 20 • C and zero pressure, α denotes the thermal expansion coefficient chosen to be 2.7 × 10 −5 K −1 in the crustal layer and 3 × 10 −5 K −1 within the lithospheric and asthenospheric mantle and K is the bulk modulus (Table A1). We also compared it to the geoid estimate from the simulation using a layered/radial viscosity structure (Steinberger and Calderwood, 2006) for all depths (see Supplement Fig. S2c) resulting in a somewhat lower correlation of 0.82. In Fig. 2d, we show profiles of the estimated effective creep viscosity for continents and oceans within the upper 300 km using olivine parameters modified after Hirth and Kohlstedt (2004) (Shown in the Appendix, Table A1). Fig-ure 2d shows how a laterally averaged (dependent on depth only) asthenospheric viscosity decreases with increasing water content (i.e., 100, 500, 1000 H/10 6 Si). The viscosities are averaged separately across the continental and oceanic regions ( Fig. 2d, dashed versus solid lines). The average oceanic viscosity profiles give lower values than the respective average continental viscosity (η eff ), within the depth range of 100 ± 60 km. Seismological studies (e.g., Schaeffer and Lebedev, 2013;Kawakatsu et al., 2009;Fischer et al., 2010;Rychert et al., 2005) show this as a seismic wave velocity drop (∼ 5-10 %), and as a transition between the lithosphere and asthenosphere corresponding to the low viscosity channel (Fig. 2d). Figure 2d shows tractions causing stresses and topography in the lithosphere from the simulation using the creep parameters that correspond to the green effective viscosity profile in Fig. 2d, which was used to model the dynamic geoid and plate motions. Here, at all plate boundaries we have used a friction coefficient µ = 0.02 within the crust and lithospheric layers to generate the global plate velocities in a no-net-rotation (NNR) reference frame shown in Fig. 2b with a RMS of 3.5 cm yr −1 . Since the focus of this study is to investigate the effect of the upper mantle lateral density variations on the horizontal stress field and dynamic topog- Tutu et al.: Mantle heterogeneities, lithospheric stress and dynamic topography raphy, an assessment of the influence of the plate boundary friction and water content in the asthenosphere on plate velocities has been carried out in a separate study (Osei Tutu et al., 2018). In the present work, we therefore constrain our resulting creep viscosity with a cutoff for extreme viscosity values in the upper mantle by setting permissible minimum and maximum viscosity values similar to  and Osei Tutu et al. (2018); this approach yields a good fit between the observed and modeled geoid.

Shallow and deep contributions to the crustal stress state
We start by examining the separate contributions of the mantle heterogeneities below (deep Earth setup) and above (shallow Earth setup) 300 km to the global lithospheric stress field and topography. To calculate the contribution of the lower domain, we use a constant lithosphere thickness (100 km) and density (3.27 kg m −3 ), with the same configuration of the mantle below as was used to derive the geoid anomaly, plate motions and shear tractions in Fig. 2a-c. The resulting maximum horizontal magnitude (SH max ) and direction of the lithospheric stress field are shown in Fig. 3a.
We have obtained compressional regimes in regions of past and present subduction. On the North and South American continents, beneath which the ancient Farallon and Nazca plates were subducted, compressive stress magnitudes reach about 40 MPa. In the Far East, downwelling flows stretching from north to south from the northwestern Pacific through Australia towards Antarctica create compressional stress regimes with magnitudes ranging between ∼ 50 and 80 MPa. These compressional regions are connecting the Arctic with Antarctica and engulf two distinct regions with extensional stress regimes centered on the Pacific and African superswell regions. The predicted SH max directions in Fig. 3a generally follow the first-order lithospheric stress pattern (Zoback 1992), similar to previous mantle flow predictions of lithospheric stresses (Steinberger et al., 2001;Lithgow-Bertelloni and Guynn, 2004;Ghosh and Holt, 2012). In the largest extensional regions such as those found in the Pacific Superswell and above deep upwellings across southeastern Africa, stresses reach magnitudes of around 30 MPa. The modeled extensional/compressional patterns in the constant lithosphere, which get smoothed out over large distances are induced by the gradient in the tractions (Fig. 2c) coming from the mantle flow.
To investigate the contribution of the upper domain (300 km) to the stress field, we calculate the magnitude and direction using model TM1 (Fig. 1d) combined with the CRUST 1.0 model (Laske et al., 2013) and disregard mantle density variations below 300 km (i.e., both horizontal and vertical tractions below the depth of 300 km are set to zero). Comparison of the lithosphere stress predictions from our shallow (Fig. 3a) and deep (Fig. 3b) Earth setups reveals notable differences in the model-based stress regimes, mag- nitudes and directions in continental regions. If stresses are generated by the upper domain only, then almost all continental regions are characterized by an extensional regime, with the largest stress magnitudes found in areas of high topography and orogenic belts, such as the Tibetan and Andean highlands. Our stress predictions from the shallow Earth setup with laterally varying crustal and lithospheric densities in Fig. 3b show stress magnitudes and patterns similar to Naliboff et al. (2012). However, as opposed to the results of Naliboff et al. (2012) we predict high compressional stress magnitudes at continental margins, which may in part originate from a finer treatment of the crust and our temperature-dependent creep viscosity. Also the high compressional stresses along the subduction margins in Fig. 3b are likely induced by the slab models included in our setup. The resulting topography beneath air (free surface) accompanying either modeled lithospheric stress field is shown in Supplement Fig. S3a-

Total lithospheric stresses and topography
Next we compute the combined effect of both the lower mantle buoyancy and the upper mantle heterogeneities on the global SH max magnitude and direction for comparison with the separate contributions discussed above and with observations. Note that this is not a linear superposition of the separate contributions, because changes in the properties of the upper 300 km lead to changes in the topography and stress caused by density anomalies below 300 km depth. The resulting SH max direction and magnitude (Fig. 4a) due to the combined contributions of the upper and lower mantle show compressional regimes in areas similar to Fig. 3a, while muting almost all strong extensional stresses predicted by our simulation with the shallow Earth setup in continents (Fig. 3b).
The predicted SH max orientation also generally follows the first-order lithospheric stress pattern (Zoback, 1992), similar to predictions based on only density anomalies below the 300 km depth (Fig. 3a), with some regional deviations. The dominance of the contribution from below 300 km to the lithospheric stress field orientation becomes apparent when looking at the similarities between the SH max directions in Figs. 3a and 4a and dissimilarities with Fig. 3b, especially within continents. Nevertheless, the contribution from the upper 300 km to the predicted stress magnitude is evident in areas with large crustal thickness in continents, such as Tibet and the Andes. The regions where extensional regimes are predicted with only the contribution from below 300 km (Fig. 3a) correspond well with the extensional stress regions in the combined model (Fig. 4a). The result is not very different, when we use the thermal density model TM2 for the total lithospheric stress field prediction. Both the predicted SH max magnitude and direction with TM1 (Fig. 4a) and TM2 (Fig. 3b) show notable similarities in oceans and continents, owing to the strong contributions from below 300 km which are similar for both models. They show relatively high compressional stress magnitudes in subduction or convergence regions such as the Mediterranean, south of the Tibetan Plateau, south of Alaska, and the northwest Pacific extending through the Sumatra subduction zone and underneath the Australian and Antarctic plates. However, the SH max compressional signal underneath North America in Fig. 4a is muted and that of the South American region turns into an extensional regime along the Andes (Fig. 4) with the inclusion of the mantle above 300 km and the crust. Similar to Fig. 3a both predictions with TM1 and TM2 ( Fig. 4a and b) show SH max extensional regimes corresponding to the regions of upwellings and/or volcanism. However, the model with TM2 generates a much higher extensional magnitude of ∼ 60 MPa in the North Atlantic region around Iceland, and around the Azores and Canary hotspots, compared to TM1. Stress magnitudes are more alike in the southern Pacific Rise and around southern Africa. Differences are in part due to the detailed and well resolved upper mantle structures in the S-wave model used to derive TM2 (Schaeffer and Lebedev, 2013), as opposed to the upper mantle structure in TM1, which is based on the seafloor age in oceanic regions  and slab temperatures from Steinberger (2000b). Also in regions where the coverage of heat flow data is poor (e.g., in South America and Antarctica, Artemieva, 2006;Pollack et al., 1993), TM2 (Fig. 3b) may give better results. TM2 predicts compressional stress under Antarctica and along the subducting Nazca Plate in South America induced by downwelling flow. In these regions there are barely any heat flow data and TM1 remains largely unconstrained. Both modeling setups with the combined effects from the crustal structure model, the upper mantle thermal-density structure (either TM or TM2) and deep mantle contributions give topography (Fig. 4c-d) comparable to actual topography (see Supplement Fig. S4ab).
Lastly, we tested the bending of stresses inside plates with calculations in which the effect of elasticity is set to zero. This was compared with similar stress calculations considering elasticity and we found that elasticity has greater influence on the modeled stress magnitude compared to the corresponding orientations as given in Supplement nental margins and at the foot of the Andes whereas they are extensional above subduction zones for example Izu-Bonin-Mariana (Supplement Fig. S5).

Modeled versus observed lithospheric stress field
We compare our predicted SH max orientation to the observational stress data. Following the stress interpolation method presented by Müller et al. (2003), we used their fixed search radius (FSR) method which uses a global weighting defined by a fixed Euclidean distance for the stress data interpolation and stress quality. The smoothed stress field orientation at a grid point is based on the dominant stress data orientation within the selected radius. For a detailed explanation on the FSR method see Müller et al. (2003). Stress data with quality A, B and C with known stress regime were considered. Since we do not consider the respective regime in our quantitative analysis, we also included the stress data with unknown style having quality A and B in our smoothing procedure to make our smooth field more robust. We smoothed the observed SH max orientation of the World Stress Map 2016 (WSM2016) (Heidbach et al., 2016), with a search radius of 270 km (Fig. 5a) on a grid interval of 2.5 • × 2.5 • . The background dot colors in the smoothed map represent the stress data regimes with red denoting normal fault, blue as thrust fault, green as strike-slip fault and black as unknown regime. For the interpolation we only took the orientation pattern of the stress data into account. We limit our comparison with modeled lithospheric stress orientation to areas with enough data for interpolation. The new WSM2016 has a relatively good coverage in some regions that were not well covered in the previous version (Heidbach et al., 2010) such as Brazil, parts of North America, eastern Russia, and central Africa. We regard it as appropriate to compare the modeled stress orientation with the smoothed observational stress data and regard deviations of actual stress from smoothed stresses as a second-order pattern.

Angular misfit between WSM2016 and modeled lithospheric stress
In Fig. 5b we have superimposed our total modeled stress fields resulting from TM1 depicted by thin bars on top of the TM2 results as thick bars. There is relatively good agreement between the stress patterns and regimes at a longer wavelength; however, the smaller-scale contribution from the upper 300 km generate regional variations, which are mainly due to density contrasts in the lithosphere or below (which are almost isostatically compensated or cause lithosphere flexure) and topography (Zoback, 1992;Zoback and Mooney, 2003;Bird et al., 2006). Compared to the observed SH max patterns and regimes (Fig. 5a) we predict similar styles in regions such as eastern Africa and Tibet with normal faulting comparable to earlier works that considered the effect of the whole mantle including lithosphere and crustal models (Lithgow-Bertelloni and Guynn, 2004;Ghosh and Holt, 2012;Ghosh et al., 2013;Wang et al., 2015). We predict normal faulting mostly in regions above upwellings (mostly extensional regions) such as the Icelandic swell, eastern African Rift, or along divergent plate boundaries, while thrust faults are mainly found in compressional regions such as subduction zones and some other tectonically active regions in continents.
To further evaluate the influence of each thermal structure we performed a quantitative comparison between modeled and smoothed observed stress orientations. The angular misfit (Fig. 6) is the minimum angle between the modeled lithospheric stress orientation (Fig. 5b for TM1 and TM2) and smoothed observed stress orientation (Fig. 5a), which ranges from 0 to 90 • . Here, an angular misfit lower than 22.5 • is regarded as representing a good agreement between the modeled and observed stress orientations, with values above 67.5 • regarded as indicative of a poor fit. The general SSW to NNE stress orientation observed over the North American Plate is matched by our model predictions with both thermal structures TM1 and TM2. The angular misfit maps over North America obtained with both thermal structures show a poor fit over Yellowstone and the Rocky Mountains extending to the Great Plains (Ghosh et al., 2013). The observed localized NW to SE stress direction deviates (Fig. 5a) from the predicted long-wavelength stress pattern (Ghosh et al., 2013;Humphreys and Coblentz, 2007). Even though the thermal model TM2 includes high-density cratonic roots, as opposed to TM1, their respective results for the angular misfits show that the North American cratonic root has a limited influence on the stress field. The two density structures TM1 and TM2 yield mean values of 22.2 • (standard deviation, SD = 19.6 • ) and 22.9 • (SD = 20.7 • ), respectively. As the upper mantle thermal structure TM1 for the South American continent is not well constrained, due to lack of heat flow data, the predicted stress field in continental Brazil gives a relatively poor fit, with a mean misfit of 37.73 • (SD = 20.24 • ). However, TM2 does not perform much better resulting in a mean misfit of 33.79 • (SD = 21.9 • ). Both models fail to match the observed stress field in the Andes, where dominant localized N-S orientation is predicted, mainly as a results of the high topography and large crustal thickness compared to either Fig. 8c or 8d without the crustal contribution or stress field due to mantle below 300 km (Supplement Fig. S6a). In the African continent, predicted N-S stress orientations along the eastern African Rift from either model match the observed stress quite well with TM1 fitting observations much better compared to TM2 around the Ethiopia-Somalia-Yemen region but both fail over the Congo craton and the South African Plateau.
It has been suggested that the stress field in western Europe is influenced by the North Atlantic Ridge (NAR) push in the west and possibly by the far-field slab pull from the northwestern Pacific subduction zones, while in the south, the driving forces are induced by the convergence of the African and Eurasian plates, with Africa subducting under Eurasia in the Mediterranean (Zoback, 1992;Müller et al., 1992;Gölke, 1996;Heidbach and Höhne, 2007). Schiffer and Nielsen (2016) emphasize the importance of the anomalous mantle pressure underneath the North Atlantic lithosphere for generating the dominant first-order NW-SE stress pattern. In our study, due to mantle contribution > 300 km, we could match the NW-SE stress orientation almost perfectly, with the model using TM1 (Fig. 7a) showing small regional deviations, while the use of TM2 (Fig. 7d) results in larger deviations from this NW-SE pattern in some regions.
These regional pattern deviations between modeled and observed orientations are mainly induced by differences in the upper mantle density structures and topography (Heidbach and Höhne, 2007) (compared to Fig. 3a). The high density of heat flow data (Pollack et al., 1993;Artemieva, 2006) in continental western Europe (TM1) improves the fit to the observed stress field compared to the thermal structure based on S-wave velocity (TM2) yielding mean misfit values of 18.30 • (SD = 22.67 • ) and 19.9 • (SD = 22.64 • ), respectively. Similarly, the large amount of heat flow data in the Australian continent, improves the fit of the predicted intra-plate stress to the WSM2016 (Fig. 7c, mean = 23.07 • Figure 6. Angular misfit between the observed (WSM 2016) and total modeled stress directions with (a) TM1 and (b) TM2 upper mantle thermal and density structures. and SD = 19.4 • ) compared to TM2 (Fig. 7e, mean = 32.7 • and SD = 24.22 • ). It has been argued that the stress pattern in Australia is mainly driven by plate boundary forces (Reynolds et al., 2002), but based on the lithospheric and crustal structures used we show here that crustal and sublithospheric heterogeneities have a certain degree of influence.
In the Tibetan region, the collision of India and Eurasia leads to a complex crustal and lithospheric deformation (van Hinsbergen et al., 2011;Gaina et al., 2015) generating NE-SW compressional stresses. The SH max predictions with TM1 (Fig. 7c) fit better the stress pattern over the Tibetan Plateau with a mean misfit value of 28 • (SD = 23 • ) compared to TM2, where a predicted E-W direction results in a misfit of ∼ 50 • (Fig. 7f). Both models perform relatively poorly over parts of China, when compared to the observed stress field. The comparatively large angular misfit from the modeled stress field with only crustal and mantle contributions above the 300 km shows how much stress field is influenced by the mantle below (Supplement Fig. S6b). This is also supported by the fact that considering different thermal structures and/or correcting for the continental depletion does not seem to significantly improve the pattern of angular misfits (Supplement Fig. S7a-c).

Modeled "dynamic" topography
Following the above prediction of lithospheric stress field, we repeated the two simulations to compute the topography, but this time without crustal thickness variations (Fig. 8ab) to distinguish isostatic contributions from non-isostatic contributions. The corresponding stress magnitude and orientation from TM1 (Fig. 8c) and TM2 (Fig. 8d) without the crustal contribution are quite similar to the respective previous results that include the crustal contribution but show some regional differences, such as the N-S predicted stress orientation in the Andes in Fig. 4a-b compared to Fig. 8cd. Here, the resulting topographies with TM1 (Fig. 8a) and TM2 (Fig. 8b) show similar amplitudes due to the seafloor cooling and thickening along the ridges in the Atlantic, Indian and Pacific oceans, peaking above ∼ 1.5 km. With TM1, which explicitly contains subducted slabs, narrow, deep trenches are computed above subduction zones, such as in the northwestern Pacific and at the west coast of South America. Also the negative topography in the plate boundary south of Indonesia is reproduced well with the TM1 model reaching a value ∼ −1.8 km. Based on tomography (model TM2) the computed topographic lows are wider and less prominent.
Predicted topography with TM2 is higher in eastern Africa (2 to 2.5 km), and highly elevated regions are more extensive. Figure 8a with TM1 (based on seafloor age) shows relatively low topography amplitudes in the northwest of the Pacific Plate around Hawaii and towards the Mariana Trench compared to Fig. 8b with TM2 (based on the S-wave model SL2013sv) corresponding to a mean regional temperature difference of about ∼ 200 • C between TM1 and TM2 ( Fig. 1d-e). The "dynamic" topography with TM2 replicates nearly all island chains associated with hotspots in and around the African Plate, in the Pacific and along the Atlantic opening. In the North Atlantic, the positive topography (Icelandic swell) due to the Iceland plume-lithosphere interaction (Rogozhina et al., 2016;Schiffer and Nielsen, 2016) is more pronounced in Fig. 8b with TM2 based on the tomography of Schaeffer and Lebedev (2013). Here the heights exceed 2 km as compared to Fig. 8a with TM1 based on the ocean floor ages of Müller et al. (2008), showing values slightly below 2 km. The high isostatic topographic amplitudes along the mid-ocean ridges (MORs) as a result of high temperatures beneath these spreading centers where new seafloor is created are generally more pronounced in the TM2 model simulation than in the TM1 experiment. Despite the striking differences between topographic amplitudes in Fig. 8a and b along the MORs, the resulting modeled stress orientations (Fig. 8c-d) are very similar in these regions.
Also the large negative topography amplitude in cratons observed in dynamic topography with TM2 compared to TM1 does not readily translate into similarly large variations in the respective predicted SH max orientation (Fig. 8cd), showing that cratonic roots have less influence on the lithospheric stress field (Naliboff et al., 2012). Low temperatures as shown in the thermal model TM2 (Fig. 8b) translate into strong negative topographic anomalies, which are due to the conversion from seismic models to temperature and density, with the assumption that all seismic velocity anomalies are due to thermal variations only. This produces unrealistically strong density anomalies and hence, large negative topography in cratons (Forte and Perry, 2000), if correction due to the chemical depletion in the mantle lithosphere is not considered. Previous studies of cratonic mantle depletion in relation to density and temperature inferred from Swave models (for example, Cammarano et al., 2011) identified composition as the key dominant agent for the lowamplitude topography. They showed that a 100 K hotter mantle combined with lateral variations in composition resulted in a density of about 0.1 g cm −3 lower compared to models assuming pyrolitic composition. In contrast, Djomani et al. (2001) found that the depletion-related density drop in cratons is age-dependent and increases from 30 to 80 g cm −3 (i.e., 0.03 to 0.08 g cm −3 ) for Phanerozoic through Protozoic to Archean platforms. Here we aim at a qualitative first order analysis and therefore apply a density drop of 0.04 g cm −3 (modeled as an equivalent temperature increase of about 300 K) as correction in TM2 cratons. Also, following the realistic compositional correction in cratons by Cammarano et al. (2011) we adopt two additional thermal structures from different seismic tomography models SAW24B16 (Mégnin and Romanowicz, 2000) and S20RTS (Ritsema et al., 2007) with corrections applied to the depleted mantle based on the thermodynamic model PerpleX (http://www.perplex.ethz.ch, last access: 15 July 2017; Connolly, 2005) (See Supplement Fig. S8a-b) and compare with our results.

Comparing the modeled dynamic topography to the observation-based residual topography
Here, we compare our modeled dynamic topography to two independent observation-based residual topography fields (Hoggard et al., 2016;Steinberger, 2016). Residual topography gives a convenient way to constrain both isostatic and non-isostatic contributions to the modeled dynamic topography (Crough, 1978;Gurnis et al., 2000;Wheeler and White, 2000;Becker et al., 2014;Heidbach et al., 2016;Steinberger, 2016). This is done with the assumption that if topography is perfectly compensated isostatically within the upper mantle at depths within the range of 100-150 km, the integral of density with depth, as a function of crustal thickness and density to the Moho depth and of seafloor age will be the same everywhere for the chosen depth. The observation-based model by Hoggard et al. (2016) is derived from ocean seismic surveys (in situ) in oceanic regions and free-air gravity anomaly data in continents (Fig. 9a), while the residual topography model of Steinberger (2016) (Fig. 9b) is derived with the CRUST 1.0 model (Laske et al., 2013). These two models are comparable in most oceanic regions, but give large mismatches in continents. For example, the subducting plate under South America induces a negative anomaly in Fig. 9b but in the same region there is a positive anomaly in Fig. 9a due to the free-air gravity data used across continents. Hence, we perform a regional quantitative comparison for oceans and continents separately. To compare the modeled dynamic topography from TM1 and TM2 simulations ( Fig. 8a and b) to the observation-based fields ( Fig. 9a and b), we first remove the height due to ocean floor cooling. This is done by subtracting the height estimates from seafloor age  from the modeled dynamic topography, using the relation H topo = 3300 m × (1 − age 100 Ma ). Here we assume a half-space cooling for the sea floor with age. For a smooth transition of topographic height from ocean to continent and to avoid large jumps we nominally assume a 200 Ma lithosphere age for continents following the approach of Steinberger (2016). The resulting modeled dynamic topography fields ( Fig. 9c-d) corrected for the effect of the seafloor cooling with age with locations of active hotspot volcanism (Steinberger, 2000a) plotted as green dots show to which extent each of the models is able to predict the positive topographic amplitudes due to upwellings induced by plume heads pushing the lithospheric base.
A visual comparison of the two observation-based residual topography fields (Fig. 9a-b) with the modeled topography ( Fig. 9c-d) shows some features that are well reproduced such as the Pacific swell and the Hawaiian plume track, while the Canary Island plume, and the heights around southeastern Africa are much better reproduced by the TM2-based dynamic topography (Fig. 9d). Removing the height due to ocean floor age results in either zero or negative topographic amplitudes along MORs in the Atlantic and Indian oceans in the TM1-based dynamic topography (Fig. 9c), giving correlation of 0.323 and 0.198 (Table 1) (2016) with modeled dynamic topography using TM1 (c) and TM2 (d) upper mantle thermal density structures with the effect of seafloor cooling with age removed. (e) Similar modeled dynamic topography using TM2 upper mantle thermal density structures with constant temperature (300 K) added in cratons. Green dots with black circles around them show locations of major hotspots (Steinberger, 2000a). tively. This model uses the thermal density structure derived from the ocean floor age in the upper 300 km; hence, when this contribution is removed, only the lower mantle contribution remains. In contrast, the TM2 model still gives smallscale topography anomalies (Fig. 9d) due to density anomalies other than from the seafloor cooling at depths above (300 km), which are resolved by the seismic model used to derive TM2, thereby giving relatively higher correlation with S2016 and H2016 of 0.348 and 0.284 in oceans, respectively. To estimate the separate regional ratio between the modeled and observation-based residual topographies for continents and oceans, we assigned the continental mean value in continental areas to estimate the degree by degree ratio for oceans only (Fig. 10b) and vice versa for oceanic regions to estimate continents ratio (Fig. 10a).
In continents, the TM1 model (Fig. 9c) is similar to the residual models (Fig. 9b), exhibiting a correlation of 0.481 and a ratio of 0.98 (Fig. 10a) up to the spherical harmonic degree 30. The TM2 model gives similar ratio and correlation, but at degrees lower than 15 the TM2-induced modeled dynamic topography is about twice the amplitude of TM1 (Fig. 10a). Over the African continent with far less heat flow data used to derive TM1, this thermal density structure gives a large continental uplift up to about 2 km, similar to parts of Antarctica (Fig. 9c). In TM2 (Fig. 9d) this uplift is less extended, better resolving the negative topography of the Congo craton but reaching a height above 2 km over the East Table 1. Correlation between the modeled dynamic topography and the observation-based residual topography models (Steinberger, 2016;Hoggard et al., 2016) for continents and oceans.

Modeled topography
Steinberger (  African swell similar to S2016 (Fig. 9b). Many of the remaining continental regions, however, show large negative topographic magnitudes of −2 km and more, resulting from neglecting the compositional effects in cratons (e.g., Eurasia, Australia and North America). The wide range of variations shown in degree 1 to 2 ratio for continents ( Fig. 10a) is due to the strong contributions coming from the different cratonic structures in each thermal model. To further evaluate the impact of accounting for the correction due to chemical depletion in cratonic regions on the stress field and the dynamic topography, we have assumed an additional 300 K converted to a negative density as a compositional contribution in all cratons to the depth of 100 km for TM2; this is as opposed to the more realistic treatment of compositional effects as done for SAW24B16 and S20RTS, with the method from Cammarano et al. (2011). The modeled topography shows improvements in cratonic regions but there is almost no change in the resulting lithospheric stress field (Supplement Fig. S8c). The correlation with S2016 increases to 0.512 for TM2 (with an assumed 300 K compositional effect) in continents. SAW24B16 and S20RTS give much higher correlation 0.653 and 0.718 in continents (Table 1), which could be the result of a more realistic treatment of cratonic regions but also of using different seismic tomography models. Steinberger (2016), for example, used a similar simple procedure to convert seismic velocities from different tomography models to density and still obtained a rather high correlation of 0.64 in continents. The assumed compositional correction is not very large giving about a 100 m reduction in the cratonic negative anomaly (Fig. 9e) compared to the case without correction in continents (Fig. 9d). This in part supports the proposed treatment of the upper mantle thermal density structure with joint petrological and seismological constraints (Forte and Perry, 2000;Forte et al., 2010;Cammarano et al., 2011), which is outside the scope of our studies. The residual topography of Hoggard et al. (2016) shows positive amplitudes over the Eurasian craton due to the free-air gravity data used, while the other residual (Fig. 9b) and all modeled dynamic topography models give negative values, resulting in a low correlation with H2016 on continents for all models.

Conclusions
The aim of our study is to identify and quantify the influence of density anomalies and rheology in the crust and mantle on the present-day lithospheric stress field and dynamic topography. The focus is on anomalies and rheology above 300 km depth; therefore we use a number of different density structures, and nonlinear temperature and stress dependent rheology above 300 km. Our first upper mantle thermaldensity model (TM1) is based on heat flow data on continents (Artemieva, 2006) and seafloor age  in the oceans; while the upper mantle second thermal-density model TM2, and several alternative models considered, are based on seismic tomography (Schaeffer and Lebedev, 2013;Ritsema et al., 2011;Mégnin and Romanowicz, 2000). In contrast, only one density structure, based on the SMEAN (Becker and Boschi, 2002) tomography, and a radial viscosity structure (Steinberger and Calderwood, 2006) is used below 300 km depth. A key feature that distinguishes our work from previous studies is the use of a coupled code (Sobolev et al., 2009) that considers density heterogeneity in the entire mantle, along with a realistic lithosphere with free surface, such that lithosphere stresses are computed with a fully three-dimensional, rather than a thin-sheet approach.
Resulting lithosphere stresses are rather similar, both among the different models we consider, and to previously published results. They are also similar to the case where only the contribution from the mantle below 300 km is considered, showing that a larger portion of the contribution to the lithospheric stress field originates from mantle flow driven by density anomalies below 300 km depth (Steinberger et al., 2001;Lithgow-Bertelloni and Guynn, 2004). Only in some regions, particularly those with large and variable crustal thickness, such as Tibet, or the Altiplano, shallow contributions are dominant. A lower mantle stress contribution is dominated by very large-scale structures, with stress directions remaining similar over thousands of kilometers. It is related to very large scale mantle structures, which are well imaged by seismic tomography, causing overall similarity between our models and published ones. However, the modeled stress magnitudes coming from the mantle below 300 km or the total contributions (i.e., crust, lithosphere and the mantle below 300 km), are influenced by the respective density structures.
We compare computed directions of maximum compressive stress with the World Stress Map, and find a rather good overall agreement, confirming previous comparisons. However, regional comparison highlights those areas where the fit remains poor: these include the Colorado Plateau, the Altiplano, parts of Brazil, the Congo craton, and parts of China, highlighting regions on which future studies could focus. Computed stresses based on heat flow (Model TM1) compare more favorably to observations in those regions where heat flow coverage is good (e.g., western Europe), whereas the stresses computed from tomography (Model TM2) give a better fit for regions of poor heat flow coverage, such as South America.
In contrast to stress field, density anomalies above 300 km depth contribute dominantly to dynamic topography. Therefore, dynamic topography is more variable among the different models we consider and differs more strongly from published models. Dynamic topography also has a larger contribution at smaller scales. Some of these contributions can be related to subducted slabs or mantle plumes. Confirming previous results, we find that negative topography in cratons is too large, unless a correction for the depletion of cratonic lithosphere is considered. The best fit can be obtained, if the method of Cammarano et al. (2011) is used to convert seismic tomography models to temperature structures, taking chemical depletion in cratonic areas into account. The best agreement is found with residual topography on continents that considers crustal thickness variations based on CRUST1.0 (Laske et al., 2013) rather than deriving it from the gravity field. In order to fit either observable (stress or topography) attention has to be mostly paid to a detailed treatment of the Earth's parts (deeper or shallower) that give the largest contribution.
Code availability. The coupled global numerical code used to generate the results in this study builds on an in-house SLIM3D code (Popov and Sobolev, 2008) and a spectral code based on Hager and O'Connell (1981) and is available upon request from the main author. The coupling between the lithosphere and the mantle in our model allows for an implementation of realistic rheological parameters in both model domains. In SLIM3D, the stressand temperature-dependent rheology is implemented according to an additive strain rate decomposition into the viscous, elastic and plastic components: where G denotes the elastic shear modulus, Q = τ II is the plastic potential function,τ ij is the objective stress rate,γ denotes the plastic multiplier, τ ij = σ ij + P δ ij is the Cauchy stress deviator, P = −σ ii /3 is the pressure, τ II = (τ ij τ ij ) 1/2 is the effective deviatoric stress, and η eff is effective creep viscosity derived by combining the diffusion and dislocation creep mechanisms, as follows: The effective scalar creep strain rates are given by Kameyama et al. (1999): where the symbols A, E and V denote the experimentally prescribed pre-exponential factor, the activation energy and the activation volume, respectively, R denotes the gas constant, T is the temperature, n is the power law exponent, d is the grain size, and p is the grain size exponent, C H 2 0 is water content in ppm H/Si and r diff and r disl are the water content exponents. Along plate boundaries we account for the brittle deformation, with the yield stress defined according to the Drucker-Prager criterion based on the dynamic pressure: where c is the cohesion and µ is the coefficient of friction. Following Sobolev et al. (2009), andOsei Tutu et al. (2018) we use reduced friction coefficient values at the predefined plate boundaries (Bird, 2003) treated as narrow zones in the crustal and lithospheric layer in the depth range 0-80 km, and high friction coefficient of 0.6 in all lithospheric materials outside of the plate boundaries. Table A1. The upper mantle creep viscosity is calculated using olivine parameters from the axial compression experiments of Hirth and Kohlstedt (2004). Crustal rheology is taken from (Wilks, 1990). The rheological parameters used in this study with varying olivine water content of 100, 500, 1000 ppm H/10 6 Si in the weak asthenospheric mantle with dry lithosphere material. For more details regarding the formulation of the physical model and numerical implementation the reader is referred to Popov and Sobolev (2008