The role of sediment subduction and buoyancy on subduction dynamics and geometry

18 Subducted sediments are thought to lubricate the subduction interface and promote faster 19 plate speeds. However, global observations are not clear-cut on the relationship between 20 the amount of sediments and plate motion. Sediments are also thought to influence slab 21 dip, but variations in subduction geometry depend on multiple factors. Here we use 2D 22 thermomechanical models to explore how sediments can influence subduction dynam23 ics and geometry. We find that thick sediments can lead to slower subduction due to an 24 increase of the megathrust shear stress as the accretionary wedge gets wider, and a de25 crease in slab pull as buoyant sediments are subducted. Our results also show that larger 26 slab buoyancy and megathrust stress due to thick sediments increase the slab bending 27 radius. This offers a new perspective on the role of sediments, suggesting that sediment 28 buoyancy and wedge geometry also play an important role on large-scale subduction dy29 namics. 30 Plain Language Summary 31 At subduction zones, an oceanic plate dives into the mantle below another plate. 32 The downgoing plate is usually covered by sediments. These sediments can be carried 33 down to depth along the interface and/or scraped off the top of the downgoing plate and 34 appended to the edge of the upper plate, forming an accretionary wedge. Sediments sub35 ducted to depth act as a lubricant, influencing the shear resistance of the interface, and 36 in turn, downgoing plate speed. However, natural data show that slow subduction can 37 be associated with thick sediments. Sediments are also thought to affect the dip angle 38 of the downgoing plate, but subduction geometry is also influenced by other factors. We 39 conducted a numerical modeling study to understand the effect of sediment thickness 40 and density on the downgoing plate speed and dip. We observe that thick sediments on 41 the downgoing plate lead to a slower subduction and a shallower dip, due to the decrease 42 in slab pull and increase of stress along the contact interface associated to a bigger ac43 cretionary wedge. Our findings suggest that the effect of sediments might be not lim44 ited to the lubrication of the contact interface, but buoyancy and accretionary wedge size 45 also play a role. 46

X -2 : Text S1. Numerical method We use the two-dimensional, continuum, visco-elasto-plastic, seismo-thermo-mechanical version (van Dinther et al., 2013) of the code I2ELVIS (Gerya & Yuen, 2007), which uses an implicit, conservative finite difference scheme on a fully staggered Eulerian grid in combination with a Lagrangian marker-in-cell technique.
The code solves for the conservation of mass, momentum and energy with a visco-elastoplastic rheology. Lagrangian markers advect physical properties (e.g., viscosity, stress, plastic strain, and temperature) according to the velocity field interpolated from the Eulerian grid (Gerya & Yuen, 2007).
The continuity (Eq. 1), and momentum (Eq. 2 -3) equations are solved to obtain the horizontal and vertical components of velocity, v x and v z , and pressure P (defined as mean stress): where ρ is density, σ ij are the deviatoric stress tensor components, and g = 9.8 m/s 2 is the gravitational acceleration.
The continuity equation (Eq. 1) assumes an incompressible medium, i.e., Poisson's ratio where C p is the isobaric heat capacity, q x and q z are the horizontal and vertical heat flux, H a is the internal heat generation due to adiabatic (de)compression, H s is the shear heating due to anelastic deformation, and H r is the radioactive heat production. q x , q z , H s , and H a are defined as follows: where k is thermal conductivity,ε ij,vp is the visco-plastic component of the deviatoric strain rate tensor, and α is the thermal expansion coefficient. H r is constant for each rock type (Table S1). A constitutive equation (Eq. 9) relates visco-elasto-plastic deviatoric stresses and strain ratesε ij by applying linear elasticity and non-Newtonian viscosity (Gerya & Yuen, 2007): where η is the effective viscosity, G is the shear modulus,

Dσ ij
Dt is the objective co-rotational time derivative solved using a time explicit scheme, χ is a plastic multiplier connecting plastic strain rates and stresses, σ II = σ 2 xx + σ 2 xz is the second invariant of the stress tensor, and σ yield is the local plastic yield strength.

X -4 :
Plastic failure is simulated with the Drucker-Prager plastic yielding criterion. At each Lagrangian marker, yielding occurs when the second invariant of the deviatoric stress tensor σ II reaches the local plastic yield strength σ yield defined as where C is cohesion, µ s is a constant static friction coefficient, P is pressure, and λ is the pore-fluid pressure factor (P f luid /P ). and wet olivine rheology (Ranalli, 1995). We impose a constant velocity of 7.5 cm/yr within a small region of the subducting plate (see Fig. 1) until 300 km of the slab is subducted into the mantle. This corresponds to 4 Myr of subduction. After this kinematically prescribed phase, the pushing velocity is removed, and subduction is selfdriven.
Text S3. Estimation of the percentage of accreted and subducted sediments.
To estimate how sediments are partitioned between the accretionary wedge and the subduction channel, we use compositional maps. At each time step, we contour all sediment rock types (i.e., incoming plate and pre-existing wedge) starting from the trench (Fig. S1a). We then compute the area of this contour polygon, which is a measure of the total amount of sediments. We then discriminate between sediments that are above and below the continental Moho (∼40 km depth). The amount of subducted sediments d ss is computed as the area of sediments below the continental Moho (Fig. S1b). To estimate X -6 : the amount of accreted sediments, we consider the area of sediments above the continental Moho and we identify the pre-existing wedge sediments (Fig. S1c). The amount of accreted sediments is computed as the difference between the area of sediments above the continental Moho and the area of the pre-existing wedge sediments (Fig S1d). The percentages of accreted and subducted sediments are computed considering the area of the total amount of sediments.  Table S1. Rheological parameters. ρ 0 is the reference density, E a is the activation energy, V a is the activation volume, n is the stress exponent, η 0 is the reference viscosity, H r is the radiogenic heat production, G is the shear modulus, µ s is the static friction coefficient and λ is the pore-fluid pressure factor (P shear stress F sl is computed as the sum of σ || within a 3 km-wide polygon that extends from the trench to the 450°C isotherm (red rectangle) times the length of the polygon.
For simplicity, the length of the polygon is approximated to a straight line.
June 5, 2021, 8:10pm X -14 : Figure S5. a) Average shear stress of the plate interface as a function of sediment thickness d sed and density ρ sed . The average shear stress is the average of the second invariant of the deviatoric stress tensor σ || in a 3 km wide region that extends from the trench to the 450°C isotherm (see also Fig. S4). b) Relationship between integrated shear stress of the megathrust F sl and interface length for models with different sediment thickness d sed and density ρ sed .