Jonathan R Weiss
email: j.r.weiss@leeds.ac.uk

Garrett Ito

Benjamin A Brooks

Jean-Arthur Olive

Gregory F Moore

James H Foster

Formation of the frontal thrust zone of accretionary wedges

The combination of recent megathrust earthquakes that have ruptured the frontal portions of accretionary prisms and historical evidence for large wedge-front, fold-and-thrust belt events has sparked interest in the development of wedge-front fault systems. Here we explore the formation of these faults using two-dimensional finite-difference models of spontaneous thrust fault formation and wedge accretion. A reference model with uniform material strength parameters predicts wedge-front thrusts to nucleate near the surface of the sedimentary layer and propagate down towards the basal décollement. The formation of a new frontal thrust is preceded by incipient faulting over a zone several kilometers wide forward of the wedge toe, resembling the protothrust zone (PTZ) imaged in seismic reflection data from the front of submarine accretionary prisms. Similar behavior, with a new frontal thrust forming from top to bottom, occurs if the fraction (l) of pore-fluid pressure relative to the total (dynamic) pressure is uniform and < 0.8 throughout the sediment section, or if l increases with depth and is < 0.7 at the base of the décollement. For greater values of l, models predict more frontal thrusts to form from the bottom to top. These results are reasonably well explained by a simple elastic stress analysis. Models that simulate cohesive brittle strength that is negligible in the shallow sediments and increases with depth also predict new thrust faults to form from top to bottom. Models in which viscous flow at shallow depths represents creep of sediments predict bottom-up thrust nucleation when viscous stresses due to plate convergence remain lower than the brittle yield stress in the top ~15-50% of the sediment section. Top-down nucleation is otherwise favored when viscous stresses exceed the brittle yield strength throughout most of the sediments. Our study illuminates key controls on the mechanical and temporal links between the PTZ, the forward propagation of slip along the underlying décollement, and the formation of a new frontal thrust.

accretionary prism; fold-and-thrust belt; protothrust zone; thrust fault formation, numerical modeling; wedge-front deformation 1. Introduction

In the aftermath of recent large earthquakes associated with low-angle megathrust faults that underlie accretionary wedges, our understanding of how and where wedges rupture has begun to change. Previously, the frictional properties of the shallow up-dip portions of megathrusts were thought to be more likely to facilitate aseismic rather than coseismic slip [START_REF] Scholz | Earthquakes and friction laws[END_REF]. The 2010 Mw 7.8 Mentawai, 2011 Mw 9.0 Tohoku-Oki, and the 1999 Mw 7.6 Chi-Chi earthquakes [START_REF] Fujiwara | The 2011 Tohoku-Oki Earthquake: Displacement Reaching the Trench Axis[END_REF][START_REF] Hill | The 2010 Mw 7.8 Mentawai earthquake: Very shallow source of a rare tsunami earthquake determined from tsunami field survey and near-field GPS data[END_REF][START_REF] Yue | Structural geology of a classic thrust belt earthquake: the 1999 Chi-Chi earthquake Taiwan (Mw=7.6)[END_REF], however, are just a few examples of recent events characterized by coseismic slip that advanced to the wedge toe and, in the case of the Japan and Sumatra earthquakes, resulted in devastating tsunamis. These events demonstrate the need to better understand the mechanics of the wedge toe including the formation of the décollement and steeply dipping thrust faults that facilitate surface uplift (e.g. [START_REF] Hubbard | Coseismic slip on shallow décollement megathrusts: implications for seismic and tsunami hazard[END_REF].

Our knowledge of wedge-front structures comes primarily from seismic reflection and well-log data, and in the case of fold-and-thrust belts, outcropping faults, folds, and geomorphic indicators. Such observations are often interpreted with geometric and kinematic models (e.g. fault propagation folding, trishear) that assume thrusts initiate near the tip of the basal décollement, propagate upward towards the surface, and are linked to a décollement at depth (e.g. [START_REF] Allmendinger | Estimation of fault propagation distance from fold shape: Implications for earthquake hazard assessment[END_REF][START_REF] Hughes | Insights into the mechanics of fault-propagation folding styles[END_REF][START_REF] Suppe | Geometry and kinematics of fault-propagation folding[END_REF]. In contrast, an observation that has received relatively limited attention is the faulting associated with incipient deformation forward of the main frontal thrusts. High-resolution seismic reflection data outboard of accretionary prisms including Cascadia, Nankai, and New Zealand, for example, image a "protothrust zone" (PTZ), characterized by numerous cross-cutting faults midway through the sediment section that often offset the seafloor but do not clearly intersect the décollement (e.g. [START_REF] Barnes | The role of protothrusts in frontal accretion and accommodation of plate convergence, Hikurangi subduction margin[END_REF]Karig, 1986;MacKay, 1995) (Fig. 1).

This paper investigates the mechanics and formation of the wedge-front thrust fault system using a combination of a simple stress analysis and finite-difference models of spontaneous fault nucleation and growth. Our results show that in the simplest case of a sediment layer with homogeneous material strength parameters, new wedge front faults initiate at shallow depths and then propagate down to the décollement. We then vary material properties to further map out the range of conditions needed to produce top-down versus bottom-up fault propagation. One group of models explores different uniform fractions (l) of pore-fluid pressure relative to the total pressure, another group includes a more realistic increase in l with depth, another considers an increasing cohesive brittle strength with depth. Finally, we explore models in which viscous creep is possible at shallow depths in the crust, and the strain rate-dependent strength increases with depth. We conclude with a discussion of the implications for natural systems.

Numerical Approach

We model the formation and evolution of orogenic wedges using SiStER (Simple Stokes solver with Exotic Rheologies, [START_REF] Olive | The role of elasticity in simulating long-term tectonic extension[END_REF], a two-dimensional MATLAB ®-based code that uses finite differences on a fully staggered mesh with the particle-in-cell method [START_REF] Gerya | Introduction to Numerical Geodynamic Modelling[END_REF] to solve for conservation of mass and momentum in an incompressible visco-elastic-plastic continuum. At each time step, Picard iterations are used to solve the non-linear terms in the governing equations based on updated estimates of the strain rate-dependent rheology, computed on the finite difference grid. This process is repeated until the L2-norm of the strain rate field residual is £10 -3 , or a maximum of 75 iterations are performed. Between time steps (50% of the Courant condition), Lagrangian tracer particles passively track material properties and stresses as they are advected in the Eulerian velocity field using fourth-order Runge-Kutta. Material properties are passed between nodes/cells and particles using bilinear interpolation [START_REF] Gerya | Introduction to Numerical Geodynamic Modelling[END_REF]. Cells have dimensions of 125 m × 125 m for all model runs unless otherwise noted in Table S1, which lists all important model parameters.

The Maxwell model is used for the visco-elastic deformation, but the stress is limited by the plastic strength, which leads to the formation of localized shear bands that represent faults. Plastic strength is defined by the Drucker-Prager form of the Mohr-Coulomb failure criterion in which brittle strength (as the second invariant of the 2-D deviatoric stress tensor) increases with total pressure P, friction angle f, and cohesion C, 𝜎 "" = 𝑃𝑠𝑖𝑛𝜙 + 𝐶𝑐𝑜𝑠𝜙.

(1) Brittle failure is computed on the mesh nodes rather than particles. Brittle strength reduction occurs as accumulated plastic strain ep increases from 0 to 0.05 (e.g. [START_REF] Lohrmann | The impact of analogue material properties on the geometry, kinematics, and dynamics of convergent sand wedges[END_REF] by decreasing C linearly from C0 = 20 MPa at ep = 0.0 to Cmin = 0.01 MPa at ecrit = 0.05. In addition, fault healing is simulated by reducing ep over a time scale (t ≈ 63 kyr) that is small compared to the duration individual faults are active. Shear bands typically localize to widths of 2-4 grid cells. Thus, with 125 m grids, the critical plastic strain required for full fault weakening roughly corresponds to 10-20 m of slip.

A traction free upper surface on the sediments is simulated by covering the crustal layer with a "sticky air" layer of low viscosity (ha=10 17 Pa s) and near-zero density. Topography of the sediment-air interface is tracked through time using a chain of markers that are advected passively by the velocity field. A crude representation of erosion and mass wasting using linear topographic diffusion (e.g. [START_REF] Simpson | Mechanics of non-critical fold-thrust belts based on finite element models[END_REF] with a low diffusivity coefficient of k = 5×10 -9 m 2 s -1 preserves much of the faulted topography while adequately smoothing unrealistically sharp or overhanging features.

A contracting wedge configuration (Figs. 2a,S1) is created with the right side of the model box being free-slip and having an imposed leftward velocity of Vx = -10 mm yr 1 . The impenetrable bottom boundary is also imposed to move leftward at Vx = -10 mm yr -1 . The left side "backstop" is rigid and impenetrable, and the top (of the sticky air) is free-slip and open. This layer has an initially uniform thickness of H = 4.5 km, and overlies a basal décollement having an initial thickness of Hd = 0.5 km. The décollement has the same material properties as the overlying crust except C0 = 0.01 MPa and the friction angle is lower (fb = 15°). The evolution of the numerical wedge begins with the development of two prominent shear bands (i.e., thrust faults) near the left-hand-side of the model domain that form a v-shaped conjugate pair, intersecting at the basal décollement (Fig. S1). Shortly thereafter, the forwardverging (dipping to the left) thrust dominates, behind which topography develops and the wedge forms. Growing topography increases horizontal compression in the incoming crustal layer and eventually new fault formation occurs ahead of the frontal thrust leading to wedge widening. The cycle of wedge thickening, with slip on the frontal and internal thrusts, and wedge widening, with sequential formation of new frontal thrusts repeats. In this manner, the wedge grows and preserves an average slope consistent with critical taper theory [START_REF] Davis | Mechanics of fold-and-thrust belts and accretionary wedges[END_REF], as shown in prior numerical (e.g. [START_REF] Simpson | Mechanics of non-critical fold-thrust belts based on finite element models[END_REF][START_REF] Stockmal | Whence the Mountains? Inquiries into the Evolution of Orogenic Systems: A Volume in Honor of Raymond A. Price[END_REF] and laboratory modeling studies (e.g. [START_REF] Gutscher | Episodic imbricate thrusting and underthrusting: Analog experiments and mechanical analysis applied to the Alaskan accretionary wedge[END_REF]Koyi, 1995) (Figs. 2a, S1).

Reference model with uniform material properties and no

Insight to where and how new frontal thrusts form is revealed by the pattern of plastic strain as well as the differential stress (i.e., difference between the least s1 and most compressive s3 stress, where compressive stress is negative) normalized by the local brittle strength,

𝐹 = (0 1 20 3 ) [2(0 1 60 3 )789:6;<=>7:] (2)
where the term in the denominator is the Drucker-Prager failure criterion. Note also that cohesion C in Eq. ( 2) can locally decrease where plastic strain accumulates and can be ≤ C0. The actively straining shear bands beneath the wedge are by definition at failure (F = 1). Prior to the formation of a new frontal thrust, a ~10-km-wide region in which F > ~0.95 forms near the sediment surface outboard of the active wedge front (Fig. 2b). As time proceeds, this high-F region enlarges and a series of incipient thrusts form at conjugate angles within this shallow region (Fig. 2c). While this is occurring, slip on the décollement has propagated forward well beyond the main frontal thrust (t = 2.167 Myr; Fig. 2c), indicating that the material above it is moving toward the backstop slower than the convergence rate (see also [START_REF] Morgan | Effects of cohesion on the structural and mechanical evolution of fold and thrust belts and contractional wedges: Discrete element simulations[END_REF]. Hence this material is shortening horizontally, and this shortening is accommodated in the incipient protothrust zone (PTZ). With continued contraction, a new frontal thrust emerges from the PTZ at or near the surface and progressively lengthens down-dip toward the décollement (Fig. 2d). When the new thrust fault meets the décollement, slip on the décollement ceases beyond the intersection location (t = 2.214 Myr; Fig. 2e). Varying the basal friction does not change the basic behavior of the new frontal thrust propagating from top to bottom.

Ideal elastic stress models

To understand the cause of the top-down propagation we analyze the state of stress in the incoming sediment layer. Far forward of the deformation front, including near the right side of the model domain, the state of stress is nearly lithostatic (Figs. 3a,3b,S2). However, within a distance L (~25 km) ahead of the frontal thrust, the stresses are perturbed by the wedge. Inclined contours of horizontal normal stress sxx indicate that sxx becomes increasingly compressive both downward and towards the wedge; roughly horizontal contours of sxz indicate a downward increase (in absolute value) from zero at the surface to the frictional strength of the décollement at depth.

The prominent characteristics of this stress pattern can be understood with a simple model of force balance on a purely elastic body of rectangular cross-section having height H and width L. The top of the rectangle is stress free; the right boundary has a horizontal normal traction equal to the lithostatic pressure; the bottom has a uniform, right lateral shear traction tb = rgHtanfb (g is gravitational acceleration); the left side has a rightward horizontal normal traction in excess of lithostatic that, when integrated along the side, balances the net leftward force due to the basal shear (Figs. 3d,S2a). The solution is found using the Airy stress function (also used by [START_REF] Goff | Stresses beneath a ramping thrust sheet[END_REF], see Supplement for derivation), 𝜎 @@ = -

B C D 𝑥 -𝜌𝑔𝑧, 𝜎 JJ = -𝜌𝑔𝑧, 𝜎 @J = B C D 𝑧. (3)
Here, z is depth (>0 downward), and the horizontal coordinate x = 0 is on the right boundary at a distance L forward of the wedge toe and increases to the left (i.e., x = L at the wedge toe).

The distance L forward of the wedge over which stresses are significantly perturbed from the background lithostatic state can be estimated by approximating sxx beneath the surface intersection of the frontal thrust as being at critical failure (Dahlen, 1990). Given that the active frontal thrust is fully weakened, the differential stress profile is approximated by the Mohr-Coulomb failure limit without cohesion

0 1 20 3 ; = -K 0 1 60 3 ; L 𝑠𝑖𝑛𝜙. (4)
If the least compressive stress is approximately vertical, s1 » szz, and if the most compressive stress is approximately horizontal, s3 » sxx, then 𝜎 @@ ≈ -K N6789:

N2789: L 𝜌𝑔𝑧.

(5)

Integrating Eq. ( 5) over depth gives the net force of the wedge pushing the block to the right. Balancing this with the leftward shear traction along the base, integrated across the distance L yields

𝐿 = K N6789: N2789: L K PQD R ;B C L = K N6789: N2789: L K D ;ST9(: C ) L. (6)
To verify this prediction, we estimate L in the SiStER models as the distance beyond the frontal thrust over which the model shear stress at a given depth is approximately constant. The SiStER values of L for different décollement friction angles fb are well-explained by Eq. ( 6) (Figs. 4,S3). Thus, the elastic stress solutions in Eq. ( 2) are applied over the domain of x = 0 to L (again x = L is at the hypothetical wedge frontal thrust) (Figs. 3c,d). The solutions qualitatively match the general stress pattern predicted by SiStER, specifically the nearly horizontal invariance of |sxz| and its increase with depth toward the décollement where the shear stress is imposed on the layer, as well as the approximately linear increase in |sxx| with depth and toward the deformation front. To illuminate the cause of the top-down, frontal thrust propagation, the stresses given by Eq. ( 3) are used to compute the differential stress (s1-s3) as a function of x and z for comparison with the corresponding brittle failure envelopes. Depth profiles at various locations of x show that (s1-s3) exceeds the brittle strength near the surface first, and the thickness of the failure zone increases right-to-left toward the hypothetical wedge front (Fig. 3f). This behavior resembles that of the SiStER models (Fig. 3e).

The simple elastic and reference numerical models predict top-down faulting for two main reasons: (1) the depth-dependent failure law leads to decreasing brittle strength with proximity to the surface and (2) the differential stress decreases more gradually toward the surface. The theoretical tendency for top-down fault propagation was one of numerous aspects discussed in prior analytical studies (e.g. [START_REF] Eisenstadt | Alternative model of thrust-fault propagation[END_REF][START_REF] Goff | Stresses beneath a ramping thrust sheet[END_REF], and is often apparent in numerical and analogue modeling studies (e.g. [START_REF] Bernard | Kinematics of fault-related folding derived from a sandbox experiment[END_REF][START_REF] Del Castello | Underthrusting-accretion cycle: Work budget as revealed by the boundary element method[END_REF][START_REF] Dotare | Initiation of a thrust fault revealed by analog experiments[END_REF][START_REF] Morgan | Effects of cohesion on the structural and mechanical evolution of fold and thrust belts and contractional wedges: Discrete element simulations[END_REF]. This result contrasts with the notion that new frontal thrusts often nucleate from the décollement tip and propagate upwards (e.g. [START_REF] Allmendinger | Estimation of fault propagation distance from fold shape: Implications for earthquake hazard assessment[END_REF][START_REF] Chester | Fault-propagation folds above thrusts with constant dip[END_REF][START_REF] Tavani | A review of deformation pattern templates in foreland basin systems and fold-and-thrust belts: Implications for the state of stress in the frontal regions of thrust wedges[END_REF] and thus motivates further exploration.

Effects of pore-fluid pressure 4.1 Constant pore-fluid pressure

What qualities are needed to produce top-down, versus bottom-up fault propagation? We first explore the importance of brittle-strength reduction due to pore-fluid pressure [START_REF] Dahlen | Noncohesive critical Coulomb wedges: an exact solution[END_REF][START_REF] Ruh | Effect of fluid pressure distribution on the structural evolution of accretionary wedges[END_REF]. We do this in the numerical models by assuming pore pressure is a constant fraction l of the total pressure (i.e. l = pf/P) in the Drucker-Prager brittle failure law 𝜎 "" = 𝑃(1 -𝜆)𝑠𝑖𝑛𝜙 + 𝐶𝑐𝑜𝑠𝜙.

For all models we apply the same l value to both the sediment layer and décollement, but C=Cmin (0.01 MPa) in the décollement. We note that while we simulate the effects of pore fluid on the frictional strength the material, we do not incorporate the density contrast between the fluid and solid when solving the momentum equation.

Numerical model predictions

Numerical simulations predict that the behavior of new frontal thrusts can change due to l. With l < ~0.8 the new frontal thrusts grow much like those from the reference case (l = 0): from top to bottom. Recall that for l = 0, the failure ratio F forward of the frontal thrust is highest in a zone near the surface, and it is in this zone that the new frontal thrust initiates before propagating downward (Fig. 2). With increasing l the depth extent of the zone of high F increases. For example, with l = 0.5, a zone of high F forms approximately midway through the sediment layer but as the zone broadens, failure still initiates in the upper half of the layer and the frontal thrusts propagate top-to-bottom (Fig. S6, Anim. S4-S8). When l = ~0.8, however, a different behavior emerges. Stresses forward of the frontal thrust first become largest in a localized zone near the décollement. With continued convergence, the zone of high F expands up and forwards (Fig. 5).

The associated PTZ is characterized by incipient shear bands forming near both the surface and the base of the sediment section. One of these shear bands becomes the new, dominant frontal thrust fault, sometimes forming nearly simultaneously across the whole layer but often propagating from the décollement upwards (Fig. 5, Anim. S8).

Elastic stress models

We now consider the effects of pore-fluid pressure on our simple model of elastic stresses. The solutions are still given by Eq. ( 3), but, pore-fluid pressure also reduces the basal shear stress according to 𝜏 X = 𝜌𝑔ℎ(1 -𝜆)𝑡𝑎𝑛𝜙 X .

(8) The second adjustment is to L, the domain of x over which this basal traction acts and the elastic solutions are computed. To compute L, the approximate failure stress near the active thrust (Eq. 4) becomes

\ (0 ] 1 20 ] 3 ) ; \ ≈ -K (0 1 60 3 ) ; + 𝑝 _ L 𝑠𝑖𝑛𝜙 (9)
where 𝜎 `N = 𝜎 N + 𝑝 _ and 𝜎 `a = 𝜎 a + 𝑝 _ are the least and most compressive effective principal stresses based on Terzaghi's principle [START_REF] Skempton | Terzaghi's discovery of effective stress[END_REF]. With Eq. ( 8), Eq. ( 5) becomes

𝜎 @@ ≈ -K N6789:(N2;b) N2789: L 𝜌𝑔𝑧, (10)
and Eq. ( 6) becomes

𝐿 = K N6789:(N2;b) N2789: L K PQD R ;B C L = K N6789:(N2;b) N2789: L K D ;ST9(: C )(N2b) L. (11)
Thus, although increasing pore pressure reduces the horizontal compressive force near the frontal thrust (Eq. 9), the effects are largest at the base, which controls tb and so L increases with increasing l (Fig. 4).

The elastic stress solutions Eq. ( 3) are again used to predict the differential stresses and failure envelopes, but this time with tb and L given by Eqs. ( 8) and ( 11). The resulting depth profiles qualitatively reproduce the behavior of the SiStER models that increasing l diminishes the difference between the failure and differential stress curves (Fig. 6). For low l, when the failure strength still increases with depth much more quickly than the differential stress, both models predict failure to initiate near the surface. For higher l's, the two curves approach each other, and in the SiStER model they become approximately parallel at l = 0.8. At this value of l, the SiStER model initiates brittle failure over a greater thickness of the sediment layer, including near the décollement (Figs. 5,6). Compared to the SiStER model, the elastic model under-predicts the downward increase in differential stress and predicts the curves to be parallel only when l = 1 (instead of 0.8). We conclude that top-down thrust propagation should dominate for low-tomoderately high pore-fluid pressure ratios (l < 0.8), whereas faults can initiate over the full range of depths, including near the décollement only at very high ratios (l ≥ 0.8) 4.4 Depth-dependent pore-fluid pressure Recent work suggests that more typically, pore-fluid pressure ratios increase with overburden pressure (i.e. depth) (e.g. Moore and Saffer, 2001;Zhang, 2013). Often the gradient in porepressure is nearly hydrostatic near the surface and increases towards the lithostatic gradient at depth, leading to an effective brittle strength that increases non-linearly and tends to flatten out with depth in the crust [START_REF] Suppe | Fluid overpressures and strength of the sedimentary upper crust[END_REF]. Recently, Ruh (2017) used numerical models of wedges to characterize how changes in l influence the overall wedge geometry. They found that steady increases in l can cause a décollement to form spontaneously (i.e., the décollement needs not be imposed as in our models) and the suggestion is that new frontal thrust faults may also initiate at depth although this topic is not explicitly described.

We take an approach similar to [START_REF] Ruh | Effect of fluid pressure distribution on the structural evolution of accretionary wedges[END_REF] and simulate depth-varying l. Here, the porepressure gradient is hydrostatic (l = rw/rs = 0.37; see Table S1) at shallow depths, and then shifts to a lithostatic gradient, continuing to the base of the sediment section (Fig. 7) [START_REF] Suppe | Fluid overpressures and strength of the sedimentary upper crust[END_REF]. We ran 6 models in which the pore-fluid pressure ratio reaches different values (l3.5) at a depth of 3.5 km, which corresponds to the initial base of the décollement. The friction angle of the décollement was decreased with increasing l3.5 to preserve the ratio of the effective friction of the décollement relative to the mean effective friction in of the initial sediment layer (Table S1). The maximum cohesion was set to C0 = 5 MPa, which is much lower than in previous models. This is in response to models with depth-varying l in which individual thrust faults developed unrealistically large offsets for greater values of C0 for the higher values of l3.5.

Figure 7 shows the results of a model with l3.5 = 0.8. In this case, pore-fluid pressure deviates from hydrostatic values at a "fluid retention depth" (e.g. Suppe, 2014) of zFRD ~1.1 km beneath the surface. The PTZ is characterized by the gradual formation of two zones of high failure ratio F; one zone expands upward from the décollement and the other expands downward from the surface. The two zones coalesce near z= zFRD (t = 0.768; Fig. 7, Anim. S10). This behavior contrasts with the previously discussed model of constant l = 0.8 in which the zone of failure initially nucleates only above the décollement. In the l3.5 = 0.8 model, incipient frontal thrust faults consistently form at the base of the sediment section (t = 0.831 Myr) before rapidly propagating up-dip (t = 0.834 Myr) and intersecting the surface (t = 0.858 Myr). Models with l3.5 ≤ 0.7 also exhibit the same high-F zone coalescence but the new frontal thrusts form near the surface and propagate down-dip towards the décollement (Anims. S8-9). We also note that higher values of l3.5 produce increasing dominance of forward-verging thrusts, a behavior that was first documented by [START_REF] Ruh | Effect of fluid pressure distribution on the structural evolution of accretionary wedges[END_REF].

Additional depth-dependent material properties 5.1 Depth-dependent cohesion

Closely associated with depth-dependent pore-fluid pressure is a decrease in pore-fluid volume with depth (e.g., Yue and Suppe, 2014). One likely effect of depth-varying porosity is on the bulk cohesive strength associated with brittle deformation. We explore this effect by including a depthdependence on the initial cohesion (prior to any strain-weakening) according to

𝐶 c = 𝐶 dT@ e1 -exp f-K PQ(J2J g ) PQJ h L i jk. ( 12
)
Here zs is the depth of the sediment surface and z0 is the characteristic depth over which cohesion remains low (Fig. 8). The form of Eq. ( 12) is motivated by the aim to have a smoothly varying function that has cohesion being near zero at the surface, remaining low over the characteristic depth, and increasing below this depth, approaching Cmax = 20 MPa at the base of the sediment layer. We ran four models with depth scales of z0 = 0.5, 1, 1.5, and 2 km (Figs. 8, S9).

Figure 8 shows the resulting numerical predictions for z0 = 2 km: cohesion is ≤ 2 MPa (10% of Cmax) in the upper ~1/3 of the sediment layer and is ³ 18 MPa (90% of Cmax) in the bottom third of the layer. The model predicts that the shallow, low-cohesion layer remains near failure over large distances behind, as well as in front of the frontal thrust (Fig. 8, 0.692 Myr). In front of the wedge, F is always high near the surface, which contrasts the reference case where zones of high F form for a more limited time just prior to the nucleation of a new frontal thrust. With continued contraction, the high-F zone grows down toward the décollement, and numerous, primarily forward-verging thrusts nucleate within this near-surface zone (Fig. 8; 0.941 Myr). This behavior contrasts with the reference model where only a few cross-cutting shear bands form near the surface. Soon thereafter a single, dominant fore-or back-thrust propagates downward to the décollement (Fig. 8; 0.966-0.985 Myr). Thus, although the shallow layer with low cohesion deforms diffusely most of the time, new frontal thrusts always nucleate from within this layer and propagate downward to the décollement. The models with thinner low-cohesion layers behave the same way, except that the zone of diffuse deformation and closely spaced shear bands becomes thinner with decreasing z0 (Fig. S9).

Depth-dependent viscosity: shallow viscous creep in the upper crust

To this point our models have focused on simulating elastic and brittle deformation. However, time-dependent creep or ductile deformation is known to be an important mechanism in the upper crust for accommodating compaction, folding, and shearing, often related to faulting and fracturing. Rate-dependent processes include pressure-solution creep (e.g. [START_REF] Croizet | Chapter 3 -Compaction and Porosity Reduction in Carbonates: A Review of Observations, Theory, and Experiments[END_REF][START_REF] Gratier | The Role of Pressure Solution Creep in the Ductility of the Earth's Upper Crust[END_REF] and brittle creep (e.g. [START_REF] Atkinson | Subcritical crack growth in geological materials[END_REF][START_REF] Chester | Subcritical creep compaction of quartz sand at diagenetic conditions: Effects of water and grain size[END_REF]. Here we explore the effects of strain-rate-dependent deformation by approximating the non-brittle deformation as a viscous rather than a primarily elasto-plastic process. Pressure solution creep results in viscous behavior that depends on effective stress, temperature, and porosity among other factors. Here we assume a simple linear stress dependence, neglect temperature effects, and assume porosity (i.e. effective viscosity) varies only with depth. Our objective is to define the regimes of frontal thrust evolution for different depth variations in strain-rate dependent stress.

As our models do not simulate volumetric strain, we approximate a strain-rate dependent stress with an effective shear viscosity that varies as the following function of depth z below the surface depth zs,

𝜂 = 𝜂 > exp p-𝑓 N K rQD rQ(J2J g )6_ R rQD L s t ( 13
)
where h0 is the reference viscosity, and f1, f2, and g are constant input parameters (Fig. 9, Table S2). The term in parentheses in the exponential is the "porosity function", which represents the relative change in porosity with depth (Fig. S10), displaying a rapid decrease near the surface and a more gradual decrease at greater depths as is commonly observed within accretionary prisms (e.g. Moore et al., 2001;Moore and Saffer, 2001). The exponential function arises in creep laws for shear deformation in partially molten mantle rocks during diffusion creep [START_REF] Kohlstedt | 2.18 -Constitutive Equations, Rheological Behavior, and Viscosity of Rocks A2 -Schubert, Gerald, Treatise on Geophysics[END_REF] and in pressure solution creep laws describing compaction of sediments [START_REF] Niemeijer | Compaction creep of quartz sand at 400-600°C: experimental evidence for dissolution-controlled pressure solution[END_REF][START_REF] Noort | Kinetic effects of microscale plasticity at grain boundaries during pressure solution[END_REF]. In both situations, the exponential dependence on porosity accounts for the decreased intergranular contact area (and proportional increased stress on grain contacts) with increasing porosity. By varying f1, f2, and g we change the variation in viscosity with depth. However, the more fundamental and generally applicable parameter that we vary is the ratio Fa of the approximate, ambient differential stress imposed by the converging underlying plate, relative to the yield stress,

𝐹 T = ;uv w ⁄ ;[789:PQ(J2J g )6< h =>7:] . ( 14
)
The numerator is the approximate ambient (strain-rate dependent) differential stress with U = 10 km/Myr and L = 25 km as defined by Eq. ( 6).

We ran calculations for 14 different depth variations in Fa (Fig. 9, Table S2). Results show that models in which Fa ≤ 1 in the upper 85-100% of the sediment layer never produce brittle failure, and the wedge deforms entirely viscously. Cases in which Fa > 1 at depths greater than 0 to ~10% of the starting sediment layer thickness predict brittle failure with frontal thrusts propagating from top to bottom much like in the reference model. Within the intermediate range of Fa tested, the frontal thrust initiates just above the décollement and propagates from the bottom up.

The case shown in Figure 10 (and Anim. S13) is one such case of bottom-up propagation. Here the failure ratio F in the numerical model forward of a mature frontal thrust is seen to be largest at depth where the effective viscosities are highest. Immediately in front of the wedge, the zone of high F thickens to shallower depths, but F remains highest just above the décollement (t = 1.005-1.173 Myr). Eventually, cross-cutting shear bands form at the base (t = 1.173 Myr), and a dominant forward-verging thrust fault emerges and propagates up towards the surface (t = 1.192-1.214 Myr). The strain-rate dependence in the models is essential for bottom-up propagation. For these models with intermediate values of Fa, the stress due to ambient contraction (2hU/L) is too low to initiate faulting in the upper half of the layer. Once a prominent shear band initiates near the base, strain rates and hence differential stresses become magnified near the tip of the shear band and brittle failure extends upward (Fig. 9c). This feedback between magnified strain rate and fault tip propagation continues upward either to the surface or to the level where the viscosities are too low for further brittle failure. In the model show in Figure 10, the shallowest sediment deforms by folding.

Discussion

Model shortcomings

The most significant shortcoming of our models is their simulation of incompressible deformation, whereas matrix compaction (and the concurrent dewatering) can accommodate a substantial portion of the horizontal shortening in accretionary prisms (e.g. [START_REF] Barnes | The role of protothrusts in frontal accretion and accommodation of plate convergence, Hikurangi subduction margin[END_REF]Moore et al., 2001;[START_REF] Morgan | Décollement processes at the Nankai accretionary margin, southeast Japan: Propagation, deformation, and dewatering[END_REF]. Our models therefore are likely to over-predict the amount of shearing on incipient or main thrusts for a given amount of total shortening. This shortcoming is likely to be most important at shallow depths were the initial porosities are lowest, and hence the models over-predict the tendency for protothrusts to nucleate near the surface. This shortcoming is also likely to be more significant for higher-porosity, submarine accretionary prisms than for lower-porosity mountain belts that incorporate crystalline basement rock. Our models also do not simulate fluid flow and the dynamic responses of the partitioning between volumetric versus shear strain, or changes in pore pressure. Thus, the rates of fault formation and propagation, as well as the time-integrated distribution of slip along individual faults are unlikely to be realistic.

Top-down propagation

Models that predict the frontal thrust to initiate near the surface include those with moderate-tolow pore-fluid pressure ratios. Pore pressure has the effect of reducing the rate at which brittle strength increases with depth compared to the rate differential stress increases therefore promoting greater depths of incipient shear-band formation. For the Japan Trench, a critical-taper-based, whole-wedge pore pressure ratio estimate of l = 0.5 [START_REF] Davis | Mechanics of fold-and-thrust belts and accretionary wedges[END_REF] is compatible with Nankai PTZ formation near the top of the sediment section (Figs. 1,6).

For models in which cohesion increases with depth, a shallow, low-cohesion layer displays diffuse deformation, largely accommodated by small, short-lived shear bands. The model example shown in Fig. 8 bears resemblance to both the Nankai and Hikurangi seismic reflection images, which reveal a PTZ composed of closely spaced thrusts that appear to intersect the surface but not the décollement (Fig. 1). The zone of deformation composed of fine shear-bands on the hanging wall side of the fully developed frontal thrust fault (e.g. thrust #2 in Fig. 8) is also reminiscent of small-scale, PTZ-like structures observed in the Hikurangi and Nankai examples.

The minimal size and spacing of the shear bands in the above model, however, is limited by the size of our model grids. In nature, brittle deformation, perhaps in localized bands, may in fact occur on the grain scale and thus be truly diffuse on the scale of stratigraphy or seismic reflection resolution. It is therefore appropriate to question how the model discretization impacts the results. Morgan (2015) used a discrete element numerical method to study the effects of cohesion on faulting in wedges and showed that individual thrust faults tend to be more stable (longer-lived) and more widely spaced when cohesion (uniform across the model) is largest, and thrust faults become less stable and more closely spaced with decreasing cohesion. With this insight, we can expect that as cohesion increases with depth, incipient thrusts should begin to emerge at increasing length scales and stability. Thus, although the SiStER models may not accurately represent the size and spacing of shear bands in the low cohesion layer, the nucleation of shear-bands that are sizeable compared to the scale of stratigraphy or seismic reflection resolution, the downward propagation of these bands toward the décollement, and their transformation to major frontal thrusts should be relatively robust predictions. That is, at least to the extent that prior compaction results in a state where shortening by plastic shearing is dominant.

At orogenic wedges, the depth at which thrusts initiate and the direction in which they lengthen may impact inferences about the associated kinematics and structure. Evidence for shallow nucleation and downward propagation is provided, as previously discussed, by geologic outcrops in the Appalachian Valley and Ridge province (McConnell et al., 1997), where the exposure spans a sufficient depth range to reveal a decrease in fault displacement with depth (Fig. 1), and in the Caledonides of Southern Norway [START_REF] Morley | Fold-generated imbricates -examples from the Caledonides of southern Norway[END_REF]. Also, Nivière et al. ( 2006) document anomalously high rates of local uplift in the outer Western Alps, which they attribute to shallow ramp activation (i.e. new frontal thrust formation) in the incoming sediment section well above the active décollement. In contrast, geologic reconstructions often assume that many of the exposed shallow folds and short-strike anticlines at the front of active accretionary wedges are associated with faults that extend all the way from the basal décollement to near the surface (e.g. [START_REF] Hubbard | Structure and Seismic Hazard of the Ventura Avenue Anticline and Ventura Fault, California: Prospect for Large, Multisegment Ruptures in the Western Transverse Ranges[END_REF][START_REF] Lavé | Active folding of fluvial terraces across the Siwalik Hills (Himalaya of central Nepal)[END_REF]. This assumption is used to estimate fault slip and crustal shortening rates. If instead, the folds are related to incipient faults that have not yet linked to the basal décollement, then inferred total slip and shortening would be less. We suggest that anticlinal folds with small along-strike dimensions are more consistent with shallow faults if the vertical extent of the fault is comparable to the horizontal extent.

Bottom-up propagation

Models that consistently predict the frontal thrust to propagate from bottom to top include cases in which the pore-fluid pressure ratio increases with depth to values ³ ~0.7 at the base of the sediment section. The implication for natural systems is that bottom-up propagation should occur in sediments with a low effective basal coefficient of friction and thus lower surface taper angle.

The other type of model that leads to more definitive nucleation of faults near the décollement and upward propagation are the viscous deformation models in which the converging plate creates ratios of ambient differential stress to plastic yield stress (Fa) that exceed 1 starting at depths of ~15-50% of the sediment thickness. Again, our models simulate viscous creep accommodating shear strain not volumetric strain. In natural porous sediments, substantial shear deformation by pressure-solution creep becomes apparent at temperatures of ~150°C or more (e.g. [START_REF] Fagereng | Subduction megathrust creep governed by pressure solution and frictional-viscous flow[END_REF][START_REF] Gratier | The Role of Pressure Solution Creep in the Ductility of the Earth's Upper Crust[END_REF]. These temperatures, however, are only reached near the décollement of submarine wedges (e.g. [START_REF] Fagereng | Subduction megathrust creep governed by pressure solution and frictional-viscous flow[END_REF]Moore and Saffer, 2001) suggesting that the model profiles of Fa are not likely to be achieved by viscous shear deformation in these systems. Higher temperatures are more likely at shallower levels in subaerial fold-andthrust belts, particularly in backarc settings where geothermal gradients are elevated (e.g. [START_REF] Currie | The thermal structure of subduction zone back arcs[END_REF] and in continent-continent collision zones with high radioactive heat production in the upper crust and rapid exhumation rates (e.g. [START_REF] Cattin | Modeling mountain building and the seismic cycle in the Himalaya of Nepal[END_REF]. However, if pressure-solution creep that accommodates compaction leads to similar ambient stress ratio Fa profiles in subaerial or submarine systems, then our model predictions of brittle failure initiation may be more generally applicable to both submarine and subaerial systems.

The second feature of our viscous flow models is the upward fault propagation due to the concentration of stress near the upper tip of the lengthening faults. If this is not accommodated by volumetric strain in natural systems, then we speculate that distributed brittle creep can lead to the required rate dependence on this process (e.g. [START_REF] Brantut | Time-dependent cracking and brittle creep in crustal rocks: A review[END_REF]. Indeed, the interaction of brittle and ductile deformation is known to be important in fault zones (e.g. [START_REF] Gratier | The Role of Pressure Solution Creep in the Ductility of the Earth's Upper Crust[END_REF], and may be integral to the time-dependent behavior of deep megathrusts of active subduction zones (e.g. [START_REF] Fagereng | Subduction megathrust creep governed by pressure solution and frictional-viscous flow[END_REF].

The model shown in Fig. 10 also predicts the shallowest sediment to always remain too weak for brittle failure to occur. This behavior is reminiscent of the "ductile" deformation that creates large amplitude folds above buried thrust faults and evidence for wedge-front earthquakes that do not rupture the surface but rather are attributed to blind thrusting below more distributed deformation above the fault tip (e.g. [START_REF] Ainscoe | Blind Thrusting, Surface Folding, and the Development of Geological Structure in the Mw 6.3 2015 Pishan (China) Earthquake[END_REF]Stein and Yeats, 1989). 6.4 The relationship between the protothrust zone and the frontal thrust A recent study using multi-channel seismic reflection data examined the detailed structure of the PTZ and frontal thrusts across the Hikurangi accretionary prism, offshore New Zealand [START_REF] Barnes | The role of protothrusts in frontal accretion and accommodation of plate convergence, Hikurangi subduction margin[END_REF]. The key results relevant to this study are that the protothrusts often clearly intersect the seafloor, do not extend to the décollement, and show the largest offsets in the middle to upper half of the sediment layer that decrease both up-and down-dip. Across the frontal thrust, however, displacements reach a maximum at the base of the sediment section and decrease up-section, which can be interpreted as indicating the frontal thrust initiated near the décollement and propagated up. If this is true, then these observations are not clearly reproduced our models: the cases that exhibit bottom-up propagation predict strain in the PTZ to extend across the whole sediment layer (Figs. 5,7) or be highest near the décollement (Fig. 10), not at shallow levels as imaged for the PTZ at Hikurangi. In contrast, the models that exhibit a prominent shallow PTZ always predict frontal thrusts to propagate top-down. The depth-varying l cases, with l3.5 ³ 0.7, are perhaps more consistent with Hikurangi observations as regions approaching brittle failure form both at depth and near the surface before a listric thrust propagates up from the décollement (Fig. 9). The closely spaced, shallow shear bands that typify Hikurangi (and some of our top-down models), however, do not form across the shallow, high-strain regions in these models at the current grid resolution.

We then question whether the frontal thrust necessarily grows out of an incipient protothrust. The suggestion that the protothrusts are compactive shear bands below the critical stress for Mohr-Coulomb failure [START_REF] Barnes | The role of protothrusts in frontal accretion and accommodation of plate convergence, Hikurangi subduction margin[END_REF] while the frontal thrust forms by critical Mohr-Coulomb failure provides a reason for a decoupling between the two. Our viscous deformation models of bottom-up propagation could represent such a behavior if strain-rate dependent compaction produces subcritical stresses in the shallow half the sediment layer but critical stresses below. On the other hand, the dynamical interaction between shear deformation, volumetric strain, and fluid flow on or around the newly forming frontal thrust probably leads to the large variation in slip distribution along the Hikurangi frontal thrust, and this slip distribution may not necessarily, directly reflect upon the depth of origin or sense of propagation of the fault. It may then be possible for the frontal thrust to grow from top to bottom out of the protothrust but then preferentially accumulate greater slip at depth as the frontal thrust matures.

Conclusions

We have used finite-difference numerical models of spontaneous thrust fault formation as well as a simple elastic stress analysis to investigate the formation of the frontal thrust zone of accretionary wedges. As is commonly predicted in other numerical and analogue modeling studies, our reference models with uniform material strength parameters predict new frontal thrust faults to initiate near the surface and propagate downward. Similar behavior occurs with low pore-fluid pressures and with cohesion that is negligible in the shallow sediment layer and increases with depth. These models help to explain evidence from geologic outcrops for decreasing fault throw with depth as well as high-resolution seismic reflection images of protothrust zones characterized by numerous cross-cutting faults in the shallow and/or middle portion of the sediment section that do not clearly intersect the décollement (Fig. 1).

In contrast, if pore-fluid pressure is a large fraction of the lithostatic pressure (≥ 0.8 if l is constant and ≥ 0.7 if l increases with depth), models predict faults to initiate over a full range of depths including directly above the décollement. Faults consistently nucleate near the décollement and propagate towards the surface if differential stress depends on strain-rate and increases with depth such that brittle failure occurs only in the deeper half of the incoming sediment layer. Upward propagation of the new frontal thrust occurs due to enhanced strain-rate and thus stress near the upper tip of the incipient thrust.

All of our models indicate a close mechanical and temporal link between the formation of a broad zone of incipient thrusting, the forward propagation of slip along the décollement and the formation of a new frontal thrust. We always predict the new frontal thrust to grow out of an incipient protothrust, but whether this is the case in nature depends on the relative importance of compactive deformation and the dynamical effects on fluid flow, pore-pressure variations, and plastic strength. We suggest that determining the depths at which new frontal thrusts nucleate in accretionary wedges can help provide constraints on the pore-pressures and hydrologic processes of the system, advance our understanding of ductile and brittle deformation mechanisms, and lead to a more accurate interpretation of surface features that signal the widening of an accretionary wedge.

ACKNOWLEDGEMENTS

We thank Samuel Howell, Eric Mittelstaedt, and Ben Z. Klein for their work in developing the numerous features of SiStER. We also thank Eileen L. Evans, Jessica C. Hawthorne, and Jonas B. Ruh for their comments that helped improve the manuscript. This research benefited tremendously from discussions with Julia K. Morgan about different modeling approaches as well as her reviews of multiple drafts of this manuscript. Ito's effort in our numerical code development benefited from discussions during the Melt in the Mantle program which was hosted by the Isaac Newton Institute (INI) for Mathematical Sciences and supported by EPSRC Grant Number EP/K032208/1. Ito also thanks the Simons Fellowship for helping support of his stay at the INI. This study was supported by grants NSF-EAR-1443317 (J.R. Weiss), NSF-EAR-1118481 (J.R. Weiss and J.H. Foster) and NSF-OCE-1155098 (G. Ito). J.R. Weiss also benefited from the UK Natural Environmental Research Council (NERC) support through the Earthquakes without Frontiers (EwF) project (EwF_NE/J02001X/1_1) and the Centre for Observation and Modeling of Earthquakes, Volcanoes and Tectonics (COMET, GA/13/M/031, http://comet.nerc.ac.uk). This is SOEST contribution XXXX.

REFERENCES CITED

Zhang, J., 2013. Effective stress, porosity, velocity and abnormal pore pressure prediction accounting for compaction disequilibrium and unloading. Marine and Petroleum Geology 45, 2-11. [START_REF] Barnes | Tectonic and geological framework for gas hydrates and cold seeps on the Hikurangi subduction margin, New Zealand[END_REF], and (d) Nankai Trough [START_REF] Moore | Structure of the Nankai Trough Accretionary Zone from multichannel seismic reflection data[END_REF]. Karig (1986) define the PTZ as a region forward of the main thrust characterized by disrupted reflection horizons, diffuse thickening, shortening, and porosity reduction that occurs prior to thrusting. 6) for a case where l = 0. The solid vertical black bars are the range in estimated L measured from the corresponding SiStER models (fb=5°,10°,15°, and 20°; also see Fig. S3 of the Supplementary Information). The dashed black curve is the theoretical length scale L predicted by Eq. ( 11) for a case with fb=15° and variable l. The dashed vertical black bars are the range in estimated L measured from the corresponding SiStER models (l=0,0.2,0.5,0.8). The SiStER L values are generally well explained by the theoretical predictions. C0=25 MPa; Cmin=5 MPa). Top panel shows stratigraphic markers for the model domain at ~1 Myr, after forethrust #4 has accumulated significant strain. Subsequent panels show strain rate (left) and failure ratio (right) at different time steps leading to the initiation of forethrust #4. At 0.724 Myr an area at the base of the sediment section outboard of fault #3 begins to approach failure. This zone broadens up and outward and after 0.799 Myr almost the entire protothrust region, from the surface to the décollement is close to failure. At 0.842 Myr, a high-strain shear band forms above the décollement, propagates up-dip, and intersects the surface at 0.849 Myr. High strain rates (10 -13 yr -1 ) are contoured in the bottom 3 failure panels to draw attention to new frontal thrust initiation. An animation of the associated model run is included in the Supplement. 5; see Figure S6 and supplemental animations for l=0.0/0.2/0.5 examples). The Airy profiles are taken from comparable locations (i.e. prior to brittle failure; x position labeled in plots). Both the numerical model and stress analysis predict a decrease in the absolute value of stress and a decrease in the slopes of the failure and differential stress curves with increasing l so that the wedge front region is closer to failure over a greater depth range with increasing l. Therefore, at higher values of l a new thrust can initiate anywhere in the sediment layer including near the décollement (e.g. Fig. 5).

FIGURES

Figure 7

The effect of depth-varying pore-fluid pressure ratio (l) on the formation of a new frontal thrust (l3.5=0.8; f=30°; fb=15°; H=3 km; Hd=0. 5 km; C0=5 MPa; Cmin=0.1 MPa). Top panel shows stratigraphic markers, below left is l for the frontal portion of the wedge at 0.768 Myr. To the right are profiles of fluid pressure, theoretical yield stress (cohesionless), and l for our depth-varying l models. Red profiles in each panel correspond to the model run shown in this figure where l3.5=0.8 at the top of the décollement. Black profiles are models shown in the supplement. Thin black profiles in the "pressure" panel correspond to increasing values of (constant) l from 0.4 to 1.0 in increments of 0.1. Bottom panels show the evolution of strain rate (left) and failure ratio (right) during new frontal thrust formation. High strain rates are contoured in the failure plots to draw attention to the location of new frontal thrust initiation at the base of the sediment section. See the main text and supplement for detailed descriptions of this and related models. (b) Failure ratio Fa defined as the approximate ambient (strain-rate dependent) differential stress normalized by the yield stress for each of the models shown in (a). See section 5.2 of the main text for a detailed description. (c) Yield stress (red) and theoretical differential stress (black and gray) curves for the model shown in Fig. 10. The yield stress curve is identical to that shown in the leftmost panel of Fig. 3f for the elastic stress solution. The solid black curve is the differential stress using lithostatic pressure as s1 and the product of the viscosity profile with the reference contractional strain rate (i.e., the convergence rate of 10 mm/yr divided by L=25 km where L is computed using Eq. 6) as s3. The shallower light gray curves show differential stresses resulting from strain rate increased by factors of (2, 3, 4, …, 10) as occurs with shear band localization (e.g. Fig. 10). Protothrusts initially form in the failure zone (gray) between the décollement and a depth equivalent to ~30% of the total sediment thickness (i.e. ~2 km beneath the surface for this particular model). A rapid increase in strain rate within the shear band subsequently drives fault propagation towards the surface (red arrows; see Anim. S13, S14 and the main text for a detailed description of this behavior). S3 and are essentially independent of G (also see Figure 4). All model parameters are the same as the reference case (Table S1) except the model domain height H, which is 12 km for these runs. c) so that the wedge surface taper is comparable in the 3 models, thus isolating the effect of pore-fluid pressure variations. The top panel shows stratigraphic markers, the 2 nd panel shows strain rate (2 nd invariant), the 3 rd panel shows the failure ratio, the 4 th panel shows yield stress, the 5 th panel shows the total stress (2 nd invariant), and the bottom panel shows pore-fluid pressure l. Animation S14 Close-up of the wedge front region of the depth-dependent viscosity model run shown in Fig. 10 and Animation S13 focusing on the formation of a new frontal thrust. Top panel shows passive stratigraphic markers, the middle panel shows strain rate (2 nd invariant), and the bottom panel shows the failure ratio with contours of elevated strain rate in black. A new thrust fault nucleates near the base of the sediment section, propagates rapidly towards the surface resulting in folding above the fault tip, and punches through the low-viscosity layer at ~1.214 Myr. Values chosen for depth-dependent viscosity models shown in Fig. 9.

c

Values chosen for model shown in Fig. 10. a All other parameters are identical to the reference model b Pore-fluid pressure ratio is constant in the l 3.5 = 0.37 model

  pore-fluid pressure: top-down thrust formation 3.1 Numerical model predictions In the reference model (Figs. 2, 3), the Eulerian mesh is 80 km long and 17 km high. The deformable sediment layer has a friction angle of f = 30°, a starting cohesion of C0 = 20 MPa, a density of rs = 2700 kg m -3 , a viscosity of hs = 10 25 Pa s, and a shear modulus of G = 2 x 10 9 Pa.

Figure 1

 1 Figure 1 Evidence of thrust faults initiating at or near the surface. (a) Line drawing from the Appalachian Valley and Ridge fold-and-thrust belt, showing an outcropping thrust fault where measured displacement (inset) decreases down-section into the footwall syncline (redrawn from McConnell et al., 1997). (b-d) Multi-channel seismic reflection images showing the protothrust zones (PTZ) in front of the (b) Cascadia accretionary prism (MacKay, 1995), (c) Hikurangi subduction zone[START_REF] Barnes | Tectonic and geological framework for gas hydrates and cold seeps on the Hikurangi subduction margin, New Zealand[END_REF], and (d) Nankai Trough[START_REF] Moore | Structure of the Nankai Trough Accretionary Zone from multichannel seismic reflection data[END_REF].Karig (1986) define the PTZ as a region forward of the main thrust characterized by disrupted reflection horizons, diffuse thickening, shortening, and porosity reduction that occurs prior to thrusting.

Figure 2

 2 Figure 2 Predictions of the formation of a new frontal thrust in our reference model (f=30°; fb=15°; k=5×10 -9 ). (a) Whole model domain at t = 2.182 Myr. The brittle sediment layer contains passive yellow and brown markers for visualizing accumulated deformation. The frictionally weak basal décollement is blue. Boundary conditions are illustrated and forethrusts are numbered (1-6) in order of when they first appeared. The bold line above the fault numbers indicates the theoretical surface taper angle (a=5.3°) for a cohesionless wedge (Dahlen, 1990). (b -e) Time steps illustrating the formation of a new frontal thrust where the failure ratio is (s3-s1) divided by the Mohr-Coulomb yield stress. Black contours mark regions of elevated strain rate. Panel (d) is the same time step as (a).

Figure 3

 3 Figure 3 Predicted stress state in front of the wedge for the model shown in Figure 2. (a) and (b) show sxx and sxz, respectively, predicted by SiStER just after a new frontal thrust forms (just after time step shown in Figure 2d). Black contours are in units of MPa and contours of elevated strain rate (bold white) highlight the active frontal thrust and décollement. (c) and (d) show sxx and sxz from the elastic stress model; boundary conditions indicated with thick black arrows in (d). Profiles of the (e) SiStER and (f) Airy differential stresses (solid black) at locations marked by vertical white lines in (a)-(d). Yield stress profiles are solid red, regions of predicted (incipient) brittle failure are gray.

Figure 4

 4 Figure 4 A comparison of L versus basal friction angle fb and pore-fluid pressure l from the SiStER models and simple force balance arguments presented in Section 3.2 and 4.2 of the main text. The solid black curve is the theoretical length scale L predicted by Eq. (6)

Figure 5

 5 Figure 5 Numerical model predictions of the formation of a new frontal thrust under the influence of elevated pore-fluid pressure (l=0.8; f=30°; fb=15°; H=3. 5 km Hd = 0. 5 km;

Figure 6

 6 Figure 6 Differential stress (black) and Mohr-Coulomb failure (red) profiles from just prior to the formation of a new frontal thrust for SiStER (top) and Airy elastic stress models (bottom). Pore-fluid pressure ratios (l) increase from left to right. For the SiStER models, profiles are located 26 km from the backstop at variable model times when stresses are increasing but prior to brittle failure (e.g. 0.724 Myr from the l=0.8 model in Fig. 5; see Figure S6 and supplemental animations for l=0.0/0.2/0.5 examples). The Airy profiles are

Figure 8

 8 Figure 8 Formation of a new frontal thrust for a numerical model with depth-dependent cohesion (z0=2 km). Top left panel shows passive markers at 2.48 Myr. Top right panel shows depth-dependent cohesion curves described in the text. Lower panels show strain rate (left) and failure ratio (right). High strain rates are contoured in the failure plots. Initially a broad region close to failure develops near the surface (0.692 Myr). Cross-cutting shear bands form within this zone (0.941 Myr) and then a single, dominant thrust propagates downwards (0.966 Myr) and intersects the décollement (0.985 Myr). By 1.191 Myr significant strain has accumulated on the new fore-thrust and cross-cutting shear bands have begun to form in the protothrust zone near the surface, outboard of thrust fault #3. Note the similarity between the model protothrust zone and the Nankai and Hikurangi examples in Fig. 1 (i.e. a broad, near-surface zone of cross-cutting thrusts).

Figure 9

 9 Figure 9 (a) Depth dependent viscosity curves used in the numerical models with domains for ductile deformation (i.e. no faulting), bottom-up fault propagation, and top-down fault propagation. Dashed lines are the suite of models used to define the domains. The solid red curve corresponds to the model shown in shown in Fig. 10 (f1=0.2, f2=2.5, g=1, h0=4x10 23 ).

Figure 10

 10 Figure 10 Numerical prediction of the formation of a new frontal thrust for a model with depth-dependent viscosity (f1=2.5, f2=0.2, h0=4×10 23 ). Top right panel shows the failure (red) and differential stress (black) profiles just after the new frontal thrust has nucleated above the décollement (1.173 Myr; X=25 km). Note that much of the sediment layer is near failure except for the near-surface, low-viscosity channel. The other panels are colored and contoured as in Fig. 7. A broad region approaching failure forms near the base of the sediment layer at 1.005 Myr. This region enlarges, and short, cross-cutting shear bands nucleate just above the décollement near x=25 km at 1.142-1.173 Myr. A single, dominant fore-thrust begins to propagate towards the surface at 1.192 Myr, the strain rate within the shear band rapidly increases (1.192-1.214 Myr), and the new thrust fault interacts with the surface at 1.438 Myr. By this time, significant strain has accumulated on the new forethrust resulting in fault-related topography and a new cycle of frontal fault formation begins as the stress outboard of fault #3 increases.

Figure S1

 S1 Figure S1 The top panel illustrates the boundary conditions for our numerical model. Beneath are panels showing the evolution of a representative SiStER orogenic wedge with 250 x 250 m grid cells shown at 2 Myr time steps. The strong, brittle sedimentary layer (f=30°) contains passive yellow and brown markers to aid in visualizing accumulated deformation. The weak basal décollement (fb=10°) is blue. Diffusion-based surface processes (k=5×10 -9 ) simulate erosion/deposition. Also see Animation S1.

Figure

  Figure S2 (a) Cartoon depicting the boundary conditions for the Airy stress analysis. (b) Strain rate just after the formation of a new frontal thrust (same model and time step as Figure 2d from the main text). The white box outlines the area of Figure 3(a) and (b) of the main text, and the vertical white lines are the locations of the three stress profiles also shown in Figure 3(e) of the main text. Color contour plots show (c) sxx and (d) sxz for the same SiStER model run and time step shown in (b) and in Figure 3 of the main text.

Figure S3

 S3 Figure S3 Shear stress (tb) along the décollement in SiStER models for progressively greater (top-to-bottom) basal décollement friction angle. Lmin is the region where tb is non-zero and approximately constant whereas Lmax extends out to where tb is approximately zero. This range in L is plotted as vertical bars in Figure 4 of the main text. The range in estimated L for SiStER models under the influence of pore-fluid pressure were measured in the same way.

Figure S4

 S4 Figure S4Shear stress (tb) along the top of the décollement for SiStER models with progressively increasing (top to bottom) shear modulus G. Overall wedge geometries differ only slightly, and the higher G models propagate outwards more slowly than the reference model (top). Lmin and Lmax are defined the same as in FigureS3and are essentially independent of G (also see Figure4). All model parameters are the same as the reference case (TableS1) except the model domain height H, which is 12 km for these runs.

Figure S5 (

 S5 Figure S5(Top) Theoretical vs SiStER estimates of the characteristic length scale L for models with increasing cohesion Co. Theoretical prediction (thick black curve) show an increase with L with increasing C0 on the order of 5 km for every 10 MPa increase in cohesion. Colored bars are measurements from SiStER models, which also show a detectable increase in L with increasing C0. (Bottom) Shear stress along the top of the décollement for SiStER models with increasing C0 at time steps just prior to the formation of the protothrust zone. The increase in L is most noticeable in the shear stress contours near the thrust front where |sxz| is relatively constant before slowly decreasing across a broad zone towards the right-hand side of the model domain.

Figure S6

 S6 Figure S6SiStER model snapshots illustrating the effect of constant pore-fluid pressure l on the wedge geometry (left) and the nucleation depth and location of the region that initially begins to approach failure ahead of the frontal thrust (right; x=~20 km). The circular region initiates near the surface when l=0 but gradually deepens as l increases until it nucleates just above the decollement when l=0.8. For all models, l=lb, f=30°; fb=15°; H=3. 5 km Hd = 0. 5 km; C0=25 MPa; and Cmin=5 MPa.

Figure S7 (

 S7 Figure S7 (Left) Fluid pressure (middle) yield stress curves (cohesionless) and (right) pore-fluid pressure profiles for models shown in Figs. X and S8.

Figure

  Figure S8 (a-c) Comparison of SiStER models illustrating the effect of variable pore-fluid pressure on overall wedge/thrust geometry and stress state. (a) Model with constant pore-fluid pressure (l=0.37) throughout the wedge. Corresponding pressure, yield stress, and pore-fluid pressure profiles are shown in Figure S7. f=30°; fb=14.6°; H=3.5 km, Hd=0.5 km; C0=5 MPa; and Cmin=0.1 MPa. Basal friction angle fb is varied in models (a-

Figure

  Figure S8 (b) Model with hydrostatic pressure (l=0.37) in the shallow sediment layer and increasing l with depth to l3.5 = 0.5 at the top of the décollement. Also, l increases with depth as the model wedges thicken. Corresponding pressure, yield stress, and porefluid pressure profiles are shown in Figure S7. For models S8a-c, f=30°; fb=17.9°; H=3.5 km, Hd=0.5 km; C0=5 MPa; and Cmin=0.1 MPa. Color scales are identical to Fig. S7.

Figure

  Figure S8 (c) Model with hydrostatic pressure (l=0.37) in the shallow sediment layer and increasing l with depth to l3.5 = 0.8 at the top of the décollement. Also, l increases with depth as the model wedges thicken. f=30°; fb=30°; H=3.5 km, Hd=0.5 km; C0=5 MPa; and Cmin=0.1 MPa. Color scales are identical to Fig. S7.

Figure S9

 S9 Figure S9Depth dependent cohesion model runs for varying depth scales p0. Failure ratios of 1 are contoured to highlight the individual shear bands and their variable penetration depth, which increases with increasing p0. Arrows/brackets on the righthand side of the model domain highlight the increase in the width of the zone of incipient failure with increasing p0.

Figure S10

 S10 Figure S10 Depth-dependent porosity function and viscosity curve used in the numerical model shown in Figure 10 (f1=0.2, f2=2.5, h0=4x10 23 ).

  f d 10 degrees (•) Internal friction angle in the basal layer p o 0.5-2e3ρg MPa Pressure scale that controls the chracteristic low-cohesion depth (Figs. 7, S9, and Eq. 8) a All other parameters are identical to the reference model b Pore-fluid pressure ratio is constant in the l 3.5 = 0.37 model

Table S1c

 S1c Depth-dependent pore-fluid pressure l models a

	Symbol	Value	Unit	Definition
	y max	8	km	Height of the the numerical model
	x max	50	km	Width of the the numerical model
	H a	4.5	km	Initial thicknes of the sticky air layer
	H s	3	km	Initial thicknes of the sedimentary layer
	H d	0.5	km	Initial thicknes of the basal decollement
	C o	5	MPa	Initial (maximum) cohesion of the the sedimentary layer
	C MIN	0.1	MPa	Minimum cohesion
	C d	0.01	MPa	Cohesion of the basal decollement
	l 3.5	variable		Pore-fluid pressure ratio at the base of the sediments
	f d	variable	degrees (•) Internal friction angle in the basal decollement

Table S1d

 S1d

		Depth-dependent pore-fluid pressure l models b
	l 3.5	f d (•)
	0.37	14.56
	0.5	17.94
	0.6	21.13
	0.7	25.07
	0.8	30.28
	0.9	38.19
	Table S1e Depth-dependent cohesion models	a
	Symbol	Value	Unit	Definition
	y max	8	km	Height of the the numerical model
	x max	60	km	Width of the the numerical model
	H a	4.5	km	Initial thicknes of the sticky air layer
	H s	3.25	km	Initial thicknes of the sedimentary layer
	H d	0.25	km	Initial thicknes of the basal decollement

Table S2a

 S2a

		Depth-dependent viscosity models	a,b,c
	Symbol Value Unit	Definition
	y max	10	km	Height of the the numerical model
	x max	60	km	Width of the the numerical model
	H a	5	km	Initial thicknes of the sticky air layer
	H s	4.75 km	Initial thicknes of the sedimentary layer
	H d	0.25 km	Initial thicknes of the basal decollement
	f d	10	degrees (•) Internal friction angle in the basal layer

Table S2b

 S2b Constants controlling the shape of the viscosity profiles and the total variation in viscosity (Figs. 9, S10, and Eq. 9)

						b,c
	Model f 1	f 2	h o	g	Faulting style
	1 c	2.5	0.2	1.5e18	1	bottom up
	2	2.5	0.4	1.0e18	1	no faulting
	3	2.5	0.2	2.0e19	1	top down
	4	2.5	0.4	2.0e22	1	top down
	5	2.5	0.4	2.0e22	0.5 top down
	6	2.5	0.5	1.0e18	1	no faulting
	7	2.5	0.25	1.0e18	1	no faulting
	8	2.15 0.125	1.0e18	1	top down
	9	2.5	0.21	1.0e18	1	bottom up
	10	2.5	0.2	1.0e18	1	bottom up
	11	2.5	0.22	1.0e18	1	bottom up
	12	2.5	0.23	1.0e18	1	bottom up
	13	2.5	0.2	5.0e18	1	bottom up
	14	2	0.125	1.0e18	1	top down

a All other parameters are identical to the reference model (Table

S1

) b

Table S1c

 S1c Depth-dependent pore-fluid pressure l models a

	Symbol	Value	Unit	Definition
	y max	8	km	Height of the the numerical model
	x max	50	km	Width of the the numerical model
	H a	4.5	km	Initial thicknes of the sticky air layer
	H s	3	km	Initial thicknes of the sedimentary layer
	H d	0.5	km	Initial thicknes of the basal decollement
	C o	5	MPa	Initial (maximum) cohesion of the the sedimentary layer
	C MIN	0.1	MPa	Minimum cohesion
	C d	0.01	MPa	Cohesion of the basal decollement
	l 3.5	variable		Pore-fluid pressure ratio at the base of the sediments
	f d	variable	degrees (•) Internal friction angle in the basal decollement

Table S1d

 S1d Depth-dependent pore-fluid pressure l models b Internal friction angle in the basal layer p o 0.5-2e3ρg MPa Pressure scale that controls the chracteristic low-cohesion depth (Figs. 7, S9, and Eq. 8)

	l 3.5	f d (•)
	0.37	14.56
	0.5	17.94
	0.6	21.13
	0.7	25.07
	0.8	30.28
	0.9	38.19
	Table S1e Depth-dependent cohesion models	a
	Symbol	Value	Unit	Definition
	y max	8	km	Height of the the numerical model
	x max	60	km	Width of the the numerical model
	H a	4.5	km	Initial thicknes of the sticky air layer
	H s	3.25	km	Initial thicknes of the sedimentary layer
	H d	0.25	km	Initial thicknes of the basal decollement
	f d	10	degrees (•)

Supplementary information for

Formation of the frontal thrust zone of accretionary wedges

Jonathan R. Weiss 1,2 , Garrett Ito 2 , Benjamin A. Brooks 3 , Jean-Arthur Olive 4 , Gregory F. Moore 2 , and James H. Foster 5 1. COMET, School of Earth and Environment, University of Leeds, Leeds, United Kingdom 2. Dept. of Geology and Geophysics, University of Hawai'i at Manoa, Honolulu, HI, USA. 3. U.S. Geological Survey Earthquake Science Center, Menlo Park, CA, USA. 4. Laboratoire de Géologie, Ecole Normale Supérieure / CNRS UMR 8538, PSL Research University, Paris, France. 5. HIGP, University of Hawai'i at Manoa, Honolulu, HI, USA. * Corresponding author. E-mail address: J.R.Weiss@leeds.ac.uk (J.R. Weiss)

Elastic stress analysis using the Airy stress function

The cause for the stress state observed at the front of the deformation zone in our numerical wedges is examined with a simple model of a homogeneous elastic volume with rectangular cross-section having height H and width L (Figs. 3d andS2a). In other words, the rightward force from the integrated normal traction sxx along the lefthand side of the wedge is such that is balances the leftward force from the basal shear stress tb integrated along the base. Note that for this simple elastic model the horizontal coordinate X increases from right to left towards the wedge front for algebraic convenience.

Using the standard "stress formulation" in elasticity, the equations of continuity and conservation of momentum are expressed as a 2-D biharmonic equation of the Airy stress function y (Timoshenko and Goodier, 1970)

(2) A simple set of solutions to eq. ( 2) is represented by the third order polynomial

which relates to stress according to

.

(3)

The Animation S3 SiStER orogenic wedge model run with l=0.0. Note that the model parameters differ slightly from the reference case shown in Animation S2. See Table S1 for details. The top panel shows stratigraphic markers, the 2 nd panel shows strain rate (2 nd invariant), the 3 rd panel shows the failure ratio, the 4 th panel shows the yield stress, the 5 th panel shows the 2 nd invariant of the total stress, and the bottom panel shows porefluid pressure l. Values chosen for depth-dependent viscosity models shown in Fig. 9.

Animation S4

c

Values chosen for model shown in Fig. 10.