The effect of low-viscosity sediments on the dynamics and accretionary style of subduction margins

Observations of sediments at subduction margins appear to divide them into two classes: accretionary and erosive. Accretionary margins are dominated by accretion of thick piles of sediments (>1 km) from the subducting plate, while tectonic erosion is favored in regions with little or no sedimentary cover (<1 km). The consequences of the two styles of margins for subduction dynamics remain poorly resolved. In this study, we used 2-D numerical simulations of subduction to investigate how low-viscosity sediments influence subduction dynamics and margin type through plate coupling. We vary the thickness and viscosity of the sediment layer entering subduction, the thickness of the upper plate, and the driving velocity of the subducting plate (i.e., kinematic boundary conditions). Diagnostic parameters are extracted automatically from numerical simulations to analyze the dynamics and differentiate between modes of subduction margin. We identify three margin types based on the extent of viscous coupling in the sediment layer at the subduction interface: (a) tectonic coupling margin, (b) low-angle accretionary wedge margin, and (c) high-angle accretionary wedge margin. In the tectonic coupling case – analogous to an erosive margin – high-viscosity or thin-layer sediments increase coupling at the interface. On the other hand, when the viscous coupling is reduced, sediments are scrapped off the subducting slab to form an accretionary wedge. Models that develop tectonic coupling margins show small radii of curvature, slow convergence rates, and thin subduction interfaces, while models with accretionary margins show large radii of curvature, faster convergence rates, and dynamic accretionary wedges. These diagnostic parameters are then linked with observations of present-day subduction zones.

1 Introduction

Sediment subduction at convergent plate boundaries has long been recognized to play an important role in the dynamics of our planet as it can provide direct feedbacks between plate tectonics, climate, and life. Quantifying the sediments mass flux through subduction zones is important for understanding (i) generation of arc magmas and the problem of petrogenesis of continental crust (e.g., Plank and Langmuir, 1998; Kelemen and Manning, 2015), (ii) whether large volumes of existing continental crust are ever recycled back into the mantle over long periods of geologic time, and (iii) cycling of volatiles from Earth's crust and atmosphere to its deep interior (e.g., Hawkesworth et al., 1997; Plank and Langmuir, 1998; Dasgupta and Hirschmann, 2010). Regarding the latter, carbon and water global cycles in particular depend greatly on the amount of subducted sediment (e.g., Plank and Manning, 2019; Dutkiewicz et al., 2018; Merdith et al., 2019), which in turn have important implications for climate stability (Kasting, 1989), biogeochemical cycles (Husson and Peters, 2017), and the rheology of the mantle (Hirth and Kohlstedt, 1996).

Sediments are also fundamental to the dynamics of subduction zones and plate tectonics. The lubricating effect of sediments at the plate interface (referred to here as the domain between the subducting and upper plates, as defined in Agard et al., 2018) was recognized more than 30 years ago (e.g., Shreve and Cloos, 1986) and is critical for the mechanism of plate tectonics (e.g., Tackley, 2000; Bercovici, 2003; Bercovici and Ricard, 2014; Sobolev and Brown, 2019). On a regional scale, sediments influence patterns of deformation by controlling the morphologies of subduction interfaces, accretionary prisms, and forearc basins (Clift and Vannucchi, 2004; von Huene and Scholl, 1991b; Melnick and Echtler, 2006; Simpson, 2010). Moreover, sediments occupying the shallow seismogenic subduction interface, for example, appear to influence seismic coupling and the frequency of megathrust earthquakes (e.g., Moore and Saffer, 2001; van Rijsingen et al., 2018; Heuret et al., 2012; Brizzi et al., 2020; Bangs et al., 2020).

However, the relative role of sediments in subduction dynamics and evolution remains unclear. The force balance during subduction includes the potential energy change of the negative buoyancy of the subducting slab, viscous dissipation in the mantle, bending of the lithosphere, and slab–upper plate interface (e.g., Conrad and Hager, 1999). Considerable effort in subduction dynamics studies focused on quantifying dissipation due to slab bending (e.g., Conrad and Hager, 1999; Becker et al., 1999; Capitanio and Morra, 2012; Garel et al., 2014a). That is because slabs were assumed to be strong (3000–5000 times stiffer than the mantle) and sediments weak, implying a low coupling degree (low shear stresses) at the interface (Conrad and Hager, 1999; Duarte et al., 2015; Billen and Hirth, 2007). Laboratory and numerical models instead suggest that slabs may be weaker (only 150–500 times stiffer than the mantle) (Funiciello et al., 2008; Moresi and Gurnis, 1996; Zhong and Davies, 1999), implying a more prominent role for plate interface dissipation than previously thought. Recent simplified energy balance calculations by Behr and Becker (2018) also suggest that sediment subduction may modulate plate speeds, with sediment-lubricated plates subducting faster than slabs with metabasaltic (i.e., exposed mantle rocks) interfaces.

The lubricating effect of sediments has largely been considered an implicit assumption in previous large-scale subduction dynamics studies. The subduction interface in numerical models is typically implemented by imposing either a fixed interface layer (i.e., subduction channel or weak fault) or a layer at the top of the subducting plate (i.e., weak crust) that is advected with the flow and continuously entrained into the decoupling region. In this way, the weakening effect of sediments or any other deformation-localizing processes such as damage, grain size reduction, and fabric development, is parameterized by setting a low viscosity or low friction coefficient of the interface material (Gerya et al., 2002). This is a necessary model component for stable asymmetric subduction (Petersen et al., 2017; Crameri and Tackley, 2015; Gerya, 2009; Sandiford and Moresi, 2019) and has become an increasingly common strategy in the last decade (Babeyko and Sobolev, 2008; Capitanio et al., 2010; Magni et al., 2012; Chertova et al., 2012; Cizkova and Bina, 2013; Garel et al., 2014b; Pusok and Kaus, 2015; Agrusta et al., 2017; Pusok and Stegman, 2019).

Indirect observations suggest that subduction interfaces are not discrete isosurfaces separating two plates but are rather exemplified by melange zones (Shreve and Cloos, 1986; Vannucchi et al., 2008; Agard et al., 2018) that incorporate material from the subducting plate, the accretionary prism, and the upper plate (Menant et al., 2020; Angiboust et al., 2021). Moreover, when sediments are considered at the trench, convergent margins appear to fall into one of two classes: accretionary and erosive (i.e., Clift and Vannucchi, 2004, Supplement, Fig. S1). Accretionary margins develop from an accumulation of material from the subducting plate being transferred onto the overriding plate, either by frontal off-scraping at the trench axis or by underplating of the forearc wedge above the decollément at greater depths (Angiboust et al., 2021). Erosive margins develop from a strong coupling between overriding and subducting plates that results in erosion of the underside of the upper plate, as indicated by margin truncation and forearc subsidence (von Huene and Scholl, 1991a; Clift and Vannucchi, 2004; Straub et al., 2020). Accretionary margins are dominated by accretion of thick piles of sediments (>1 km), while tectonic erosion is favored when the sedimentary cover is <1 km at the trench.

Both accretionary and tectonic erosion margins are widely distributed. Clift and Vannucchi (2004) classified the global subduction zones in the two categories based on their dominant mode in the last 10 Myr and found that 43 % of global subduction margin is accretionary and 57 % is tectonic erosion. The implications of this equally distributed duality in margin type to the global subduction system has not been investigated. Convergence rates for the two types of margins seem to correlate well with sediment thickness, the taper angle, and radius of curvature (i.e., Supplement Fig. S1, replotted data from Clift and Vannucchi, 2004, and Wu et al., 2008). De Franco et al. (2008) also observe a correlation between the margin type and upper plate strain (i.e., proxy for back-arc extension). Lamb and Davis (2003) went further to suggest that the type of margin can affect mountain building, with tectonic erosion producing a higher degree of coupling between the subducting and upper plate. They argue that changes from a sediment-rich to sediment-starved subduction regime during Cenozoic climatic cooling may have been responsible for the rise of the Andean mountain belt.

It is clear that the consequences of the two styles of margins for large-scale subduction dynamics remain poorly resolved and have not been explored extensively with numerical models. Accretionary margins have been investigated in more detail, following the theory for critical Coulomb wedges by Dahlen (1984) and Dahlen et al. (1984) (i.e., in analogue models – Lallemand et al., 1994; Gutscher et al., 1998, and numerical models – Beaumont et al., 1999; Selzer et al., 2008; Ruh, 2017; Menant et al., 2020). However, accretionary margins have generally not been studied in the same framework as erosive-style margins. Thus, we identify a number of outstanding questions regarding the influx of sediments to trenches and the style of margin that could be addressed with numerical models. Why do some margins accrete sediments while others do not? What is the feedback between sediment fluxes and subduction dynamics? How much sediment material gets subducted into the mantle?

In this study, we run systematic 2-D numerical simulations of ocean–ocean subduction to investigate how low-viscosity sediments influence subduction dynamics and the plate coupling. We aim to understand what causes convergent margins to either accrete material delivered by the subducting plate or, alternatively, to subduct the trench sediment pile and even erode the basement of the overriding plate. The purpose of the present work is not to model the dynamics of accretionary or erosive margins in detail but rather, by carrying out numerical experiments on the effect of sediments in geometrically simple configurations of subduction, to further understand the occurrence and evolution of each style of margin.

We begin our investigation with a discussion of the numerical setup and diagnostics. We then extract automatic diagnostics from numerical results to evaluate end-member regimes and compare to available observations in natural subduction zones. In particular, we considered a range of dependent and independent variables from statistical analyses of present-day subduction zones (Sect. 2; Clift and Vannucchi, 2004; Lallemand et al., 2005; Wu et al., 2008; De Franco et al., 2008; Heuret et al., 2012) to constrain and validate the results of numerical models. This study intends to consolidate insights from numerical models of subduction with an integrated set of global observations.

2 Methods

Numerical models presented below are purely mechanical. The lithosphere and mantle are assumed to behave as a slow-creeping fluid deforming over a timescale of millions of years (Turcotte and Schubert, 2014). Conservation of mass and conservation of momentum for the fluid motion, assuming incompressibility and neglecting thermal diffusion, are given by

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{v}=\mathrm{0},\text{(2)}& {\displaystyle }-\mathrm{\nabla }P+\mathrm{\nabla }\cdot \mathit{\tau }+\mathit{\rho }\mathit{g}=\mathit{0},\end{array}$$

where τ is deviatoric stress tensor, P is pressure, ρ density, g is the gravity vector, and v is velocity. We use a variable viscosity constitutive relationship ${\mathit{\tau }}_{ij}=\mathrm{2}\mathit{\eta }{\dot{\mathit{\epsilon }}}_{ij}$, where η is the Newtonian viscosity that is constant for each material phase, ${\dot{\mathit{\epsilon }}}_{ij}=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}\right)$ is the deviatoric strain rate tensor, and i and j represent spatial directions following the Einstein summation convention.

The above equations are solved using the parallel 3-D finite-difference code LaMEM (Lithosphere and Mantle Evolution Model) (Pusok and Kaus, 2015; Kaus et al., 2016). We use a pseudo 2-D Cartesian domain in an approach similar to Pusok and Stegman (2019, 2020), meaning we consider infinite extension in the y direction. A Lagrangian marker-in-cell method (Harlow and Welch, 1965; Gerya, 2009) is used for accurately tracking distinct material domains (Pusok et al., 2017) as they undergo extensive deformation due to creeping flow. We also employ an internal free surface using the “sticky air” approach (Schmeling et al., 2008; Crameri et al., 2012), with a free surface stabilization algorithm (Kaus et al., 2010) that allows for the development of topography.

2.1 Model setup

We performed 2-D numerical simulations of ocean–ocean subduction (Table 1, Fig. 1) to investigate the role of sediments in slab dynamics and topography. The model domain is 6000 km wide and 1500 km deep. The computational domain has a variable grid spacing, with higher resolution in the upper mantle asthenosphere and close to the subduction trench (minimum and maximum grid spacings: $\mathrm{\Delta }x\in [\mathrm{1.73},\mathrm{15.62}]$ km, $\mathrm{\Delta }z\in [\mathrm{2},\mathrm{16}]$ km). Free-slip boundary conditions are imposed on all boundaries and a 60 km layer of “sticky air” on top of the plates. Previous studies have shown that using a free surface (i.e., sticky-air method in this case) instead of a free-slip top boundary dramatically changes subduction style (Kaus et al., 2010; Crameri and Tackley, 2015). Here, the rock–sticky air interface represents an internal free surface formulation, from which topography is calculated.

The initial model setup and material parameters are similar to the ones used in Pusok and Stegman (2019, 2020). These models were used to investigate the dynamics of same-dip double subduction systems and to explain dynamical processes leading to the fast convergence history between India and Eurasia in the Cretaceous. Here, we model a single subduction setting (Fig. 1). The model consists of an oceanic plate subducting beneath another oceanic upper plate (i.e., ocean–ocean subduction). The length of both plates is 2500 km and they are not attached to the side walls, mimicking mid-ocean ridges at the trailing edge of the plate away from the trench. Additional experiments with a fixed upper plate to the wall are presented in the Supplement (Fig. S8).

The subducting plate lithosphere has a thickness of 80 km with a 20 km thick core and 15 km combined weak crust and sediments. Behr and Becker (2018) estimate that metasediments and serpentinites can be more than 2 orders magnitude weaker than the reference asthenospheric mantle viscosity, while metabasalts are predicted to be of the same order of magnitude viscosity at temperatures 600–800 ^∘C. In our model, we assume that the weak crust is formed by a layer of highly fractured metabasalts (dark green lithology in Fig. 1a, b) overlain by a layer of metasediments (light green lithology in Fig. 1a, b). Thus, the viscosity structure for the combined weak crust and sediments represents a parameterization of the strength weakening with depth due to hydration and weak sediment cover. The combined initial thickness of sediments and crust remains the same in all simulations.

Material parameters are the ones used in Enns et al. (2005) and Pusok and Stegman (2019); the upper mantle asthenosphere has a reference density (ρ₀=3300 kg m⁻³) and viscosity (${\mathit{\eta }}_{\mathrm{0}}=\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{20}}$ Pa s), and the plates are 85 kg m⁻³ denser and have a variable viscosity structure (lithospheric mantle: 500×η₀, strong core: 5000×η₀ and weak crust: η₀). The transition to the lower mantle is marked by a viscosity jump of 50 in the reference models, consistent with previous estimates of the viscosity jump in the Earth's mantle (Quinteros et al., 2010; Rudolph et al., 2015). For initial conditions, we impose a slab radius of curvature of 150 km and a subduction depth of 200 km, which is enough to initiate subduction.

2.2 Input and diagnostic parameters

Previous studies investigated the role of upper plate and subduction plate parameters, such as thickness and strength (i.e., Holt et al., 2015; Brizzi et al., 2020). Here, we focus on factors acting directly on the subduction interface. In particular, we use a result from Currie et al. (2007) and Cizkova and Bina (2019), which found that the effects of sediment buoyancy and viscous entrainment by the subducting plate are the main factors controlling the behavior of slab and subducted sediments. Thus, we vary the properties of the sediment layer (viscosity, thickness) and thickness of the upper plate as explained below and in Table 1. All other parameters are kept the same among simulations. By not changing the density of the sediments or the slab geometry, the magnitude of initial slab-pull force is the same among simulations. We also extract automatic diagnostics from our numerical models that can be compared to parameters available for the natural subduction systems. Parameters discussed in this section are listed in Table 2.

Input parameters. For the Newtonian rheology used, the input parameters for each phase (viscosity, thickness, density) are categorized into three groups: (1) subducting plate (slab) parameters, (2) interface (sediments) parameters, and (3) upper plate parameters (Table 2). The following parameters were varied: thickness and viscosity of sediments ($h_{\text{sed}}=[\mathrm{5},\mathrm{10}]$ km, ${\mathit{\eta }}_{\text{sed}}=[\mathrm{0.01}{\mathit{\eta }}_{\mathrm{0}},\mathrm{0.1}{\mathit{\eta }}_{\mathrm{0}},{\mathit{\eta }}_{\mathrm{0}}]$ Pa s) and thickness of the upper plate (h_UP = [50, 80, 100, 150] km). The sediment viscosity is consistent with estimates from Behr and Becker (2018), and the variation in upper plate thickness mimics variable plate ages.

Global compilation studies show that sediment thickness goes from 0–12 km (Laske et al., 2012; Dutkiewicz et al., 2015; Straume et al., 2019). Sediment thickness variation in our model setup is thus on the upper limit. When the sediment viscosity is equal to the weak crust viscosity (η_sed=η₀), it is analogous to a no-sediment case. The thickness and viscosity of sediments will control the subduction interface shear stresses, while the upper plate thickness will control the interface length.

Diagnostic parameters. The diagnostics or system response parameters, calculated from each simulation, are also categorized into three groups: (1) subduction (slab) diagnostics (including convergence rate, radius of curvature, and trench velocity), (2) interface (sediments) diagnostics (including wedge angle and width of the accretionary wedge), and (3) upper plate diagnostics (including trench depth as well as maximum and mean topographic amplitudes in the upper plate). These variables are compatible with parameters derived from statistical analyses of present-day subduction zones (i.e., Clift and Vannucchi, 2004; Lallemand et al., 2005; Wu et al., 2008; De Franco et al., 2008; Heuret et al., 2012), which will be used to constrain and validate the results of numerical models. Diagnostic parameters are calculated at every time step. Their evolution (steady-state or transient) will constitute the basis of our parameter analysis in Sect. 3.3.

A schematic of how slab parameters, such as plate convergence (u₀), trench retreat (u_T), and radius of curvature (R_C), are calculated is shown in Fig. 1b. The convergence velocity is calculated as the horizontal motion between the subducting plate and upper plate (i.e., motion between Marker 1 and Marker 2 in Fig. 1a). The trench retreat is calculated as $u_{\mathrm{T}}=\mathrm{d}{x}_{\mathrm{T}}/\mathrm{d}t$, where x_T is the trench position and t is time, starting from an initial trench position x_T=0. In all simulations, the trench is retreating, which is common to ocean–ocean subduction in nature (Schellart et al., 2007).

The radius of curvature is one of the parameters that requires more careful inspection. The radius of curvature is generally calculated from a circular fit to the available earthquake data, which for natural subduction systems can be noisy, incomplete, or subjective (Buffett and Heuret, 2011; Lallemand et al., 2005; Wu et al., 2008). For example, some analyses fit earthquake data for the first 100 km or less because the plate interface is clearly marked at high resolution. However, there are limits on the length of the arc of a circle that can reasonably fit a unique circle (i.e., see discussion in Val and Willenbring, 2020). For this reason, slab dip angle is often used as a metric for slab orientation. However, radius of curvature is more appropriate to describe slab deformation with depth, while slab dip represents only the tangent to curvature close to the surface.

We calculated the radius of curvature after Petersen et al. (2017) (Fig. 1c) in the following way: extract the upper surface of the core of the slab (black) and fit a circle to an arc defined by the inflection point where the plate starts bending (red point below A) and the point on the surface corresponding to 150 km depth (red point next to C). The slab core is the most appropriate feature for the fitting algorithm to calculate the radius of curvature of the slab, as its strength controls the bending of the slab. The upper layers (weak crust and sediments) may deform strongly during subduction and introduce noise into the circle-fitting algorithm. This algorithm remains robust throughout the evolution of a simulation (see movies in the data repository). Therefore, our calculations of the radius of curvature are approximately 30 km less than the total radius of curvature, which includes the crust and sediments.

Sediments reaching the trench may either subduct into the mantle or accumulate into an accretionary wedge. We quantify this behavior using two diagnostics for the accretionary wedge: the angle (α_wedge) and the width (W_wedge) (Fig. 1c). These parameters are not equivalent to the ones calculated in the taper-wedge theory (Dahlen, 1984; Dahlen et al., 1984) because surface topography variations in the wedge are not resolved for the current numerical resolution.

The algorithm to calculate both the wedge angle and width is the following: (1) isolate the sediment markers (grey material in Fig. 1c) and (2) determine wedge points (A, B, C): point A is the inflection point of the slab at the surface, point B is the end point to the right of the surface of sediments, and point C is at the base of the upper plate. We then connect ABC into a triangle and calculate ${\mathit{\alpha }}_{\text{wedge}}=\textcolor{red}{}\mathrm{ACB}$ and $W_{\text{wedge}}=\stackrel{\mathrm{‾}}{\mathrm{AB}}$. The Supplement and movies in the data repository show examples of the algorithm. It is important to note that the wedge angle is nonzero in all simulations, as there is a finite thickness of the sediments (i.e., points A and B will not overlap).

The effect of sediments on topography is also investigated. Lamb and Davis (2003) suggested that sediment-starved subduction may have been responsible for high topography in the Andes. They argue that tectonic erosion favors more coupling with the upper plate, while accretionary wedges favor decoupling and thus lower topographic amplitude. We investigate these hypotheses by extracting three diagnostics related to topography: trench depth (h_trench), maximum topographic amplitude in the upper plate (h_max), and mean topographic amplitude in the upper plate (h_mean). The choice of last two is motivated by the study of Pusok and Kaus (2015), which showed that the two parameters can describe a number of topographic expressions for convergent margins.

Figure 1 Model setup. (a) The model consists of an oceanic plate (SP) subducting beneath another oceanic upper plate (UP). Both plates are 2500 km long and have an upper mantle lithosphere component with a 20 km strong core. The subducting plate contains a weak crust and sediments (combined 15 km; weak crust – dark green, sediments – light green). Parameter values are listed in Table 1 and are relative to the reference density (ρ₀=3300 kg m⁻³) and viscosity (${\mathit{\eta }}_{\mathrm{0}}=\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{20}}$ Pa s) of the mantle asthenosphere. The red markers (Marker 1 and Marker 2) are located in the strong cores of the subducting and upper plate to track the convergence of the plates. The domain has free-slip boundary conditions on all sides and a 60 km “sticky air” layer on top of the plates. Subduction is initiated by prescribing an initial slab depth of 200 km and an initial radius of curvature of 150 km. The dynamics are entirely controlled by the negative buoyancy (slab pull) of the subducting slab. (b) Zoom area on the subduction interface illustrating the convergence rate (u₀), trench rate (u_T), and the radius of curvature (R_C). (c) Phase markers used to calculate the radius of curvature (core markers – black) and wedge properties such as angle and width (sediment markers – grey).

Table 1 Simulations performed in this study. Margin type: TC – tectonic coupling, low-AW – low-angle accretionary margin, high-AW – high-angle accretionary margin, u-AW – unstable accretionary margin. Reference viscosity is ${\mathit{\eta }}_{\mathrm{0}}=\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{20}}$ Pa s (upper mantle viscosity). Diagnostics (mean values): u₀ – convergence rate, R_C – radius of curvature, u_T – trench retreat, α_wedge – wedge angle, W_wedge – wedge width, h_trench – trench depth, h_max – maximum topographic amplitude in the upper plate, h_mean – mean topographic amplitude in the upper plate. The evolution of each diagnostic parameter is shown in the Supplement.

Table 2 Input and diagnostic parameters. Input parameters are prescribed at the beginning of each simulation and remain the same throughout the evolution, while diagnostics are system-response parameters and are calculated during model evolution. Subscripts represent the following: SP – subducting plate, sed – sediments, UP – upper plate. Input parameters: ρ – density, η – viscosity, h – thickness. Diagnostics: u₀ – convergence velocity between subducting plate and upper plate, u_T – trench motion, R_C – radius of curvature, α_wedge – angle of accretionary wedge, W_wedge – width of accretionary wedge, h_trench – trench depth, h_max – maximum topography in the upper plate, h_mean – mean topography in the upper plate.

3 Results

The 2-D numerical experiments below aim to understand what causes convergent margins to either accrete material delivered by the subducting plate or, alternatively, to subduct the trench sediment pile and even erode the basement of the overriding plate. In the first part of the results, we describe end-member models of margin styles (accretionary and tectonic coupling) and the corresponding reference model results. In the second part, we analyze results from a large suite of numerical models and investigate parameter correlations using the diagnostics presented above.

3.1 Margin styles and reference models

The outcome of each simulation is classified into three regimes: tectonic coupling (TC), low-angle accretionary wedge (low-AW), and high-angle accretionary wedge (high-AW) (Fig. 2, Table 1). The end-member division was done both qualitatively (i.e., formation of the accretionary wedge as seen in Fig. 2, left column) and quantitatively using the evolution of diagnostic parameters, which is shown in Fig. 3. We explain the individual classification below. Details of subduction dynamics for each model evolution are included in the Supplement, as they have been extensively investigated and described in previous studies (i.e., adjustment of the model to initial conditions and development of slab curvature, formation of accretionary wedge, interaction of slab with the lower mantle).

We focus here on characterizing overall model outcomes from incipient subduction to slab consumption. The time taken for subduction to consume the slab varies for each model because the characteristic velocity (i.e., convergence rate) is controlled by interface processes (see dependence of average time step size on convergence rate in the Supplement, Fig. S7). For this reason, we use a characteristic time ($t_{\text{char}}=t/t_{\text{final}}$) to highlight the entire model evolution in several figures (i.e., Fig. 3). Initial model time corresponds to t_char=0, while final time t_final, corresponding to ∼2500 km slab consumption, becomes t_char=1. Every model simulation has an initial adjustment period of t_char∼0.1 in which the subduction system acquires a natural curvature and a final stage t_char∼0.9 in which the slab is consumed (grey intervals in Fig. 3). These initial and final condition stages are excluded in our calculation of diagnostic parameters.

For the rest of the study, we adopt a color code for each regime: purple for tectonic coupling, yellow for low-angle accretionary wedge, and orange for high-angle accretionary wedge.

Tectonic coupling margin (TC). Figure 2a–c shows a typical simulation outcome of a tectonic coupling margin (model SubdSed03, with thin cover h_sed=5 km and high viscosity of sediments η_sed=η₀). Evolution of diagnostic parameters is shown in Fig. 3, and additional model evolution snapshots are shown in the Supplement in Fig. S3. We classified simulations as TC when the evolution of the radius of curvature, convergence rate, wedge angle and width, and trench rate remain constant (steady-state) throughout the simulation (purple lines, Fig. 3 and Supplement Figs. S13–S15). The means during the entire evolution of the model (dotted purple lines) closely follow the diagnostics.

The low convergence velocity and trench retreat rate maintained throughout the simulation (Fig. 3b, e) suggest a high degree of coupling between the subducting slab and upper plate. Figure 2b shows that the motion between the subducting and upper plate in the TC margin is accommodated in the middle of the sediment layer, a region of high strain rates. Instead, very low strain rates just below the slab suggest a rigid-body rotation of the slab in order to maintain a constant radius of curvature. The radius of curvature remains small as seen in Fig. 2c, with a steeply dipping slab.

Entrainment of sediments within upper plate material at the interface (Fig. 2a) is indicative of some erosion of the upper plate. All simulations with high-viscosity sediments show this behavior (Supplement, Figs. S13–S15, cases with η_sed=η₀), which could be regarded as having a stronger interface (i.e., a more mafic cover and/or lack of weak unconsolidated sediments).

Topographic signals in TC models (trench depth, maximum and minimum topography in the upper plate) show more variability in Fig. 3, which is discussed later.

Accretionary wedge margin (AW). When the viscous coupling is reduced, sediments are scrapped off the subducting slab to form an accretionary wedge (Fig. 2d–i). We identify two types of accretionary margins: low- and high-angle accretionary wedges, primarily controlled by the thickness of sediments in cases of low-viscosity sediments. Evolution of diagnostic parameters is shown in Fig. 3 (yellow and orange lines), and additional model evolution snapshots are shown in the Supplement in Figs. S4–S5.

The distinction between the two cases comes from the behavior of the slab; low-angle AW margins have increasing radii of curvature and wedge properties but a roughly constant convergence rate (Fig. 3b), while high-angle AW margins result in flat subduction with large radii of curvature and irregular behavior of the convergence rate. In high-angle AW simulations, in a first stage of evolution, accumulation of sediments lubricates the interface and promotes fast convergence rates, but once the wedge reaches a critical size and slab curvature is too large (i.e., subduction needs to accommodate horizontal slab motion), the subduction rate is inhibited. When this stage is reached, plate convergence may happen at slower rates than in TC simulations.

Accretionary margins models are favored by lower sediment viscosities (${\mathit{\eta }}_{\text{sed}}=\mathrm{0.01}\times {\mathit{\eta }}_{\mathrm{0}}$). By increasing the thickness of the sediments, more sediment is available to create a thicker wedge (high angle and width). The larger the wedge angle, the larger the radius of curvature, suggesting that wedge geometry has a control on slab bending (Fig. 3a–d).

The interface dynamics are also different with increasing availability of weak sediments. The velocity field at the subduction interface suggests internal counterclockwise flow inside accretionary wedges, detached from both corner flows in the mantle (Fig. 2e, h). Plate motion is also accommodated at the base of the sediment wedge, indicated by high strain rates. Despite the unrealistically large geometry of the wedge, Fig. 2h shows that the accretionary wedge can impede motion within the upper plate. However, we are not concerned here with further details of internal wedge dynamics compared to numerous previous studies because we lack the numerical resolution required (i.e., Ruh, 2017; Menant et al., 2020).

A highly unstable accretionary wedge margin is shown in the Supplement (e.g., results of SubdSed04_100 in Figs. S6, S13–S15). In these cases, the accretionary wedge reaches a critical angle, and instead of moving laterally, material is being expelled down the subduction channel. The wedge will continue to deform and grow again until it reaches a new critical angle. This unstable mode occurs in simulations with a thick upper plate, which acts as a deformable backstop, in combination with accumulating weak sediments. This margin style is a consequence of the density model chosen here, in which sediments have the same density as the rest of the lithosphere. Clearly, this is an overestimated effect.

Figure 2 Margin style end-member results. (a–c) Tectonic coupling margin (TC, reference model SubdSed03), (d–f) low-angle accretionary margin (low-AW, reference model SubdSed01), and (g–i) high-angle accretionary margin (high-AW, reference model SubdSed04). The left column shows model snapshots at 27 Ma. The middle column shows the second invariant of the strain rate, ${\dot{\mathit{ϵ}}}_{II}=({\dot{\mathit{ϵ}}}_{ij}{\dot{\mathit{ϵ}}}_{ij}{)}^{\mathrm{1}/\mathrm{2}}$, and velocity field as arrows in an enlarged area of the subduction interface. The right column shows core and sediment markers with calculated radii of curvature, wedge angles, and widths. The TC margin shows low-angle wedges and a small radius of curvature. When an accretionary wedge forms in cases of low-viscosity sediments, the wedge angle and width increase over time, together with the radius of curvature. In high-angle AW margins (g–i), episodes of flat subduction with strong plate bending may occur. The time evolution of diagnostic parameters for the reference cases is shown in Fig. 3.

Figure 3 Evolution of diagnostic parameters for end-member reference models: tectonic coupling – SubdSed03 (purple), low-angle accretionary margin – SubdSed01 (yellow), high-angle accretionary margin – SubdSed04 (orange). (a) Radius of curvature, (b) convergence rate, (c) wedge angle, (d) wedge width, (e) trench rate, with positive values indicating trench advance and negative values indicating trench retreat, (f) trench depth, (g) maximum topographic amplitude in the upper plate (h_max), and (h) mean topographic amplitude in the upper plate (h_mean). Continuous lines represent model data, while dotted lines represent the mean over the entire simulation time. Characteristic time t_char is defined in the main text. Grey portions represent initial and final conditions, corresponding to the system forming its natural slab curvature and to the last stage of slab consumption. The two stages are excluded from calculating the means.

3.2 Major impacts on subduction dynamics

Although material within the weak layer at the plate interface (i.e., between the subducting and upper plates) is a volumetrically insignificant component of the larger plate–mantle coupled system, we observe that this small feature can exert a profound influence on the emergent regional-scale subduction dynamics. Figures S10 and S11 show that models with the strongest sediment layers (η_sed=η₀) have sub-vertical slab morphologies, smaller values of R_C, and slower trench retreat rates than similar models with weaker sediment layers (${\mathit{\eta }}_{\text{sed}}=\mathrm{0.01}\times {\mathit{\eta }}_{\mathrm{0}}$). Stronger sediment layers also stabilize the subduction system as seen in Fig. S12, in which much larger variations in R_C occur for models in column A than column C, where values remain approximately constant. The steady-state values of R_C can vary by more than a factor of 2 due to the viscosity of the sediments, with model SubdSed06_50 having 400 km, while models SubdSed04_50 and SubdSed05_50 both evolve such that R_C exceeds 800 km. For these models, a similar increase of more than a factor of 2 can also be observed in convergence rate (Fig. S14) and trench motion (Fig. S15).

It is not just the low strength of the sediment layer that influences the system, but also the thickness of the upper plate. The dynamics of the plate interface depend on the total length of contact area between the two plates as well as the thickness and strength of the material between them (Beall et al., 2021). Comparisons of models that are otherwise identical except for having different upper plate thicknesses (Fig. S11, column A) exhibit more shallow-angled and variable slab morphologies for those models with thinner upper plates.

3.3 Parameter analysis

In this section, we investigate correlations between the means of the diagnostic parameters for all simulations (Figs. 4 and 5). The diagnostic parameters can then be compared to similar parameters observed in the global subduction system.

Figure 3 shows that the evolution means characterize the margin style of each simulation. For each diagnostic parameter, we calculate the mean value (Fig. 3, dotted lines) and the variability during evolution (minimum and maximum values). The means in TC models remain close to the evolution curves (i.e., steady state with the mean close to the min/max values). In AW models, the means differ significantly from the evolution curves. Therefore, the min/max values during time evolution reflect this larger variability (grey bars in Figs. 4 and 5). Diagnostic means are given in Table 1, and evolution curves are shown in the Supplement (Figs. S13–S18).

Figure 4 shows correlations between subduction and sediment diagnostic parameters, while Fig. 5 shows correlations between sediment and topography diagnostic parameters. In both figures, each colored point represents the mean value in a simulation, and the grey bars represent variability intervals (min/max values). The color scale represents the margin type, identified as in the previous section. TC simulations (purple) have smaller variability bars, emphasizing steady state. On the other hand, AW models have larger variability bars (largest for high-angle AW), suggesting transient evolution for that diagnostic.

We find strong correlations between sediment parameters and subduction parameters in 2-D numerical models (Fig. 4). Convergence rate correlates well with radius of curvature, wedge angle, and width (panels a–c). TC models have low convergence rates and radii of curvature, but also small accretionary wedge properties (angle and width). As seen in previous section, with increasing sediment availability (thickness) and/or decreasing sediment viscosity, wedges form more readily. Figure 4 shows that AW models register higher convergence rates, higher radii of curvature, and larger wedge properties. This indicates that sediments lubricate the subduction interface, promoting faster convergence rates. The rate of trench motion is also influenced by the dynamics of subducting sediments, with sediment-rich AW trenches retreating faster (panel d).

Figure 5 suggests weaker correlations between sediment and topography parameters. In panels (b)–(j), TC models tend to produce higher maximum topography and deeper trench depths (i.e., extreme amplitudes), while AW models produce lower extreme signals but higher mean values. The weaker correlation of sediments diagnostics with topographic signals is likely due to the nature of the upper plate (i.e., oceanic or continental, free or attached plate to side walls). In our set of simulations with ocean–ocean subduction and an unattached upper plate, stresses accumulated at the plate interface are accommodated in trench motion instead of upper plate deformation. Topographic signals for simulations with a fixed upper plate are shown in Fig. S9, which are higher and more distinct for end-member models, as subduction is now accommodated more in deformation of the upper plate and less in trench motion. Topography also builds faster in continental lithosphere compared to oceanic lithosphere due to less resistance against gravity (Pusok and Kaus, 2015). Investigation of the effect of sediments on topography signals with a continental upper plate is reserved for a future study.

To summarize, TC models suggest a stronger coupling at the plate interface, yielding lower radii of curvature, slower convergence rates, and higher topographic signatures. TC margins also retreat slower compared to AW margins. These correlations and simulation snapshots suggest that dynamics of the wedge control the bending of the slab and the radius of the curvature.

Figure 4 Subduction versus sediment diagnostic parameters. Each point represents the evolution mean of the given parameter in a simulation, and the grey bars represent the variability intervals (see Fig. 3 for how the mean is calculated). Initial and final stages of the evolution are removed from calculations. TC models (purple) have lower variability in diagnostics than AW models (yellow and orange). TC margins have low convergence velocities, small radii of curvature, and wedge properties. AW models have higher radii of curvature, faster convergence velocities, and larger wedge properties.

Figure 5 Sediment versus topography diagnostic parameters. Each point represents the evolution mean of the given parameter in a simulation, and the grey bars represent the variability intervals (see Figs. 3 and 4). TC models yield smaller trench depths and mean topography, as well as higher maximum topography in the upper plate. The weaker correlations are due to the nature of the upper plate, which is considered oceanic and unattached to the right boundary; any topographic signal is smaller than it would be for continental upper plate lithosphere. Topographic signals for simulations with a fixed upper plate are shown in Fig. S9, which are higher and more distinct for end-member models.

4 Discussion

Our results show three modes of subduction interface: tectonic coupling margin, low-angle accretionary wedge margin, and high-angle accretionary wedge margin. We obtain a diverse response in subduction geometry to just a few varied input parameters: sediment viscosity and thickness, as well as upper plate thickness. The focus of our analysis is solely on the effect of low-viscosity sediments on large-scale subduction dynamics. Other parameters have also been shown to be important (i.e., density of sediments, thermal structure of plates) that will be discussed below.

The viscosity of sediments represents the critical parameter, and thickness is a secondary parameter in our simulations. High-viscosity sediments lead to TC margin, while low-viscosity sediments lead to AW margins (Table 1). Sediment thickness controls the availability of sediments to be accumulated in accretionary wedges. A detailed analysis of Figs. S10–S18 shows that the thickness of the upper plate plays an important role in determining the subduction interface length and the depth at which sediments can be locked into an accretionary wedge. For the same thickness of sediments, when the upper plate is thinner, the accretionary wedge is volumetrically smaller but wider and with a higher angle, while for a thick upper plate, sediments are distributed across the entire interface into a thinner layer indicative of tectonic coupling (i.e., compare results of SubdSed02_50 for a 50 km thick upper plate and SubdSed02_150 for a 150 km thick upper plate in Fig. S11).

4.1 Parameter correlation and observations

A common approach to constrain the effect of different parameters that control subduction dynamics has been done considering the statistical analysis of present-day subduction zones (Supplement, Fig. S2) (Clift and Vannucchi, 2004; Lallemand et al., 2005; Wu et al., 2008; De Franco et al., 2008; Heuret et al., 2012). Our diagnostic parameter analysis in Sect. 3.3 builds a bridge between these studies and numerical models, as diagnostics from numerical models are often not comparable with those from statistical analyses. A number of fundamental differences remain between our results and these studies. Present-day subduction parameter correlations represent only current snapshots and do not always account for the evolution stage of a given subduction system. Our analysis considers the entire evolution of a subduction system. Moreover, natural subduction zones have variable sediment influxes, changing boundary conditions, and multiphase physics. On the other hand, the numerical model setup is ideal and simplified, which is further discussed below.

Despite these differences, we succeed in obtaining subduction margins that accumulate sediments in accretionary wedges, as well as ones that do not and instead conserve the subduction interface geometry. The low accretion margin obtained in numerical results for tectonic coupling (TC) is assumed to be analogous to the tectonic erosion observed in nature. The limitations of this assumption are discussed in the next section. Our results then agree with findings in Clift and Vannucchi (2004). Accretionary margins form in simulations with thicker and weaker sediment cover entering the trench, while tectonic coupling margins form in simulations with less sediment cover (high-viscosity sediments are representative of a stronger crust and mantle component at the slab surface). Replotted data from Wu et al. (2008) and using the margin classification from Clift and Vannucchi (2004) (Supplement, Fig. S2) show that accretionary margins also tend to have larger radii of curvature, which is consistent with our numerical results. A recent study obtained similar correlations between margin style and radii of curvature (Brizzi et al., 2021).

Our findings using slab curvature are consistent with observations on slab dip. Diversity of subduction zones is generally investigated through the perspective of slab orientation (i.e., Beall et al., 2021; Riel et al., 2018). Many analogue and numerical studies used slab dip as the preferred diagnostic for slab orientation instead of slab curvature (i.e., Heuret et al., 2007, Sect. 2.2). They suggest that slab dip depends on properties of the subducting plate, mantle, overriding plate, and the coupling between the subducting and overriding plate (e.g., Bellahsen et al., 2005; Heuret et al., 2007; Billen and Hirth, 2007; Schellart et al., 2007; Funiciello et al., 2008; Babeyko and Sobolev, 2008; Duarte et al., 2013; Riel et al., 2018). We find that coupling at the plate interface due to sediment influxes can strongly influence bending of the slab, exemplified by the radius of curvature. TC models have small R_C (large dip), while AW models have large R_C (small dip). The nature of the upper plate (oceanic or continental), however, will further influence these correlations by changing the load on the subduction interface, which should be investigated further.

Topographic signals have been incorporated less in both statistical and numerical studies of subduction zones, despite the fact that topography is a direct and easily acquirable observable (i.e., Pusok and Kaus, 2015; Riel et al., 2018). In higher-resolution models, Menant et al. (2020) found that deep accretion processes influence forearc topography over time but did not compare the effect in both accretionary and erosion margins. In our ocean–ocean subduction, we see a lesser control of sediments on topographic diagnostics, with TC margins yielding deeper trenches. However, the effect increases when the upper plate is attached to the domain wall (Supplement, Fig. S9). This is because the subduction interface stresses are transferred in the upper plate (i.e., topography build-up) rather than accommodated in trench retreat.

4.2 Convergence rate and margin type

Results from numerical models suggest strong correlations between margin type and convergence rate (Fig. 4). Convergence rates in AW models are faster than in TC models, as sediments help lubricate the interface and reduce coupling between the subducting and upper plate. An exception is the high-angle AW model results (Fig. 3b), wherein the convergence rate is faster in the growing stage of the accretionary wedge, but once it reaches a critical value, the convergence rate becomes slower than in the case of tectonic coupling. Behr and Becker (2018) obtain similar predictions, in which weaker (lower viscosity) sediments promote faster convergence rates. However, observations of the global subduction system suggest an inverse correlation between convergence rate and margin type (Supplement, Fig. S2), with tectonic erosion margins subducting faster than accretionary wedge margins. Brizzi et al. (2021) obtained lower convergence rates for accretionary margins but analyzed numerical results only at a time instance in model evolution. The inverse correlation between convergence rate and margin type in the global subduction system may also suggest that we do not capture all complex processes happening at the subduction interface (von Huene et al., 2004).

The dominant mechanisms governing accretionary and tectonic erosion margins are different. The style of deformation within accretionary wedges is thin-skinned; that is, associated thrusts and folds are uncoupled from the underlying oceanic crust by a basal decollément with a large amount of displacement (Weiss et al., 2018; Angiboust et al., 2021). The dynamics of accretionary wedges are more complex if one considers multiple types of sediments, fluid pressure, and deformation mechanisms (Ruh, 2017; Menant et al., 2020). In accretionary margins, ocean sediments are available to lubricate the interface, while little sediment enters tectonic erosion margins. However, in tectonic erosion margins, the subducting plate may erode the basement of the upper plate and create further debris along the plate interface (von Huene et al., 2004). In any given system, both processes may be occurring simultaneously, either in time and space or at the same time in different parts of the subduction zone (Clarke et al., 2018; Comte et al., 2019; Ducea and Chapman, 2018; Straub et al., 2020). At the physics level, these interface processes are governed by hydraulic fracturing, melting, and deformation of multiphase media, which are not accounted for in our numerical models.

Moreover, it is unclear whether sediment influxes affect convergence rate, convergence rate affects sediment accumulation at the trench, or both. Accretion is generally favored by slow convergence (<7.6 cm yr⁻¹) (Silver et al., 1985; von Huene and Scholl, 1991a; Clift and Vannucchi, 2004; Syracuse et al., 2010), while fast convergence favors larger volumes of sediment being dragged down at the interface, thus lubricating the interface. Erosive margins tend to occur in regions where the convergence rate exceeds 6 cm yr⁻¹ and where the sedimentary cover is thin (Straub et al., 2020). Testing this hypothesis requires different experimental setups: one in which convergence rate evolves dynamically in response to sediment input and a second in which the convergence rate is prescribed and sediments deform accordingly. Here, we tested the first scenario, in which the convergence rate is the dynamic response of the system due to incoming sediments at the trench. We briefly tested the second scenario in Supplement Figs. S8–S9 with kinematic boundary conditions. Imposed convergence with a fixed upper plate reduces the amount of trench retreat, and the variation of diagnostic parameters is smaller among end-member cases, but the correlations remain valid. Menant et al. (2020) investigated the dynamics of accretionary margins as a function of convergence rate to find that more sediment and crustal material are subducted at higher convergence rates.

4.3 Sediment fluxes to trench and to depth

Sediment fluxes to depths below the lithosphere influence the quantity of volatiles recycled into the mantle (Plank and Langmuir, 1998; van Keken et al., 2011; Plank and Manning, 2019). The type of margin may also affect how much sediment gets subducted into the mantle. It has been proposed that tectonic erosion margins can subduct higher percentages of sediment influx; however, large volumes of continental crust are subducted at both erosive and accretionary margins (von Huene and Scholl, 1991a; Clift and Vannucchi, 2004). Clift and Vannucchi (2004) calculated that accretion is a relatively inefficient process for cleaning sediment off the oceanic basement and that 70 % of the sediment column is likely subducted to great depths below the forearc.

Here, we calculate the volume fraction of sediments accreted and subducted in the mantle below the lithosphere in our reference models (Fig. 6). We find that TC margins subduct a higher percentage of influx sediments than AW margins. The total percentage of sediments subducted right before the slab was consumed in TC margins was 60 %, while in AW margins the percentage remains at ∼10 %. However, considering that TC margins have a smaller sediment cover (and influx), accretionary margins may in total subduct a larger volume of sediments. For example, 50 % of a 1 km column of influx sediments is equivalent to 10 % of a 5 km column of influx sediments, which gives a 500 m column of subducted sediments. We suggest that both tectonic erosion margins and accretionary wedges can subduct a large volume of sediments, but at different rates relative to influx material.

This is, however, a simplistic view of sediment subduction to depth. As sediments approach the trench, they can be accreted in the trench, subducted into the mantle, accreted structurally to the bottom of the upper plate, or returned to the upper plate via magmatism, diapirism, or some other form of partially molten upward transport (Gerya and Yuen, 2003; Currie et al., 2007; Tian et al., 2019). Modeling all these processes requires future development of multi-physics models that is beyond the scope of this study.

Figure 6 Evolution of the total volume fraction of sediments accumulated in the accretionary wedge and subducted to depth in the mantle for the three reference cases: (a) tectonic coupling (SubdSed03), (b) low-angle accretionary wedge (SubdSed01), (c) high-angle accretionary wedge (SubdSed04).

We also considered a constant sediment flux to the trench. Sediment subduction is neither a steady-state nor a globally averaged process (Plank and Manning, 2019) and can have major implications for subduction dynamics. Modern oceanic sediments cover 70 % of the planet’s surface, but sediment distribution and lithology occur in drastically different proportions globally (Dutkiewicz et al., 2015, 2018). Moreover, the oceanic lithosphere is covered by various sediment types depending on the depth, proximity to continental margins, and interactions with the oceanic currents and biosphere. For example, an abundant carbonate cover is subducted at the Central American margin, while little sedimentary carbonate is subducted along the Tonga, central Aleutian, and Kuriles–Kamchatka trenches (Plank and Langmuir, 1998; Plank and Manning, 2019). The global sedimentary cover also varies in both space and geological time, with the greatest volume being geologically recent and decreasing exponentially with increasing age (Peters and Husson, 2017).

We expect the abundance and lithology of sediments at trenches to influence the mode of occurrence of margin styles in space and geologic time. The subduction interface structure and properties are sensitive to the composition of the material that is being subducted (Behr and Becker, 2018). Our results broadly cover the end-member scenarios. The degree to which variability in these influxes impacts long-term subduction dynamics remains debated (Cloos and Shreve, 1988; Duarte et al., 2015; Behr and Becker, 2018) and should be studied with further numerical modeling.

4.4 Model limitations

In order to be able to address the points above, a number of model improvements are needed. The models shown here provide an initial experiment on the viscosity of sediments, wherein we considered constant sediment fluxes at the trench. Future work should explore the effect of other sediment properties (density, lithology) and active surface processes such as erosion and sedimentation for potential delivery of continental sediments to the trench. The effect of sediments is overestimated in some of our model outcomes, especially in models with unrealistically large wedges. Thinner and lighter sediment cover requires higher resolutions at the trench, which could also help connect results with those from high-resolution accretionary wedge models (Menant et al., 2020).

Despite our simplified model, mechanical coupling between plates should not be investigated only in variations in subduction velocity or dynamics, but also in the type of margin: accretionary or erosion. We propose that the definition of subduction interface in numerical models could be relaxed. A recent methodological analysis (Sandiford and Moresi, 2019) investigated the emergence of significant spatial and temporal thickness variations within the interface zone, with the sole focus of maintaining a constant thickness through time. Instead, Beall et al. (2021) and Riel et al. (2018) also investigated variations in interface thickness as a potential factor controlling the dynamics. All these studies, in fact, highlight the tendency for the subduction interface to develop spontaneous thickness variation as the models evolve. The interface widens near the trench, building a prism-like complex, and thins at depths beyond the brittle–ductile transition. This pattern was also noted in the boundary element models of Gerardi et al. (2019), who attributed a down-dip thickness variation to lubrication layer dynamics. However, natural observations suggest that variable sediment fluxes enter the subduction trench, questioning this approach of a constant-thickness interface layer.

5 Conclusions

Systematic 2-D numerical simulations of ocean–ocean subduction are run to investigate how low-viscosity sediments influence subduction dynamics and plate coupling. The aim is to understand what causes convergent margins to either accrete material delivered by the subducting plate or, alternatively, to subduct the trench sediment pile and even erode the basement of the overriding plate. We obtain end-member cases that are governed primarily by sediment viscosity and thickness: accretionary and tectonic coupling margins. Results from a large suite of numerical simulations represent a continuum between these end-member models. We find that the properties of the sediment layer modulate the extent of viscous coupling at the interface between the subducting and upper plate. When the viscous coupling is increased, a tectonic coupling style of margin is favored. On the other hand, when the viscous coupling is reduced, sediments are scrapped off the subducting slab to form an accretionary wedge. The geometry of the wedge controls the bending of the slab and the radius of the curvature. We perform an automated analysis of diagnostic parameters to differentiate between the two end-member modes of margin type and to better understand fundamental differences between them. Strong correlations between sediment, subduction diagnostics, and margin type are observed. Tectonic coupling margins have smaller radii of curvature, wedge parameters, and slower convergence rates, while accretionary margins are dominated by larger sediment wedges that can strongly influence subduction dynamics. With some limitations, tectonic coupling margins obtained here are analogous to tectonic erosion margins observed in nature.

A more careful examination of sediment properties and other complex processes happening at the subduction interface in both accretionary and erosive margins is necessary. The margin type is intimately linked to earthquakes. The amount of sediment filling the trench was proposed to facilitate seismic rupture (Heuret et al., 2012; van Rijsingen et al., 2018; Brizzi et al., 2020). Subduction zones with large amounts of trench sediment positively correlate with the occurrence of great interplate earthquakes.