Numerical modeling of tectonic underplating in accretionary wedge systems Open Access

Accretionary wedge systems result from scraping of clastic sedimentary sequences off a descending oceanic plate at active subduction zones (Morley et al., 2011). Sediments covering the incoming, subducting oceanic plate may be accreted at the wedge front, forming typical imbricate fans (Fig. 1A). On the other hand, buried parts of the stratigraphic sequence may underthrust the wedge body and subsequently accrete (underplate) at its base (Fig. 1A). Underplating may also occur deeper in subduction zones and concern rocks other than sedimentary rocks, i.e., during basal accretion of continental basement or slicing of oceanic crust along the subduction interface (Angiboust et al., 2014; Calvert et al., 2006; Ruh et al., 2015).

Underplating in accretionary wedges has been intensively investigated during the past four decades. A certain amount of mass transfer from downgoing plates to the base of overriding accretionary wedges was inferred by increased vertical uplift and the occurrence of metapelitic blueschist terranes landward of regions of frontal accretion (Cowan and Silling, 1978; Watkins et al., 1981). The lack of coherent reflectors from a seismic survey offshore of Mexico was interpreted as underplated pelagic and hemipelagic sediment that experienced extremely complex deformation, potentially representing a tectonic mélange (Moore et al., 1982). Physically and numerically derived material fluxes of underplating of such mechanically homogeneous weak mélanges depend on the occurrence of a rigid backstop and show typical return flow or corner flow patterns (Cloos, 1982; Cowan and Silling, 1978). On the other hand, layered seismic reflectors across the eastern Aleutian arc-trench corridor (Alaska, USA) indicate underplating in the form of tectonic antiformal stacking of coherent thrust sheets (Moore et al., 1991) leading to more complex material fluxes (Ellis et al., 1999).

In both of the above-mentioned cases, tectonic underplating requires a step-down of the major décollement horizon that allows underthrust material to accrete to the upper plate (e.g., Kimura et al., 1996). This shift in subduction interface localization within an accretionary wedge may be a result of weaker lithologies deeper in the section, dynamic variations in pore-fluid pressure, mechanical forcing by a rigid buttress, or thinning of the subduction channel toward the hinterland of subduction zones (Cloos and Shreve, 1988a, 1988b; Ernst, 2005; Menant et al., 2019; Meneghini et al., 2009; Shreve and Cloos, 1986). Long-term tectonic underplating of strata leads to thickening of the proximal part of a wedge and uplifts overlying sequences. When uplifted rocks get exposed, erosion removes material from the top and thins the affected part. Furthermore, underplating-related steepening of the overall surface taper may eventually result in gravitational collapse and extensional tectonics above the affected region (Menant et al., 2020; Platt, 1986; Ruh, 2017). In this sense, both underplating and surface removal of material contribute to the overall wedge taper, which mechanically maintains a value between its critical minimum and maximum (Dahlen and Suppe, 1988). Depending on the relationship between rate of basal accretion (ġ) and rate of material removal (k̇) at the rear of wedges, three different flow fields have been described, resulting in different orientations of the principal stress field (Feehan and Brandon, 1999; Ring and Brandon, 2008): (1) a thinning flow field and vertical contraction occurs when k̇ < ġ; (2) a thickening flow field accompanied by horizontal contraction results from k̇ > ġ; and (3) a mixed flow field with horizontal contraction deeper in the wedge and vertical contraction at shallow levels usually occurs where underplating is the controlling factor and surface extension a direct result of it (Platt, 1987). Ongoing underplating, antiformal stacking, rear extension, and potential surface erosion allow the exhumation of formerly underthrust material (Glodny et al., 2005; Ring and Brandon, 2008; Willner, 2005). If an accretionary wedge is large enough, such exhumed material may have undergone high-pressure, low-temperature metamorphism (Platt, 1986, 1987). Figure 1B illustrates a simplified schematic clockwise pressure-temperature (P-T) path for material following a typical particle path of shallow-dipping underthrusting and subvertical uplift (e.g., Likhanov, 2020). However, P-T paths of deeply subducted rocks commonly indicate more complicated multiphase burial and exhumation histories (e.g., Angiboust et al., 2018; Sizova et al., 2019). Sediments accreted at the wedge front may also experience elevated pressures and temperatures when further incorporated in the wedge body, but they show a distinct difference in the relative timing between thrusting and peak metamorphism (van Gool and Cawood, 1994): basal accretion occurs at peak metamorphic conditions, while peak metamorphism postdates thrust imbrication related to frontal accretion (Fig. 1B).

Analogue and numerical experiments of wedge systems demonstrate that the implementation of several potential décollement layers leads to underplating at the rear of the wedge if the basal one is at least as strong as the upper one (Bonnet et al., 2007; Kukowski et al., 2002; Ruh et al., 2012; Stockmal et al., 2007). Underplating may also be enhanced by variable plate coupling along the shallow subduction interface (Menant et al., 2019) or forced by the presence of a rigid backstop (Kukowski et al., 2002). Furthermore, surface processes, i.e., erosion and sedimentation, play a key role in the morphological and structural evolution of such systems, and depending on their rates, complex particle paths may be tracked (Bonnet et al., 2008; Konstantinovskaya and Malavieille, 2011). Décollement strength, mechanical forcing in form of a rigid backstop, and intensity of surface erosion all affect the dynamics of underplating in accretionary wedge systems. However, most existing studies only focused on the effect of either surface processes (e.g., Bonnet et al., 2007), several décollements (e.g., Stockmal et al., 2007), décollement strength (e.g., Kukowski et al., 2002), or backstop geometry (e.g., Rossetti et al., 2002). Therefore, a general study is needed that incorporates different features affecting underplating to assess their dynamic interaction and the relative impact of each process.

In the present study, mechanical numerical experiments are conducted to determine and quantify the main features enhancing or hampering tectonic underplating and subsequent uplift in accretionary wedge settings. A further goal is to evaluate and understand particle paths tracked in these experiments. The numerical model includes a downgoing elastic plate, an overlying sedimentary sequence, and a rigid backstop resulting in sediment accretion. In contrast to earlier studies, different parameters such as décollement strength, surface erosion, elastic stiffness of the downgoing plate, and geometry of the rigid backstop are simultaneously varied to test their relative effect on the tectonic evolution of such systems. Resulting wedge geometries and particle paths are compared to results of previous modeling studies and to several natural examples for a better understanding of the importance of underplating for the structural evolution of accretionary wedges.

Underplating is reported to occur in many subduction zones around the world. In a more general form, it refers to all material that is detached from the downgoing plate and accreted to the overriding plate related to a step-down of the subduction interface. In this sense, deep underplating (>30 km depth; e.g., Angiboust et al., 2018; Escuder-Viruete et al., 2011) is responsible for the majority of exhumed low-temperature high-pressure terranes (Agard et al., 2018; Guillot et al., 2009; Platt, 1993), alongside other viable processes such as serpentinite diapirism (Maekawa et al., 1993). However, here I focus on underplating in accretionary wedges at shallow depths without the implication of involvement of either oceanic crust or continental basement in the deformation process. In the following paragraph, several examples of accretionary wedges exhibiting important underplating are introduced.

The Kodiak accretionary complex (∼200 km wide) is the result of northwestward subduction along the Aleutian trench that has been active since the Early Jurassic (Sample and Fisher, 1986). From northwest to southeast, it can be divided into a Jurassic metamorphic terrane, an Upper Cretaceous mélange, Upper Cretaceous turbidites, a Paleocene volcanic and turbiditic unit, and Eocene turbidites (Moore et al., 1983). Seismic reflection data reveal layered reflectors beneath the Mesozoic core of the accretionary complex at depths of ∼9–35 km (Moore et al., 1991). Those reflectors resemble typical nappe structures and were interpreted to represent underplating that occurred during the last 50 m.y. as a response of lateral wedge growth to maintain the overall critical taper (Fig. 1C).

The Barbados Ridge accretionary complex is associated with the Lesser Antilles subduction zone, where the North and South American plates have been subducted toward the west under the Caribbean plate since ca. 45 Ma (Deville and Mascle, 2012; Westbrook, 1982). South of the Tiburon Rise (15°–16°N), the complex measures >250 km across strike and exhibits a seaward-dipping backstop, a wedge body of thrust imbricates, and several kilometers thick forearc deposits on top of the overriding crust (Fig. 1D). Seismic reflection data indicate that parts of the incoming sedimentary sequence underthrust the accretionary wedge toe (Peter and Westbrook, 1976; Speed and Larue, 1982). Underplating and duplexing near the backstop (Fig. 1D) may be related to increased sediment input from the Orinoco Delta on the Atlantic seafloor during late Miocene to Pliocene (Deville et al., 2015).

The active Nankai Trough accretionary prism describes the frontal part of the larger Japan forearc region. Along the Muroto transect, it consists of Pliocene to recent accreted terrigenous sediments scraped off the Philippine Sea plate subducting toward the northwest beneath the Eurasian plate (Moore et al., 2001). Underthrusting of an ∼150-m-thick section, observed from seismic reflection data, is explained by different mechanical properties due to undercompaction in contrast to the overlaying, imbricate wedge body. A step-down of the décollement layer occurs 25 km landward from the wedge toe, enhancing duplexing of the underthrust sequence (Bangs et al., 2004). Underplating between 25 and 40 km from the toe is accompanied by a local increase of surface slope and activity of large, deep-rooted out-of-sequence faults (Fig. 1E).

The Makran (Iran and Pakistan) is one of the largest active accretionary complexes, with a width of ∼350 km from rear to toe (Burg, 2018). It results from the northward subduction of the Oman oceanic plate beneath Eurasia since the Late Cretaceous. A large part of the wedge emerges above sea level. In the very north, ophiolites of an ancient Jurassic ocean (Hunziker et al., 2015) and imbricate slivers of plutonic rocks and deep-sea sediments form a hinterland-dipping backstop (Fig. 1F). The low total taper (∼3.5°) indicates the presence of a weak active décollement, potentially resulting from increased fluid overpressure (Davis et al., 1983). Large-scale underplating has taken place since at least the mid-Miocene, based on observed sedimentary stratigraphic contacts (Platt et al., 1985). It is suggested that ∼50% of the incoming ∼7-km-thick sedimentary sequence is underthrust and stacked further to the rear of the wedge (Fig. 1F).

The presented experiments are conducted with a finite-difference numerical code using the marker-in-cell method to track material properties and large strain (see Ruh, 2017). The numerical code applies a visco-elasto-plastic/brittle rheology. The simplified temperature field is based on a non-dynamic vertical thermal gradient of 30 °C/km throughout the entire wedge, which is a good average for slow accretionary wedges (see Underwood et al., 1993, and references therein).

The Eulerian box of all experiments measures 150 × 25 km in the x- and y-directions, respectively, with 1501 × 251 nodal points, resulting in a grid resolution of 100 × 100 m (Fig. 2). Every nodal cell contains nine Lagrangian markers that define the initial material distribution and advect according to the Eulerian velocity field. From bottom to top, markers initially prescribe a 10-km-thick oceanic crust that acts as a mechanically passive body allowing for basal flexure with a fixed viscosity of 10²³ Pa·s, not allowing any internal deformation. The oceanic crust is overlain by a 5-km-thick sedimentary sequence with an initial viscosity of 10²⁴ Pa·s, a friction angle of 30°, and a cohesion of 5 MPa. The applied fluid pressure ratio in the sedimentary sequence is 0.7, except for two 300-m-thick décollement layers at its very bottom and 2 km above its bottom. Fluid pressure ratios for the basal (λ_b) and intermediate décollement (λ_i) vary between 0.85 and 0.95 for different experiments (Table 1). On the left side of the model domain, a backstop overlays the basal décollement layer. Backstop geometries vary between dipping toward the toe, dipping toward the hinterland, or vertical (Fig. 2). Backstops have an initial viscosity of 10²⁴ Pa·s and are stronger than the sedimentary sequence, with a friction angle of 30°, a cohesion of 20 MPa, and no fluid pressure (λ = 0). Between x = 10 km and x = 50 km, slope sediments smooth the topographic change (of 5 km) from backstop to sedimentary sequence. These slope sediments exhibit the same rheological parameters as the sedimentary sequence. All rock types have the same density of 2500 kg/m³. Above the sedimentary sequence lays a 3-km-deep water body with a density of 1000 kg/m³. The top 7 km of the Eulerian domain is defined by a density of 1 kg/m³. Together with the water body, these low-viscosity layers (η = 10¹⁷ Pa·s) act as “sticky air”, assuring low shear stresses along the rock-air interface (Crameri et al., 2012).

Lateral velocity boundary conditions resemble those of typical sandbox models, where the lower plate is pulled out below a rigid backstop. However, the implementation of the elastic flexure complicates this otherwise simple setup. The right side boundary prescribes an incoming velocity of –1 cm/yr in the horizontal direction and free slip in the vertical direction (Fig. 2). The left side has no cross-boundary velocity and no slip condition above the basal décollement, while the oceanic plate exits the model domain in a direction parallel to the elastic flexure at x = 0 with a total velocity of –1 cm/yr. This description of the left boundary prohibits sediments from exiting the model domain. The bottom boundary condition defines a total left- and downward-directed velocity of –1 cm/yr parallel to the overlaying elastic flexure (Fig. 2). This assures that the volume of oceanic crust entering and exiting the model domain is equal and that no major deformation takes place within this oceanic crust. The top boundary prescribes free slip in the horizontal direction and upward vertical velocity preserving the conservation of mass (–0.1 cm/yr).

The applied velocity boundary conditions imply that all material below the lower décollement (oceanic crust) exits the model domain, while all material above the lower décollement is accreted and incorporated into the accretionary wedge. This ignores the potential subduction of sediments observed in large-scale experiments without lateral boundary constraints (e.g., Currie et al., 2007; Gerya and Meilick, 2011).

In the following, I present eleven experiments that were conducted to test the effects of backstop geometry, flexural rigidity of the oceanic lithosphere, décollement strength, and surface erosion on the structural evolution of accretionary wedge systems involving underplating. A twelfth experiment investigates the influence of a sudden olistostrome emplacement on the system. All experiments deviate in one or two above-mentioned parameters from the reference experiment M1 (see Table 1).

Reference experiment M1 is defined by a backstop dipping toward the wedge toe, a basal flexural rigidity of D = 3 × 10²³ Nm², décollement fluid pressure ratios of λ_b = 0.9 and λ_i = 0.95, and an erosion diffusion coefficient of κ_ero = 10⁻⁶ m²/s (Fig. 3). After the initiation of shortening, deformation within the sedimentary sequence localizes close to the backstop. At 3 m.y., the upper part of the sequence forms thrust imbricates with piggyback basins in front of where the slope sediment pinches out (wedge toe), while the lower sequence is scraped off the basal décollement along its contact with the backstop (Fig. 3A). The vertical thickening related to underplating pushes up the upper sequence and slope sediments, which reach water level. Consequently, the forearc-type basin toward the hinterland of this popup zone is filled with sediments. The weight of the thickened sequences and the forearc sediments results in a deflection of the wedge base of ∼4 km at the left boundary (Fig. 3A). With ongoing shortening, the wedge toe defined by the frontal thrust within the upper sedimentary sequence propagates away from the backstop, resulting in horizontal wedge growth (Fig. 3B). This is accompanied by further basal accretion and stacking of underplated material. Ongoing stacking and vertical growth results in further deflection of the base and subsidence of the overfilling forearc basin. Toward the hinterland, the stacked pile of underplated material migrates upwards along a steep, toeward-dipping backthrust reaching from the surface down into the lower décollement, separating the undeformed forearc sediments from the popup structure (Fig. 3B). After 9 m.y., an ∼100-km-wide wedge is established with a lower sedimentary sequence that underthrusts the wedge for ∼50 km before it gets underplated (Fig. 3C). Earlier underplated material gets uplifted by ongoing nappe stacking at the wedge base. Most of the vertical wedge thickening is accommodated along the main backthrust that separates the now ∼5-km-thick forearc sediments from the popup structure above the underplated region with highest topographic elevation (Fig. 3C). After 12 m.y. of shortening, initially underplated slivers are almost exhumed to the surface due to ongoing uplift and subaerial erosion (Fig. 3D). The accretionary wedge of ∼15 km vertical thickness can be divided into three structural sections, from left to right: (1) forearc basin with a thickness of several kilometers of sedimentary strata overlaying the backstop/overriding crust without major internal deformation; (2) an ∼50-km-wide zone of intense underplating accompanied by vertical uplift and surface topography; and (3) an ∼80-km-wide forward-propagating imbricate thrust system detached along the intermediate décollement layer with a stable surface taper (Fig. 3D).

Pressure-temperature (P-T) paths of four clusters (each 1 km²) of Lagrangian markers are tracked throughout the shortening process (Figs. 3E–3H). Pressure here denotes tectonic pressure (mean stress), and temperature is defined using depth below the surface and a thermal gradient of 30 °C/km. The applied thermal gradient is relatively warm, but because it is purely dependent on depth, the form of the resulting P-T paths is not affected. The effects of arbitrary geothermal gradients can be tested by varying the temperature axis in the P-T diagram. The four clusters are selected to illustrate an entire P-T cycle, including before entering the wedge, underthrusting, basal accretion, antiformal stacking, and uplift (Figs. 3A–3D). The leftmost cluster is immediately affected by underplating and then migrates upward without undergoing any further visible deformation. Its P-T path shows initial temperature increase due to burial and simultaneous pressure drop (Fig. 3E). Further uplift is characterized by a linear pressure and temperature decrease. The second cluster from the left initially undergoes a typical prograde P-T path increasing in pressure and temperature (Fig. 3F). An isothermal pressure peak occurs at 9 m.y. (blue point) before a retrograde path establishes with decreasing pressure and temperature. The third cluster from the left initially underthrusts the evolving wedge, which results in pressure and temperature increase due to thickening of the overlaying upper sequence (Fig. 3G). After the peak is reached, a reduction in pressure at isothermal conditions takes place. The rightmost cluster shows a similar P-T path, with peak pressure and temperature of 250 MPa and 250 °C, respectively, occurring at 9 m.y. when it detached along the basal décollement (Fig. 3H). Further incorporation into the underplated stack is accompanied by a decrease in tectonic pressure.

Rock composition and the second invariant of the strain-rate tensor are illustrated for two experiments (M2 and M3) after 12 m.y. of shortening, only varying in backstop geometry from the reference experiment M1 (Fig. 4; Table 1). Strain rates of the reference experiment show a steep active backthrust along which underplated material migrates upward (Fig. 4A). In the top 3 km, this backthrust deflects into a shallower backthrust and a very shallow forethrust along the uplifted intermediate décollement. Between x = 30 km and x = 80 km, active underplating takes place, whereas stacked material above is uplifted along the backthrust without any further internal deformation (Fig. 4A). Imbrication, i.e. thrusting, of the upper sequence occurs between x = 80 km and x = 140 km with local out-of-sequence thrusting further toward the backstop.

Experiment M2, with a vertical backstop geometry, exhibits a backthrust that roots at ∼5 km below the the surface, where it splits into a diffuse fault pattern (Fig. 4B). Underplating occurs between x = 10 km and x = 70 km. Underplated nappe stacks bounded by active thrusts contain several imbricates with little internal deformation. Rock composition illustrates the occurrence of a large-offset out-of-sequence thrust cropping out roughly where the initial slope sediments pinch out (x ≈ 30 km), pushed by the underplated, uplifted material (Fig. 4B).

If the backstop initially dips toward the hinterland (Fig. 2), no major backthrust develops to facilitate uplift of underplated material (Fig. 4C; experiment M3). Rock composition illustrates that underplated material is pushed below the backstop, which then gets bent upward (counterclockwise rotation). However, strain rates do not indicate extensive plastic-brittle deformation of the backstop. Similarly to the above-described experiments, several underplated thrust sheets are stacked together to build nappes that are separated by active thrusts of high strain rates (Fig. 4C). A major out-of-sequence thrust reaching from the surface down to the basal décollement occurs at the base of the backstop, overthrusting the underplated, antiformal stacks.

Varying the flexural rigidity of the base while keeping other parameters equal (experiments M4 and M5; Table 1) strongly affects the morphological evolution of accretionary wedge systems (Fig. 5). In case of a strong lower plate (D = 10²⁴ Nm²; experiment M4), the left side of the experiment subsides ∼3.5 km as a response to the additional load following vertical structural thickening of the wedge and deposition of forearc sediments (Fig. 5A). Relatively little basin-scale subsidence during wedge growth quickly leads to subaerial exposure and surface erosion. Removal of overlying material in the region of active underplating facilitates upward motion of stacked nappes along the major backthrust. The major backthrust is steeply dipping below a depth of ∼3 km but diverges into a shallow-dipping fault toward the surface, overthrusting the forearc sediment (Fig. 5A).

An intermediate flexural rigidity (D = 3 × 10²³ Nm²; reference experiment M1) results in ∼6 km of vertical deflection on the left side (Fig. 5B). The stratification of forearc sediments illustrates that subsidence is initially faster and balances out between newly deposited forearc sediments and isostatic vertical deflection. The structural development of this experiment M1 was discussed above. In contrast to experiment M4 with a stronger lower plate, more intense underplating and vertical stacking takes place away from the major backthrust, forming a buried wedge taper between x = 40 km and x = 70 km (Fig. 5B).

A relatively thin elastic plate with a rigidity of D = 10²³ Nm² (experiment M5) exhibits a vertical deflection of almost 10 km after 9 m.y. of shortening (Fig. 5C). The larger volume of forearc sediments due to subsidence enhances further vertical deflection. Two major backthrusts are visible, both rooting down to the basal décollement where the backstop pinches out. The slightly shallower-dipping backthrust does not show any offset in the uppermost 2 km of the forearc deposits, indicating that it ceased thrusting while the steeply dipping backthrust emerged to the surface (Fig. 5C).

Figure 5D illustrates surface topography and the geometry of the wedge base for the above-discussed experiments with different flexural rigidity after 9 m.y. A more rigid lower plate leads to wider wedges with higher topography. However, the difference in topographic elevation is less than the difference in vertical basal deflection, indicating more intense erosion of the rigid experiment M4 due to earlier surface exposure. Topography profiles furthermore show that the parts of the wedges affected by underplating are slightly steeper than the imbricate thrust systems detached along the intermediate décollements (Fig. 5D). For a comparison, theoretical surface tapers (Fig. 5E) for these three experiments (M1, M4, M5) are derived from the observed basal dip angles (Fig. 5F) and the internal and basal strength of the wedges (Davis et al., 1983). These analytical solutions show steeper expected surface tapers if derived using the strength of the basal décollement (solid lines) instead of the intermediate décollement (dashed lines). Surface tapers of the frontal imbricate fans retrieved from the experiments are steeper for flatter wedge bases (Figs. 5A–5C). However, they are ∼1° steeper than the analytically predicted ones for the frontal parts of the wedges (Fig. 5E, dashed lines).

Marker composition and tectonic overpressure are shown for experiments with different basal and intermediate décollement strength imposed by varying fluid pressure ratios after 12 m.y. of shortening (experiments M6 and M7; Fig. 6; Table 1). Tectonic overpressure denotes the deviation of effective (tectonic) from lithostatic pressure. If positive, a rock body is under compression, while negative values relate to extension. Experiment M6 with λ_b = 0.85 and λ_i = 0.95 shows that underplating initiates at x ≈ 60 km and immediately leads to antiformal stacking and uplift along the major backthrust (Fig. 6A). The intermediate décollement forms a toeward-dipping roof normal fault covering the duplex structure (x = 30–60 km). At the wedge front, a typical imbricate fan develops above the intermediate décollement. For décollement strengths given by λ_b = 0.9 and λ_i = 0.95 (experiment M1), the underplated region is wider (x = 30–80 km) but antiformal stacking occurs only further toward the rear of the wedge at x ≈ 45 km (Fig. 6B). If the basal and intermediate décollements exhibit the same strength given by λ_b = λ_i = 0.9 (experiment M7), imbrication of the upper sequence and underplating both occur at the wedge front (Fig. 6C). Underplating occurs at the very rear along the backthrust, enhanced by surface erosion.

Absolute values of overpressure are similar for experiments with different décollement strength (Fig. 6). In the frontal part, the largest overpressures occur along the wedge front in the upper sequence and in undeformed underthrust parts of the lower sequence. Along the rear, the largest overpressures occur where underplated material is detached from the basal décollement and forced upward along the backthrust. All three experiments show negative overpressures (extension) in the uppermost 5–7 km of the region affected by underplating as well as in form of patches along the décollement layers (Fig. 6).

Changing the diffusion coefficient for subaerial erosion strongly affects the structural appearance of numerically modeled accretionary wedges after 12 m.y. of shortening (experiments M8 and M9; Fig. 7; Table 1). Slower erosion (experiment M8; κ_ero = 10⁻⁷ m²/s) leads to increased vertical growth of the wedge, which results in ∼7.5 km of subsidence at the left side compared to the initial state (Fig. 7A). Underplated material is stacked from x = 70 km leftward and migrates upward along the main backthrust. Experiment M1 with an intermediate erosion diffusion coefficient (κ_ero = 10⁻⁶ m²/s) shows similar subsidence (∼7 km), but underplated material is exhumed closer to the surface, indicating lower topographic relief (Fig. 7B). Very intense surface erosion (experiment M9; κ_ero = 10⁻⁵ m²/s) leads to very low topographic elevation and prohibits significantly vertical (and as a result horizontal) growth of the wedge, hampering subsidence (Fig. 7C). Underplating initiates only near the backstop (x ≈ 20 km) and immediately migrates upward. In general, lesser surface erosion after the same amount of shortening implies more accreted material, hence a wider and vertically thicker wedge. This is illustrated by the location of the active frontal thrust as well as by the position where underthrust sediments get underplated (Fig. 7).

Vertical velocities show localized uplift in the form of a triangular popup structure in different parts of the wedge for all experiments (Fig. 7). Uplift velocities of structures accreting along the intermediate décollement, such as in frontal accretion, and basal accretion away from the backstop show values of ∼1 mm/yr for all experiments (note different scales for uplift velocity in Fig. 7). However, upward displacement along the major backthrusts differs, from ∼1 mm/yr for experiment M8 to ∼2 mm/yr for experiment M1 and ∼4 mm/yr for experiment M9 (Fig. 7). Along the rear of the wedge, major backthrusts reach from the surface down to the basal décollement, while their conjugate landward-dipping forethrusts migrate through the uplifted material and merge into the roof faults of the underplated antiformal stack. Experiments M1 and M8, with intermediate and little erosion, respectively, show little activity of the intermediate décollement above the underplated thrust sheet package (Figs. 7A, 7B), whereas the intermediate décollement is reactivated as a backthrust if erosion is intense (experiment M9; Fig. 7C).

Figure 8 illustrates marker distribution and vertical velocity of three additional experiments with a hinterland-dipping backstop varying from experiment M3 (Fig. 4c) in erosion, basal flexure, and additional sedimentary load in form of an olistostrome (experiments M10–M12; Table 1). After 11 m.y. of shortening, experiment M10, with intense erosion (κ_ero = 10⁻⁵ m²/s), displays an ∼50-km-wide region of active underplating, with a major out-of-sequence thrust toward the rear of the underplated sequence that roots into the basal décollement (Fig. 8A). Localized uplift of underplated and overlying strata is illustrated by vertical velocity reaching upward velocities of ∼2 mm/yr (Fig. 8A; upward velocities are denoted by negative values). Uplift at the rear occurs in form of a popup structure and is significantly faster than at the wedge toe, similar to experiment M9 with intense erosion and a toeward-dipping backstop (Fig. 7C). However, the experiment with a toeward-dipping backstop develops a much narrower area of underplating and almost twice as fast uplift rates compared to the hinterland-dipping backstop experiment (Figs. 7C, 8A).

Experiment M11, with a rigid elastic flexure (D = 10²⁴ Nm²), results in a wedge similarly thick as for intense erosion (Fig. 8B). However, underplating is more distributed (∼70 km wide) and there is no major fault crosscutting the hinterland-dipping backstop. The overall morphology and structure of the wedge, except for weight-related subsidence, is comparable to those of experiment M3 (Fig. 4C). Vertical velocities show uplift rates of ∼1 mm/yr at the rear, the wedge toe, and the front of underplating (Fig. 8B).

Experiment M12 illustrates an extreme case with strong flexural rigidity (D = 10²⁴ Nm²) and intense surface erosion (κ_ero = 10⁻⁵ m²/s) combined (Fig. 8C). Furthermore, the trench is filled with a low-viscosity (10¹⁹ Pa·s), low-density (2200 kg/m³) olistostrome after 2 m.y. of shortening (all “water” markers below y = 8 km are converted to “olistostrome” markers). After 9 m.y. of shortening, underplating occurs along a horizontal distance of ∼30 km. An out-of-sequence thrust rooting to the lower décollement is active along the rear side of the underplated area (Fig. 8C). From x = 40 km to x = 80 km, the olistostrome body covers earlier active thrust sheets, whereas a new imbricate fan developed at the wedge toe at x = 80–120 km. Several small growth synclines subsided into the mobile olistostrome. Vertical velocity shows that upward motion along the rear of ∼3 mm/yr is reached in form of a popup structure (Fig. 8C). In contrast to experiments M9 and 10, a major out-of-sequence thrust develops ∼20–30 km away from the backstop, incorporating several underplated thrust sheets into its hanging wall. The forethrust of the popup structure reaches the surface where the olistostrome starts and roots to the basal décollement (Fig. 8C).

All presented numerical accretionary-wedge experiments exhibit frontal accretion and forward propagation of deformation in the upper part of the incoming sedimentary sequence, while the lower part gets underplated. However, the results show that geometrical and frictional properties of the involved rocks and acting surface processes strongly influence where and to what extent underplating of underthrust material occurs, how it is incorporated into the wedge, and how it affects overall wedge dynamics. Several key observations on the structural evolution of accretionary wedges can be made depending on input parameters (Table 1):

• 1. Only experiment M7, where the basal and intermediate décollement exhibit equal brittle-plastic strength, shows underplating to occur at the wedge toe in form of a frontal thrust (Fig. 6C). This is because stress levels transferred across the intermediate décollement (or any fault zone) are capped by its brittle strength (see Ruh et al., 2012, their figure 10). In experiment M7, the intermediate décollement transfers exactly enough stress to fail the basal décollement as well. In contrast, weaker basal décollements in all other experiments lead to underthrusting at the front and underplating forced by indentation of the backstop.

• 2. The backstop geometry affects the shape of stacked underplated thrust slivers. While a toeward-dipping backstop remains passive and leads to a highly effective backthrust facilitating uplift of underplated material (Fig. 4A), a hinterland-dipping backstop is incorporated into the wedge and a major hinterland-dipping fault establishes (Fig. 4C). This fault acts as an out-of-sequence thrust hampering exhumation of underplated rocks, as suggested for the Franciscan Complex (California, USA) (Platt, 1993), and not as a large normal fault, which could enhance vertical uplift as proposed for the Himalayas (Royden and Burchfiel, 1987) or Corsica (Fournier et al., 1991).

• 3. Localized underplating, antiformal stacking, and related exhumation are most effective in experiments with high flexural rigidity (Fig. 5A) or intense surface erosion (Fig. 7C). In these experiments, the rear part of the wedges is kept relatively thin either by removal of material at the top, hampering subsidence, or by little subsidence, leading to early surface exposure and erosion (Perrin et al., 2013). Resulting flatter basal décollements lead to steeper surface slopes (Figs. 5E, 5F), enhancing exhumation at the rear (Davis et al., 1983). A steeper minimum critical taper in contrast to the reference experiment can also be achieved by a stronger basal décollement, which leads to more localized underplating and enhanced uplift and exhumation, although less effectively than in the examples above (Fig. 6A).

Pressure-temperature (P-T) paths of marker clusters covering different parts of the accretionary wedge illustrate a steeper prograde path followed by isothermal pressure loss and a flatter retrograde path (Figs. 3E–3H). The shape of this P-T path resembles that of published numerical and natural data from exhumed rocks within accretionary wedges and compressional orogens (Fig. 1B; e.g., Giunchi and Ricard, 1999; Platt, 1986; van Gool and Cawood, 1994). In the numerical experiments here, P denotes dynamic pressure and T is calculated using the depth below surface and a gradient of 30 °C/km. This implies that the steepness of the pro- and retrograde paths is solely defined by the deviation of effective from lithostatic pressure (Mancktelow, 2008). Experiments show that positive overpressure occurs during underthrusting and basal accretion, while negative overpressure (or underpressure) develops at shallow levels of regions affected by underplating (Fig. 6), explaining the observed P-T path. This switch of principle stress orientation describes a mixed flow field, where rates of basal accretion and surface removal are roughly in equilibrium (Ring and Brandon, 2008). P-T paths showing isothermal pressure loss reported from various mountain belts have been related to fast initial exhumation rates (Ring et al., 1999, and references therein). However, the experiments presented here indicate that this feature can alternatively be related to a switch from an over- to underpressured tectonic regime and does not necessarily require fast exhumation rates.

In this section, the temperature-time (T-t) evolution of underplated and exhumed tectonic slivers from different experiments is discussed (Fig. 9). Natural examples of accretionary wedges show a wide distribution of geothermal gradients between ∼15 and 60 °C/km (e.g., He and Zhou, 2019; Underwood et al., 1993; Yamano et al., 1982). Thermal gradients depend on sediment compaction allowing or hindering fluid flow and consequent heat convection (Le Pichon et al., 1991) and can vary strongly between the toe of a wedge to its interior (Ashi and Taira, 1993). In the experiments presented here, every marker’s temperature was calculated from its vertical depth below the surface line and a geothermal gradient of 30 °C/km that serves as an average of reported values. This static geothermal gradient is a simplification of potentially strongly varying gradients across natural systems (see the Limitations of the Modeling Approach section below). T-t paths of underplated and now-exposed markers (1 km²) in the experiments give a measure of their exhumation rate and can ultimately be compared to thermochronological data from natural accretionary wedges. Underplated and uplifted material of the reference experiment M1 shows an increase in temperature during the initial 3 m.y. related to the thickening of the overlying upper sequence (Fig. 9A). Maximum temperatures of ∼140 °C indicate that uplift initiated early and the material was never deeply buried. During the last ∼2.5 m.y., exhumation accelerated, probably due to activation of a shallow-dipping backthrust and related extension (Fig. 4A). In contrast, T-t curves for experiment M4, with a strong flexural rigidity, show higher peak temperatures (Fig. 9A), although overall the wedge remains smaller (Fig. 5A). The important feature here is the fast heating and cooling related to early initiation of erosion processes that underthrust material and then immediately exhumed it. A stronger lower décollement (experiment M6) also enhances uplift, but its rates remain similar to those of the reference experiment (Fig. 9B). Higher maximum temperatures of ∼180 °C indicate that the wedge thickened before the respective material was detached. The fastest exhumation rates are reported for experiment M9, with a high surface erosion (Figs. 7C, 9C). In that experiment, temperatures reach up to 300 °C and cool down to zero within 4–5 m.y., almost reaching the extreme cooling rates of 150–350 °C/m.y. reported from the Betic Cordilleras (Spain) related to tectonic unroofing in an extensional environment (Zeck et al., 1992).

The initial geometry of the backstop affects the upward migration of underplated material (Fig. 4). This is also true for high surface erosion, where a hinterland-dipping backstop leads to slightly slower exhumation rates (Fig. 9D). However, whereas a toeward-dipping backstop results in a major backthrust along which localized uplift occurs (Fig. 5C), experiment M10 with hinterland-dipping backstop shows a more complex structural framework (Fig. 9E). As the major out-of-sequence thrust initially follows the backstop-sediments interface (Fig. 4C), it may migrate forward and localize between underplated thrust nappes. Finally, a convex backthrust cutting the backstop forms, facilitating further uplift (Fig. 9E). Such a framework may explain local variation in observed cooling rates in ancient amalgamated cores of thrust belts (e.g., DeCelles et al., 2020; Mailhé et al., 1986).

The morphological and structural evolution of accretionary wedge systems has been investigated by a wealth of analogue and numerical modeling studies (see Biuter, 2012, and Graveleau, 2012, for reviews). Furthermore, the mechanics of wedges has been analytically described based on the balance of forces within such systems (Dahlen, 1984; Dahlen et al., 1984; Davis et al., 1983). In the following, the numerical experiments presented here are compared to previous studies that investigated tectonic underplating.

Several analogue modeling studies have tested the effect of multiple décollement layers in the deformed sand pile that ultimately leads to ramp-flat geometries and potentially to antiformal stacking (e.g., Bonini, 2001; Bonnet et al., 2007, 2008; Hoth et al., 2007; Konstantinovskaya and Malavieille, 2011; Kukowski et al., 2002; Mugnier et al., 1997; Turrini et al., 2001). Kukowski et al. (2002) demonstrated that decoupling within the stratigraphic pile followed by underplating and the development of duplexes is a result of high basal friction and the presence of weak intermediate décollement layers. In their experiments, underplating is forced by the presence of a rigid end wall, and the developing antiformal stacks uplift the overlying deformable backstop, similar to numerical experiments of this study with a hinterland-dipping backstop (Fig. 4C).

Malavieille et al. (2016) investigated parameters responsible for tectono-sedimentary reworking of oceanic crustal fragments within an accretionary wedge. They applied a very shallow toeward-dipping rigid backstop and introduced variable basal décollement strength and an intermediate décollement at different positions in the stratigraphic sequence in two of their experiments (Malavieille et al., 2016, their experiments 2 and 5). Both experiments produced stacks of underplated thrust sheets that overthrust the backstop along a major backthrust, similar to experiments with a toeward-dipping backstop presented in this study (Fig. 4A). Furthermore, Malavieille et al. (2016) reported similar effects for experiments with different décollement strengths. A wide region of underplating and an overall constant surface taper develop if both basal and intermediate décollements are weak (their experiment 2), comparable to experiment M7 in this study, with equal strength for both décollements (Fig. 6C). For a significantly weaker basal décollement, underplating is more localized along the backstop, and the surface slope can be divided into a steeper rear part and a flatter frontal part (their experiment 5), comparable to experiment M6 of this study, with an increased strength difference between the décollements (Fig. 6A).

Bonnet et al. (2007, 2008), Malavieille and Konstantinovskaya (2010), and Konstantinovskaya and Malavieille (2011) investigated the impact of surface processes on the structural evolution of accretionary wedges in general and underplating in particular. Experiments without erosion but two décollement layers develop a relatively wide area of underplating, where uplift of the antiformal stack is hampered by the overlying sequence. Additional surface erosion leads to localized underplating and faster uplift rates (Malavieille and Konstantinovskaya, 2010, their figure 3d). Intense antiformal stacking is enhanced by increased erosion and a larger strength contrast between the basal and intermediate décollement levels (Konstantinovskaya and Malavieille, 2011). Both effects are observed in numerical experiments here: (1) a weaker intermediate décollement leads to a narrow area of underplating (Fig. 6A), and (2) intense surface erosion results in localized antiformal stacking and fast uplift (Fig. 7C). Bonnet et al. (2007, 2008) presented similar results showing that a certain amount of erosion is needed to facilitate vertical exhumation of underplated material forced by a strong backstop. Gutscher et al. (1998) demonstrated that underplating may occur in order to reequilibrate a critical taper. Frontal accretion develops for a very steep initial taper in order to lower it, while underthrusting of elongated slivers occurs for shallower initial tapers. Underthrust slivers get underplated when they reach the rigid back wall and upward motion of material is accompanied by a deeply rooted backthrust (Gutscher et al., 1998, their figure 4). More sophisticated analogue experiments were presented by Lohrmann et al. (2006), in which they investigated the structural evolution of accretionary wedges affected by subduction accretion or erosion. Similar to other experiments, they showed that multiple décollements lead to separated frontal and basal accretion (Lohrmann et al., 2006, their figure 11.4). Their experiments furthermore illustrate that intense underthrusting occurs along erosive subduction channels, but without underplating. All accretive experiments developed a major backthrust related to the location of underplating, similar to experiments presented in this study (e.g., Fig. 4A).

Since the work of Cloos (1982), who calculated possible material fluxes related to return flow in large-scale accretionary wedges, many numerical modeling studies have investigated the structural development of underplating in such systems. A triad of studies tested the modes of accretion along convergent margins by means of large-scale numerical experiments including lithospheric bending and possible sediment subduction (Beaumont et al., 1999; Ellis et al., 1999; Pfiffner et al., 2000). They were able to establish a dynamic environment with an active subduction interface and an overlying conduit region where tectonically underplated material and preexisting lithospheric slivers behave episodically stationary, deform slowly, and eventually subduct into the mantle or attach to the overriding plate. Significant amounts of underplating in the accretionary wedge are achieved for experiments with lower sediment density, leading to buoyancy of underthrust material (Beaumont et al., 1999, their experiment 5), and for subduction zones with slab retreat (Beaumont et al., 1999, their experiment 10). Furthermore, spatial variation in material strength and density may lead to episodic basal accretion of small crustal terranes, resulting in crustal-scale fold nappes and related stacking of underthrust material at the base of the overriding plate (Ellis et al., 1999, their experiment 36).

Selzer et al. (2007), Stockmal et al. (2007), and Ruh et al. (2012) investigated the influence of intermediate weak décollements on the evolution of “sandbox”-type numerical experiments. Experiment IC1 of Selzer et al. (2007) exhibited two décollements with equal strength. Underplating occurs across the entire wedge, with more localized antiformal stacking and uplift along a major backthrust. This is in agreement with experiment M7 with equal décollement strength presented here, where frontal accretion of both the upper and lower sequence occurs at the wedge toe (Fig. 6C). Experiments of Stockmal et al. (2007) showed the same effect for equal-strength décollements. However, for weaker basal décollements, major underplating and antiformal stacking occurred, further enhanced by slope-dependent erosion and fast deposition of developing piggyback basins. This is consistent with results of this study, illustrating that intense erosion and increased strength difference between the décollements triggers antiformal stacking along the backstop (Figs. 6A, 7C). Ruh et al. (2012) furthermore illustrated that underplating forced by a strong basal and a weak intermediate décollement results in variable surface slopes: a steeper slope at the rear that depends on the lower décollement strength, and a frontal shallow slope resulting from the upper décollement strength (their experiment DF3). The same effect can be observed in experiment M6 here, where the underplated region shows steeper surface slopes (Fig. 6A).

Finally, Menant et al. (2019, 2020) presented mantle-scale thermo-mechanical experiments focusing on the dynamics along the subduction interface and deep basal accretion. They argued that periodic underplating and the formation of thick duplex structures develop related to fluid-driven stress changes along the subduction interface beneath the accretionary wedge (Menant et al., 2019). Although experiments of the present study lack the implementation of fluid flow and related variation of décollement strength, some first-order comparisons can be drawn: (1) a main locus of underplating is related to the strong overriding backstop; (2) the shallow part of the underplated region generates tectonic underpressure and undergoes normal faulting, similar to experiments presented here (Fig. 6); and (3) ongoing underplating may lead to the formation of a major backthrust crosscutting the hinterland-dipping backstop (Menant et al., 2019, their figure 4), comparable to experiments presented here with similar backstop geometries (Fig. 8).

Numerical results show that the structural evolution of accretionary wedges depends on the investigated parameters. However, there are many more factors playing a role in potential underplating that were ignored in this study and may vary within a wedge system from toe to rear and/or from surface to base, such as variable fluid flow (Saffer and Tobin, 2011), metamorphic fluid production (Spinelli and Saffer, 2004), or carbonate cementation (Sample, 1990), all affecting pore fluid pressure and thus brittle strength (see Ruh, 2017). Furthermore, the strength of the subduction interface potentially decreases significantly at shallow depths (10–15 km) due to the activation of pressure-solution processes (Fagereng and den Hartog, 2017). Nevertheless, some first-order observations from numerical experiments can be compared to natural examples of accretionary wedges in terms of their structural framework.

Figure 1C illustrates an interpretative profile of the Kodiak accretionary wedge in the early Tertiary, when slivers of underthrust sediments (Fisher and Byrne, 1987; Moore, 1978) were initially stacked by basal accretion (Moore et al., 1991). Since the Eocene, the Kodiak accretionary wedge has undergone substantial lateral growth due to rapid frontal accretion, accompanied by accelerated uplift of underplated material indicated by nonmarine Oligocene conglomerates overlying late Eocene submarine fan deposits (Moore et al., 1991). Furthermore, similar apatite and zircon fission-track ages (44 Ma) of underplated and now-exhumed material correspond to very fast cooling rates in a regional-scale antiform (Clendenen et al., 2003). This demonstrates how fast sediment accumulation that flattens the wedge taper triggers uplift (and potentially underplating) at the rear part to maintain the critical taper (e.g., Simpson, 2010). The state of the Kodiak wedge before intense sediment accumulation (Fig. 1C) compares to numerical experiments with little erosion and sedimentation (Fig. 7A). This experiment (M8) shows a wide area of underplating (∼50 km) that is structurally covered by younger sequences. Due to slow erosion rates, the uplift velocity and thus cooling rates are relatively low (Fig. 7A). The post-Eocene dynamics of the Kodiak accretionary wedge can to some extent be compared to experiment M12, which implicates an instantaneous emplacement of a sedimentary body together with increased erosion intensity (Fig. 8C). Instantaneous emplacement might also represent sediment accumulation rates much faster than tectonic rates, which then define the overall surface slope of a growing wedge. Numerical results show that the flattening by surface processes (fast sedimentation in the front, fast erosion at the rear) imposes accelerated uplift rates of underplated material (Fig. 8C).

The overriding oceanic crust of the Lesser Antilles subduction zone forms a toeward-dipping backstop for the Barbados Ridge accretionary wedge (Fig. 1D; Torrini and Speed, 1989). Similar to numerical experiments (Fig. 3A), a major backthrust may facilitate uplift of the outer ridge and basin strata above the overriding crust (Silver and Reed, 1988). No significant erosion and the deposition of as much as 10 km of forearc deposits (Fig. 1D) resemble experiment M5, with a low flexural rigidity, where subsidence mainly compensates for vertical wedge growth (Fig. 5C). A relatively low flexural rigidity (D = 10²³ Nm²) is in agreement with the Late Cretaceous age of the incoming oceanic crust, resulting in an effective elastic lithospheric thickness of 25–30 km (McKenzie et al., 2005). Furthermore, the toeward-dipping backstop geometry might be responsible for the localized stacking of basally accreted slivers observed in the Barbados Ridge accretionary wedge (Fig. 1D), similar to numerical experiments with such backstop geometries (Fig. 4A).

A typical structural feature of the frontal (offshore) Nankai Trough accretionary wedge along the Muroto transect is the increase in surface taper at ∼25 km and ∼40 km from the wedge toe while the oceanic crust below remains quasi horizontal (Fig. 1E). The abrupt increase in taper at ∼25 km occurs above the frontmost thrusting of the underthrust sequence, i.e., underplating (Bangs et al., 2004). The increase at ∼40 km relates to large thrust sequences rooting down to the basal décollement, forming the backstop of the underplated area (Moore et al., 2001). Sudden increase in surface slope (critical taper) may indicate that the basal décollement of the Nankai frontal wedge is significantly stronger than the intermediate one at the wedge toe (Fig. 6A). The lack of subaerial erosion of the investigated part of the wedge may explain the wide distribution of underplating relative to the width of the wedge (Fig. 7A).

The Makran is one of the largest accretionary wedge systems, with as much as 50% of the wedge consisting of underplated material (Platt et al., 1985). Underplating was interpreted to be an important factor for wedge growth due to the emergence of a low-grade metamorphic antiformal stack in the very north of the accretionary wedge (Fig. 10A; Hunziker et al., 2017). Grando and McClay (2007) showed that underplated sediments may form the footwall of large listric normal faults emerging along the coastline. Further extensive normal faulting along the Coastal Makran was reported by Normand et al. (2019). Numerical experiments of Ruh (2017) supported this idea, demonstrating that effective underplating can force the surface taper to exceed its maximum critical limit, resulting in gravitational collapse. Numerical experiments of this study are in agreement with horizontal extension above underplated areas, as illustrated by the local occurrence of tectonic underpressure (Fig. 6).

Particular features of the Makran are its vast surface exposure (Fig. 10A) and its relatively low taper (Fig. 1F; Davis et al., 1983). A 100-km-wide swath profile across the Makran shows peak elevation at the rear, a low-elevation low-relief coastal area, and a steeper offshore slope (Fig. 10B). The low onshore topography together with the low surface slope may indicate intense subaerial erosion and/or strong flexural rigidity of the downgoing lithosphere. However, experiments including hinterland-dipping backstops (as in the Makran; Fig. 1F) with intense erosion (experiment M10) and a strong basal flexure (experiment M11) cannot reproduce the significant variations between onshore and offshore surface slopes observed across the Makran (Figs. 8A, 8B, 10B; Table 1). Thus, features other than surface processes and flexural rigidity may have had a stronger effect on the morphological evolution of the wedge. In the Makran, a huge olistostrome was deposited ca. 12 Ma that emerges across large parts of the Outer and Coastal Makran (Burg et al., 2008). The implications of this olistostrome for the formation of shale-related minibasins and the overall structural evolution were severe because it lowered the surface taper drastically, enhancing out-of-sequence thrusting and underplating (Ruh et al., 2018). The implementation of a large olistostrome after 2 m.y. of shortening, filling up the trench over a distance of 100 km, leads to a surface profile resembling the natural one (Figs. 8C, 10B). I here argue that the flat surface slope of the Coastal Makran resulted from the introduction of a new, very weak décollement layer localized along the previously mentioned olistostrome (experiment M12; Fig. 8C). In fact, apatite fission-track data demonstrate that rock formations along the rear of the Makran have not been buried deeper than a maximum of 5 km, and underplated material is exposed to the surface only very locally (Figs. 8C, 10A, 10C).

The presented numerical experiments are simplified compared to natural accretionary wedges in terms of geometrical setup as well as rheological and mechanical parametrization, limiting the possibilities of direct comparison. They consist of a layered rock pile with two parallel weaker décollement layers pushed by a rigid back wall in the form of the lateral model domain boundary. With respect to other analogue and numerical studies with “sandbox”-type geometrical setups (e.g., Konstantinovskaya and Malavieille, 2011; Ruh et al., 2012), the here-introduced geometrical setup profits from an elastic flexural beam at the base of the evolving wedge, allowing for weight-related subsidence (e.g., Stockmal et al., 2007). Furthermore, limiting the modeling domain to the sedimentary wedge permits a high spatial numerical resolution (100 × 100 m) compared to mantle-scale experiments. However, several important features affecting sediment accretion along active subduction zones, and potentially underplating, are intrinsically ignored in the presented setup:

• 1. The velocity model prescribes zero outflow at the back wall of the backstop and –1 cm/yr below the basal décollement layer (Fig. 2), prohibiting any sedimentary material from exiting the model domain. This implies that all underthrust material eventually underplates and forms an antiformal stack. On the other hand, the subduction of sedimentary rocks, either undeformed or previously underplated, occurs in natural subduction zones (Bangs et al., 2003; Hilde, 1983) and large-scale numerical experiments with a dynamically evolving subduction interface (Beaumont et al., 1999; Gerya and Meilick, 2011). This inconsistency is less important for experiments with a toeward-dipping backstop, where underthrust material does not reach the lateral boundary and underplating is forced by the presence of the backstop (Fig. 4A). However, hinterland-dipping backstops are pushed upward by continued underplating, and sediments of the lower sequence would potentially exit the domain if possible (Fig. 4C).

• 2. Beaumont et al. (1999) demonstrated that subduction zone retreat and advance and affect the mode of accretion of material from the downgoing to the overriding plate. The presented geometrical and velocity conditions prohibit such a dynamic evolution of the subduction zone, which potentially influences underplating and overall mechanics of accretionary wedge systems.

• 3. The introduced sedimentary sequence hosts two mechanically weak décollement layers with static strength parameters (Fig. 2). This implementation predefines the localization of subhorizontal shear strain related to frontal and basal accretion as well as the thickness of underplated thrust slivers (Fig. 3). In contrast, natural subduction interfaces are dynamically evolving shear zones responding to local variation in lithology, fluid pressure, and stress coupling (Agard et al., 2018), which has also been demonstrated in numerical experiments (Menant et al., 2019).

• 4. In the presented experiments, purely sedimentary accretionary wedges develop, while the underlying oceanic crust exits the model domain without any internal deformation or structural incorporation into the deformed wedge (Fig. 2). In contrast, many natural examples of exhumed accretionary wedges indicate the involvement of mafic rock types (Escuder-Viruete et al., 2011; Kameda et al., 2017; Kimura and Ludden, 1995; Meneghini et al., 2009). Furthermore, seafloor roughness related to horst and graben structures in the downgoing oceanic crust can result in a patchy distribution of plate coupling, enhancing underplating (Tsuji et al., 2013).

The simplifications that were made in the experiments of this study prohibit the investigation of the above-discussed features observed in natural settings and other modeling studies. On the other hand, it allows for better control and quantitative analysis of the tested factors, such as erosion, static décollement strength, and flexural rigidity. The fact that many features related to underplating observed in nature could be reproduced by the simplified experiments underlines their validity.

Two-dimensional numerical experiments were conducted to test the influence of backstop geometry, flexural rigidity, décollement strength, and surface erosion on underplating during the growth of accretionary wedge systems. Results reveal that toeward-dipping backstops allow the formation of a major backthrust, facilitating exhumation of underplated material. In general, localization of underplating and subsequent antiformal stacking, in contrast to wide horizontal distribution of basal accretion, are enhanced by increased flexural rigidity, a larger strength difference of basal and intermediate décollements, and intensified surface erosion of the rear part. P-T paths of numerically underplated and exhumed rock slivers show typical steeper prograde and shallower retrograde evolution, as reported from natural accretionary wedges. Due to fixed geothermal gradients in experiments presented here, the P-T paths illustrate the effect of tectonic overpressure, which is positive during underthrusting and basal accretion, and negative during exhumation and at shallow levels in the rear part of wedges. Although not all potential processes important for underplating have been taken into account, some first-order comparisons to natural systems can be drawn, explaining morphological and structural observations from the Kodiak, Barbados Ridge, and Nankai Trough accretionary wedges. Furthermore, the flat onshore and steeper offshore surface taper of the Makran accretionary wedge could be explained by the deposition of a large olistostrome, further enhancing underplating by lowering the overall taper.

I thank Nathan Bangs, an anonymous reviewer, and the editors for their constructive and detailed comments that helped improve the quality of the manuscript.