The Topographic Signature of Mantle Pressure Build-Up Beneath Subducting Plates: Insights From Spherical Subduction Models

Subduction zones are associated with spatially

• Positive dynamic pressure develops beneath subducting plates in subduction models within global domains • Slab dip is inversely correlated with the magnitude of this positive pressure • This pressure produces positive oceanic dynamic topography and tilts the plate upwards toward the trench Supporting Information: Supporting Information may be found in the online version of this article.
Irrespective of coupling and resolution, sinking slabs produce a dynamic draw-down of the surface (Gurnis et al., 1996;Mitrovica & Jarvis, 1985;Ricard et al., 1993;Rubey et al., 2017), as is reflected in RT (e.g., Faccenna & Becker, 2020;Husson, 2006) and sedimentary deposits (e.g., Gurnis, 1993;Heine et al., 2008;Mitrovica et al., 1989) on the OP.Other second order or more local topographic features do however require a decoupled and well-resolved subduction zone: for example, deep trenches, outer rises, and asymmetric across-slab DT (e.g., Crameri et al., 2017;Gérault et al., 2015;Zhong & Gurnis, 1992, 1994).Regarding the latter, subduction-driven topography outboard of the trench, that is, on the subducting plate (SP) side of the slab, has received less attention than the other topography components mentioned (cf.Husson et al., 2012).This is despite having the potential to broadly affect oceanic DT due to the prevalence of wide subduction systems rimming oceanic basins.Here, I focus on this topographic component using idealized spherical subduction models that contain decoupled and highly resolved slabs.
A static force balance, between slab buoyancy and dynamic pressure in the mantle, suggests positive pressure behind the slab (McKenzie, 1969;Stevenson & Turner, 1977) which should, in turn, produce positive SP dynamic topography: (Δ) cos() + Δ = 0 (1) Here, a dynamic pressure difference across the slab (ΔP) balances the slab-normal component of negative buoyancy, which is a function of the slab-mantle density contrast (Δρ), gravitational acceleration (g), slab thickness (l), and dip (θ).The dynamic pressure difference, ΔP, connects the (indirectly) observable dip to mantle dynamics, and is defined here as the negative near-slab pressure beneath a hypothetical OP minus the positive near-slab pressure beneath the SP.At depths shallower than ∼100 km from the base of the plates, the ΔP needed to support dips <90° is likely dominated by the negative pressure associated with corner flow in the mantle wedge beneath the OP (Stevenson & Turner, 1977;Tovish et al., 1978).At mid-upper mantle depths, plate and lateral slab motions dominate (Funiciello et al., 2003;Kincaid & Griffiths, 2003), and numerical models show that ΔP becomes more evenly partitioned between positive pressure beneath the SP and negative pressure beneath the OP (Holt et al., 2017).
In summary, Equation 1 suggests that a sub-slab positive pressure signal exists for dips <90° and that ΔP can be estimated observationally (from dip and oceanic plate age/density structure).This is true providing that: (a) ΔP does not come solely from negative pressure above the slab, and (b), the slab geometry is relatively close to steady state (i.e., slab-normal shear forces in the slab are significantly less than ΔP).Positive pressure beneath the SP should, in turn, push on the surface and produce positive DT.Such a positive SP pressure/topography signal emerges in the analytical upper mantle flow models of Holt and Royden (2020) and the analog subduction models of Husson et al. (2012).This occurs due to the confinement of mantle material behind the slab, and Husson et al. (2012) suggest this may be responsible for relatively shallow seafloor outboard of the Scotia, Tonga, and Kermadec trenches.However, extrapolating these results to Earth's global range of oceanic plates sizes, and hence assessing the global relevance of this source of topography, is challenging as Husson et al. (2012) consider relatively narrow plates/trenches (1,000 km) and, like Holt and Royden (2020), confine material to the upper mantle.
Observationally isolating a relatively subtle SP DT signal (below 300 m in most of our numerical single-slab models) is challenging.Such a signal could be overpowered by higher amplitude DT contributions and/or masked by uncertainties in the lithospheric isostatic correction needed to compute RT.I therefore turn to numerical models of subduction to isolate the nature of this sub-SP positive pressure build-up and the resultant DT, in order to assess its global relevance to oceanic topography and hence guide future observational studies.Time-dependent subduction models within a global domain models are relatively rare (Chamolly & Ribe, 2021;Chen et al., 2022;Morra et al., 2009Morra et al., , 2012;;Quevedo et al., 2013)

Modeling Overview
The ASPECT finite element code (version 2.2.0) is used to construct time-evolving subduction models within 3-D global domains (Bangerth et al., 2020(Bangerth et al., , 2021;;Heister et al., 2017;Kronbichler et al., 2012).I investigate five global models, with equivalent mechanical properties but variable plate geometries.Four of the models contain a single subduction zone with variably sized square SPs (2,500-11,000 km 2 : Figure 1), chosen to represent Earth's range of major oceanic plate sizes.For our fifth model, we consider a double subduction zone (two slabs dipping outwards from the central plate).
ASPECT was used to solve the conservation of momentum (Stokes equation) and mass (continuity equation) for an incompressible viscous fluid (Boussinesq approximation).Our models are purely compositional.That is, density and viscosity are dependent on compositional fields, which are advected through the domain using the field method (Bangerth et al., 2021).There are no external forces or velocities applied to the subduction system, and the upper and lower model boundaries are free slip.

Geometric and Material Properties
I use a simplified geometrical, mechanical, and rheological setup and neglect the overriding plate (cf.Chamolly & Ribe, 2021;Chen et al., 2022;Morra et al., 2012).Subduction is initiated by allowing a notch to extend to an initial depth of 150 km with a 250 km radius of curvature (Figure 1).
There are two compositional materials: the core and non-core (upper and lower) portions of the lithosphere.The 100 km thick lithosphere contains a 25 km thick core.Both lithosphere components have an equivalent density, 75 kg/m 3 greater than the background mantle (ρ mantle = 3,300 kg/m 3 ; η mantle = 2.5 × 10 20 Pa s), and the core is rheologically stronger (= 2500η mantle ) than the non-core material (= 100η mantle ).This is consistent with a strong core sandwiched between a brittle-yielding upper and ductile-yielding lower lithosphere (e.g., Karato & Wu, 1993).We also place stiff "core" material at the trailing edges of SPs to prevent subduction initiation here.
The non-core lithosphere undergoes plastic yielding above a prescribed yield stress (τ yield ) approximated by a Coulomb friction criterion: a is the friction coefficient (0.6), b the cohesion (60 MPa), and λ a pore fluid factor (0.1).This λ value results in significant weakening yet does not weaken the lithosphere completely.τ yield is used to calculate the plastic "viscosity": where    is the strain rate second invariant.The effective viscosity (η eff ) is then a harmonic mean of η yield and the isoviscous component.This produces averaged non-yielding lithosphere viscosities of ∼200 η mantle and yielding viscosities of ∼100 η mantle .Such moderate slab strengths are consistent with low elastic thicknesses at the trench (Billen & Gurnis, 2005), tomographically imaged slab morphologies indicating low viscosity deformation (e.g., Čížková et al., 2002), and plate driving force/slab dynamics considerations (e.g., Funiciello et al., 2008;Wu et al., 2008).The background mantle viscosity increases by factor 50 at 660 km depth.

Numerical Parameters
Adaptive mesh refinement (AMR) occurs for finite elements with large gradients in viscosity and/or non-zero compositions.This enables us to highly resolve the lithosphere while capturing mantle-scale flow.The AMR parameters produce a maximum refinement corresponding to ≈5.5 km radially thick finite elements (in the lithosphere) and a minimum refinement corresponding to ≈90 km elements (lower mantle).I have conducted tests to ensure our solver tolerances are sufficiently strict, and that upper mantle pressure fields are minimally affected by lower mantle resolution which, when low, leads to pressure artifacts in the lower mantle (Figure S1 in Supporting Information S1).

Model Analysis
To isolate links between pressure, topography, and subduction properties, the following quantities are extracted from the models: Slab Dip: Dip is computed at 250 km depth using the location of the slab upper surface at 50 km above and below this depth.
Dynamic Pressure: Dynamic pressure values in the upper mantle (i.e., the non-lithostatic component of full pressure), also at 250 ± 50 km depth, are extracted from beneath both the SP and a hypothetical overriding plate (OP).At each depth, P SP is extracted at 20 km from the lower surface of the slab (i.e., positive sub-slab pressure), P OP is extracted at 20 km from the upper surface (negative mantle wedge pressure), and ΔP = P OP − P SP is computed.Our final ΔP is depth-averaged: ΔP = 0.25 (ΔP 200 km + 2ΔP 250 km + ΔP 300 km ).
Dynamic Topography: The radial normal stress at the upper model boundary is used to estimate topography (topo = σ rr /Δρg) (e.g., Zhong & Gurnis, 1994), where Δρ is 2,300 kg/m 3 (water-loaded topography).The constant "isostatic" topography component, due to the weight of the uniformly thick SP, is subtracted to leave "dynamic" topography.I extract topography at 500 km from the trench, to be close to the slab yet avoid short-wavelength flexural topography, and also at 800 km to compute the topography gradient.

General Model Evolution
The four single slab subduction models evolve similarly.Firstly, the slabs sink through the upper mantle with 6-8 cm/yr vertical sinking rates (convergence is ≈50% rollback, 50% plate motion) and 60°-70° dips.After ∼10 Myrs, the slabs impinge on the upper-to-lower mantle boundary, are bent, and then sink slowly through the lower  1).This later phase is associated with greater dips (70°-80°), particularly for the largest plate (11,000 km 2 ) case.In the double subduction model, convergence is accommodated almost purely by rollback and dips are reduced (≈45°-60°).
I focus on model times between ≈7 and ≈20 Myr.This ensures the slab is deep enough to extract a dip yet avoids the latter-most stages during which subduction sometimes initiates at the plate edge.This could be avoided by placing a flat plate adjacent to the SP (e.g., Yamato et al., 2009); Here, I simply focus on the earlier phases.

Links Between Subduction Properties, Dynamic Pressure, and Topography
Upper mantle pressure: In all models, and during most timesteps analyzed, positive pressure develops behind (i.e., beneath the SP) and negative pressure develops ahead of retreating slabs (i.e., beneath the OP) at mid-upper mantle depths (Figure 2).On the SP-side, positive pressure build-up is mainly driven by SP motion (i.e., to drive flow in the opposite direction to the SP and hence conserve mass) and slab rollback (i.e., due to compression of the sub-slab mantle material).On the OP-side, negative pressure is also induced by rollback, with an additional contribution from "corner flow" in the mantle wedge.This wedge-directed decrease in pressure produces toroidal flow around the slab edge (Dvorkin et al., 1993;Funiciello et al., 2003;Kincaid & Griffiths, 2003) and a normal force on the slab that produces dips <90° (Stevenson & Turner, 1977;Tovish et al., 1978).The dynamic pressure is normalized to produce zero average pressure at this depth and the contours correspond to 10 MPa absolute dynamic pressure.
In single slab models, the magnitude of sub-slab positive pressure increases with plate size (Figures 2; 3a; Figure S2 in Supporting Information S1).This agrees with analytical expressions linking greater sub-slab pressure with larger plates/wider slabs (Royden & Holt, 2020).Of all the models, sub-slab pressure becomes near-zero only during the final stages of the small SP (2,500 km 2 ) model's evolution; otherwise, it is positive.In the double slab model, sub-slab pressure values are ≈65% greater than the single slab model with equivalent trench width.This is due to the confinement of mantle material between the two slabs (cf.Holt et al., 2017;Király et al., 2018) and a ∼50% increase in slab rollback relative to the single slab models (Figure S2 in Supporting Information S1).
Slab dip: I now consider the imprint of this sub-SP pressure build-up on dip and topography.Equation 1 shows how across-slab pressure difference (ΔP; Figure 3d) can be linked to dip (θ) via a static force balance.For both single and double slab systems, and variable plate sizes, the pressures and dips conform to this relationship to first order (Figure 3a), albeit with a 5-10 MPa offset (likely due to non-steady state effects, i.e., internal slab shearing).That is, as plate size increases or an extra slab is inserted, the pressure magnitude increases, which supports shallower dips (cf.Holt et al., 2017).During most timesteps, positive sub-slab pressure (P SP ) contributes at least half of the across-slab pressure difference, ΔP (= P OP − P SP ), that supports slab dip.Hence, models have large-variability in P SP , which decays to near-zero toward the SP edge (Figure 2).

Dynamic topography:
Positive pressure build up produces positive SP DT in the subduction models (Figure 3b, Figure S3 in Supporting Information S1).This was explored within Husson et al. (2012)'s analog models, and can also be seen in Crameri et al. (2017)'s numerical models.In our models, the magnitude of this topography reaches ∼450 m in the double slab case.In most of the single slab models, topography peaks at 100-300 m, and is maximum at the center of the subduction system where pressure build-up is greatest (cf. Figure 2, Figure S3 in Supporting Information S1).There are only a few timesteps during which SP topography is negative (Figure 3b), during the latter stages of the smaller plate model runs.Following the negative correlation between sub-SP pressure and dip (Figure 3a), there is a similar negative correlation between DT and dip (Figure 3b).This is because positive mantle pressure dictates both properties (i.e., by pushing slabs/Earth's surface).As for pressure, SP topography decreases from the sub-slab to the trailing edge causing an upwards tilt toward the trench (Figures 3c  and 3d).This trench-ward tilt reaches ≈27 m per 100 km (or 0.015°) and extends to a few-hundred km from the trench (where flexure becomes dominant).While drawdown/tilting of upper plates is well established (e.g., Mitrovica et al., 1989Mitrovica et al., , 1996)), and often of greater amplitude (e.g., ∼0.08° of Eocene tilt for Australia: DiCaprio et al. ( 2009

Modeling Limitations
Despite enabling us to extract first-order relationships, the models are highly simplified.Most notably, the overriding plate (OP) is neglected.This is due to resolution limitations but in line with other time-dependent subduction models within a global domain (Chamolly & Ribe, 2021;Chen et al., 2022;Morra et al., 2009Morra et al., , 2012)).Including an OP would reduce the rate of trench retreat relative to SP motion (Butterworth et al., 2012;Yamato et al., 2009) and shift the mantle wedge corner (and hence negative pressure associated with corner flow) to greater depth.This would result a ΔP signal that, while dictated by slab buoyancy (Equation 1), might be more dominated by negative P OP (at the expense of positive P SP ).
Other simplifications include mantle rheology and density structure.Including the power-law component of mantle viscosity, which we neglect, would reduce the suction force associated with corner flow (Billen & Hirth, 2005;Tovish et al., 1978) and shear stresses acting on the SP (e.g., Jadamec & Billen, 2010).This produces steeper dips (Billen & Hirth, 2005) and reduced trench retreat (e.g., Holt & Becker, 2017), and so could potentially offset enhanced rollback caused by neglecting the OP.
Flow driven by non-slab density anomalies, other subduction zones, deep slab material accumulated in the lower mantle, and the motions of other plates is also neglected.The pressure patterns associated with such far-field flow will affect DT around subduction zones-particularly by dragging down the surface above accumulated slab material (e.g., Pysklywec & Mitrovica, 1998;Rubey et al., 2017;Shephard et al., 2012) and subduction properties like dip and convergence rate (e.g., Chertova et al., 2018;Ficini et al., 2017;Holt & Royden, 2020;Husson, 2012).

Subducting Plate Topography on Earth
I now investigate whether such topographic signals exist at Earth subduction zones using Holdt et al. ( 2022)'s observationally derived RT.An update of Hoggard et al. (2016Hoggard et al. ( , 2017)), this data set uses active-source seismological data sets to correct for sedimentary and crustal loading, a plate cooling model to estimate and remove ocean lithosphere loading, and comprehensively covers the ocean floor (Figure S4a in Supporting Information S1).I compare this subduction zone topography (Figure 4a) with DT predictions from two recent convection models (Figure 4b and Figures S4c and S4d in Supporting Information S1) (Davies et al., 2022;Steinberger et al., 2019).Earth subduction zones are split into ∼ 220 km-long segments, with segments associated with ridges, plateaus, or continental material omitted (cf., Lallemand et al., 2005).For Holdt et al. (2022)'s data set, I also consider a case where segments without ≥5 nearby residual depth spot measurements are thrown out.This is to check the results are independent of the interpolation strategy used to generate a continuous RT field.For each data set, I assess whether observationally derived (RT) or modeled DT topography, and/or its gradient, is positive above these SP segments.
At a global scale, I do not find evidence for preferentially positive RT above SPs (Figure S5a in Supporting Information S1).The mean RT is weakly negative irrespective of whether all or the filtered subset of subduction segments are considered (−0.16 to −0.07 km).This is unchanged when long wavelength structure (spherical harmonic degrees <10) is removed from the global grid (Figure S6a in Supporting Information S1).The mean topography is more negative in the DT models: −0.56 to −0.31 km (Davies et al., 2022;Steinberger, 2016).Inclusion of the positive SP topography that we model would thus reduce this mismatch but this is speculative as the RT/DT distributions overlap significantly.This lack of robust evidence for an absolute SP DT signal is likely due to overprinting of this few-hundred meters signal by larger amplitude contributions.Within global geodynamic models, the gradual sinking of relic slabs from previous episodes of subduction -particularly those located in the lower mantle -have been shown to draw-down subduction zones by heights on the order of a km (e.g., Conrad & Husson, 2009;Hager, 1984;Ricard et al., 1993).Our slabs are mainly contained within to the upper mantle (Figure 2) and so we do not include such effects.Also, calculating RT above old ocean floor relies on semi-empirical corrections for ocean plate loading that may also mask the modeled topographies (e.g., Holdt et al., 2022;Parsons & Sclater, 1977;Stein & Stein, 1992).The modeled topography should also be associated with a gravity signal: For 200 m of positive topography, approximated as a half-ellipsoid, simple calculations predict a maximum gravity signal of 20-35 mGal at the center of the anomaly (Figure S7a in Supporting Information S1).Such a signal is broadly compatible with positive free-air gravity observed at some (e.g., Central and HOLT 10.1029/2022GL100330 8 of 11 South America, Aleutians, Kuril), but not all (notably Japan and Sumatra), major subduction zones.In both cases, further work is needed to isolate this potential gravity signal from that due to excess slab mass.
I now investigate the modeled positive topography gradient: An upwards tilt of the SP toward the trench.Figure 4 shows the distribution of Δ x topo for Earth subduction zones, extracted from each topography map.Observational estimates, from the RT map, are more positive than the topographic gradients predicted by mantle flow models.Subduction zones within Holdt et al. (2022)'s data set exhibit a mean upwards tilt of 0.25 m/km, and this increases to 0.37 m/km when only subduction segments near ≥5 observations are included (Figure 4a).In contrast, the DT models exhibit negative/downwards average tilts of −0.68 to −0.56 m/km (Figure 4b).The average tilt in the RT data set reduces when topography is extracted further from the trench (Figure S8 in Supporting Information S1), but Δ x topo in always shifted to more positive values relative to the DT simulations.This holds when comparing Steinberger et al. (2019)'s DT to the RT generated in that paper using a global, Crust1.0-basedcrustal correction (Laske et al., 2013).Steinberger's mean tilt is ∼1 m/km more positive within the RT models, irrespective of whether ocean plate loading is removed from the RT/DT maps or not (Figure S9 in Supporting Information S1).I propose this mismatch is evidence for an upwards tilt of SPs that, despite large RT scatter (due to the dominance of other topographic sources), is not reflected in flow models without decoupled and/or high resolution slabs.The negative tilt in mantle flow models is due to the dynamic drawdown driven by dense slab anomalies, in  (Davies et al., 2022;Steinberger et al., 2019).Topography gradient (Δ x topo) is extracted 400 km outboard of the trench, and positive values represent upwards tilt toward the trench.Horizontal bars show the arithmetic means plus/minus the standard deviation.The gray bar represents the range of Δ x topo values in our subduction models (Figure 3c), scaled up according to the difference between the modeled ΔP range (∼30 MPa) and the "observed" range calculated from Equation 1 (∼75 MPa).The Steinberger (2016) topographies are expanded up to spherical harmonic degree 31 (and is their case without lateral viscosity variations) and Holdt et al. (2022) up to degree 41.Davies et al. (2022) is their case with both upper mantle structure and temperature/pressure-dependent viscosity.In the "filtered" version of Holdt et al. (2022), subduction segments are counted only if they are within 150 km distance of 5 or more residual depth measurements.contrast to the dominance of sub-SP pressure in our idealized subduction models.This drawdown is accentuated in older flow computations (Conrad & Husson, 2009;Spasojevic & Gurnis, 2012) that are more dominated by the larger-scale density anomalies present in older tomographic models (Figure S10 in Supporting Information S1).

Conclusions
I have used time-dependent, spherical subduction models to elucidate links between subduction properties, mantle pressure, and SP dynamic topography.The magnitude of upper mantle dynamic pressure, which is positive behind and negative ahead of slabs, increases with plate size, slab length, and trench velocities.This agrees with previous analytical and numerical subduction models (e.g., Royden & Husson, 2006;Stegman et al., 2006) but extends this framework globally.As across-slab pressure difference (ΔP) increases, for example, with increasing plate size, dip decreases due to an increased mantle-wedge directed pressure force.This is approximated by a force balance that shows how ΔP can be estimated from dip (cf.Holt & Royden, 2020).
The pressure distribution produced by the confinement of mantle material behind a slab -that is, positive behind the slab and decaying with slab-normal distance -produces positive SP DT and causes the plate to tilt upwards toward the trench (at distances beyond a few-hundred km from the trench; closer, flexure dominates).For a Pacific Plate sized plate (∼11,000 km 2 ), positive DT of 100-400 m is produced close to the slab.This is reduced to ≤300 m for a Nazca sized plate (∼5,000 km 2 ).In both cases, a trench-upwards tilt, of up to 0.25 m/km, is recorded.
SPs within Holdt et al. (2022)'s RT data set exhibit tilts that are, in general, more positive than DT extracted from recent convection models.I interpret this as evidence that DT simulations, which do not contain decoupled slabs, are missing a component of Earth's DT related to sub-slab and sub-SP pressure build-up.Unlike tilt, I do not find evidence for positive absolute topography above SPs at a global scale.This is unsurprising as sinking slab remnants, and/or RT uncertainties, could offset this relatively subtle (few-hundred meters) signal.Moving forward, geographically accurate global subduction models are needed to confirm and parse out the relative contribution of this source of oceanic DT.This study has benefited from the insights and help of Leigh Royden, Sam Goldberg, and Thorsten Becker.The computations used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation (NSF) Grant ACI-15485x62.This study was partially supported by NSF Grant EAR-2147997.I also thank the Computational Infrastructure for Geodynamics (geodynamics.org),which is funded by the NSF under awards EAR-0949446 and EAR-1550901, for supporting the development of ASPECT.

Figure 1 .
Figure 1.All initial model geometries (a) and temporal evolution of the single slab, 8,000 km 2 subducting plates size, model (b) and the double slab model (c).Viscosity is plotted along a vertical cross section through the center of the subduction system and overlain by velocity vectors.

Figure 2 .
Figure2.Dynamic pressure evolution within two single slab subduction models (2,500 and 8,000 km 2 subducting plates sizes) and the double slab model.Snapshots show the dynamic pressure field and mantle velocity vectors at 330 km depth, the surface location of the SP (dashed white line), and insets of the plate and slab shape (shaded by velocity magnitude).The dynamic pressure is normalized to produce zero average pressure at this depth and the contours correspond to 10 MPa absolute dynamic pressure.

Figure 4 .
Figure 4.The distribution of subducting plate topography gradient within: (a) observationally based residual topography (Holdt et al., 2022); (b) DT models from mantle flow computations(Davies et al., 2022;Steinberger et al., 2019).Topography gradient (Δ x topo) is extracted 400 km outboard of the trench, and positive values represent upwards tilt toward the trench.Horizontal bars show the arithmetic means plus/minus the standard deviation.The gray bar represents the range of Δ x topo values in our subduction models (Figure3c), scaled up according to the difference between the modeled ΔP range (∼30 MPa) and the "observed" range calculated from Equation 1 (∼75 MPa).TheSteinberger (2016) topographies are expanded up to spherical harmonic degree 31 (and is their case without lateral viscosity variations) andHoldt et al. (2022) up to degree 41.Davies et al. (2022) is their case with both upper mantle structure and temperature/pressure-dependent viscosity.In the "filtered" version ofHoldt et al. (2022), subduction segments are counted only if they are within 150 km distance of 5 or more residual depth measurements.