Controls on Continental Strain Partitioning Above an Oblique Subduction Zone, Northern Andes

Strain partitioning onto margin‐parallel thrust and strike‐slip faults is a common process at obliquely convergent plate margins, leading to the formation and migration of crustal slivers. The degree of strain partitioning and rate of sliver migration can be linked to several factors including the angle of convergence obliquity, the dip angle of subduction, frictional coupling between the plates and the strength of the upper plate, among others. Although these factors are known to be important, their relative influence on strain partitioning is unclear, particularly at natural margins where the factors often vary along strike. Here we use a 3‐D mechanical finite‐element model to investigate the relationship between continental crustal strength, the convergence obliquity angle, the subduction angle, and strain partitioning in the Northern Volcanic Zone (NVZ) of the Andes (5°N–3°S). In the NVZ the subduction dip and obliquity angles both vary along strike, weaknesses in the continental crust may be present in suture zones or regions of arc volcanism, and strain partitioning is only observed in some regions. Thus, it is an ideal location to gain insight in which of the factors have the largest influence on deformation and sliver formation in the upper plate. Our numerical experiments confirm that a moderately high obliquity angle is needed for partitioning and that a continental crustal weakness is also required for movement of a coherent continental sliver at rates similar to geodetic observations from the NVZ. In contrast, the subduction dip angle is only of secondary importance in controlling strain partitioning behavior.


Introduction
Subduction is the fundamental tectonic process recycling lithospheric material into the Earth's interior, but also one that causes shortening, extension, and lateral shearing in the upper (overriding) plate depending on the relative sense of motion of the plates and their interface. The majority of subduction margins accommodate oblique subduction (Philippon & Corti, 2016), and this oblique motion can affect deformation of the upper plate. In regions of oblique subduction, upper-plate deformation can be understood in the context of two end-member models: oblique thrusting, where the dip-slip and strike-slip components of the plate movement are accommodated by the plate-bounding fault, and partitioned strain, where convergence occurs on the plate-bounding fault and a separate strike-slip fault (or fault system) accommodates the margin-parallel residual movement (Figure 1; e.g., Fitch, 1972;McCaffrey, 1992). The partitioned strain end-member model results in separation of upper plate material proximal to the plate margin from the rest of the upper plate, and translation of a forearc, orogenic wedge, or continental sliver along strike (e.g., Fitch, 1972). In nature, many oblique subduction margins fall between the oblique-thrusting and partitioned-strain end-members, with either poorly developed margin-parallel strike-slip fault systems (e.g., Ryukyu trench, Chemenda et al., 2000; Lesser Antilles-Puerto Rico subduction system, Symithe et al., 2015;and Central Andes, Veloza et al., 2012) or incomplete partitioning of the margin-parallel velocity component onto them (e.g., North Island of New Zealand, Barnes et al., 1998; Kuril forearc sliver, Loveless & Meade, 2011;Sumatran subduction zone, McCaffrey, 2009;and Northern Andes, Nocquet et al., 2014). The distribution of deformation onto various structures in the upper plate, however, is related not only to the plate motions but also to the geometry of the subduction zone and the mechanical properties of the plate interface and upper plate (e.g., Beck, 1983;McCaffrey, 1992). Thus, understanding the partitioning of strain onto various faults in the upper plate and the motion of continental slivers requires consideration of all these components.
The most significant factors associated with the partitioning of strain in regions of oblique convergence are a high angle of convergence obliquity, strong coupling between the two plates, and the presence of weak zones in the crust of the upper plate (e.g., Beck, 1983;Bilek, 2010;Kimura, 1986). Analytical solutions for both viscous and plastic accretionary wedges have been used to estimate the conditions under which strain is partitioned in regions of oblique convergence based on minimizing energy dissipation (e.g., Beck, 1983) or a force balance (e.g., McCaffrey, 1992;Platt, 1993). In these simple analytical solutions, an increased area of the subduction interface within the brittle domain (i.e., shallow subduction dip or subduction of old oceanic lithosphere) and a high angle of convergence obliquity were both found to increase horizontal shear stresses in the upper plate, rotating the stress tensor (e.g., Enlow & Koons, 1998). This favors development of margin-parallel strike-slip fault systems and strain partitioning. Strain partitioning is also promoted by a strong coupling between the two plates along their interface, which could reflect high friction resulting from roughness of the downgoing plate (ridges, seamounts, etc.), thin sediment thickness entering the subduction zone, subduction of buoyant oceanic crust, or fault zone composition among other factors (e.g., Bilek, 2010, and references therein). In addition, geological features thought to reduce the strength of the upper plate, such as inherited faults, suture zones, or arc volcanoes (e.g., Beck, 1983;Dewey & Lamb, 1992;Diament et al., 1992;Kimura, 1986), favor strain partitioning by lowering the shear stress needed for the occurrence of strike-slip faulting in the upper plate. While all of the factors above enhance strain partitioning, it is important to note that they also interact such that individual influences alone may not be sufficient for the development of a partitioned subduction margin (e.g., Beck, 1983;Fitch, 1972). For example, a high convergence obliquity may not necessarily result in strain partitioning. The convergence obliquity for the Alpine Fault in New Zealand is~74°, yet strain partitioning only occurs in the shallow crust (depths less than 500 m), possibly because the Alpine Fault is weak with a relatively high dip angle of~50°and the crust in the upper plate has been thermally weakened by rapid surface erosion (e.g., Koons et al., 2003;Norris & Cooper, 2001;Upton et al., 2018).
Although the analytical solutions presented above and past modeling studies provide the framework needed to understand the first-order factors influencing strain partitioning, there are several simplifications that can make comparison of their predictions to natural settings challenging. By their nature, analytical solutions are typically simplified from natural systems in several ways. For example, the analytical solutions used to predict strain partitioning behavior (e.g., Beck, 1983;McCaffrey, 1992;Platt, 1993) have geometries that vary only in two dimensions with an infinite length along strike, planar subduction interfaces, and limited spatial variations in their mechanical properties. This restricts the extent to which the predictions can be compared to natural oblique subduction systems, where the angle of obliquity or subduction dip angle may vary along the length of the system. Analog modeling studies can provide a more flexible tool for exploring factors related to strain partitioning, but experimental designs are often simplified in order to explore the effects of a small number of factors on model evolution. For example, Burbidge andBraun (1998) andMcClay et al. (2004) considered the occurrence of strain partitioning only as a function of obliquity angle with homogenous crustal materials and found that strain partitioning only occurred at obliquity angles above~45°. Chemenda et al. (2000) investigated the effects of the presence or absence of weaknesses corresponding to the volcanic arc and variations in interplate friction on strain partitioning with a constant convergence obliquity, finding strain partitioning only occurred when interplate friction was high. As noted above, obliquity, interplate friction, and the presence of absence of crustal weak zones vary in natural systems. Numerical modeling studies allow for greater flexibility in terms of geometry and material behavior but are still often restricted to two-dimensional cross-sectional geometries (e.g., Braun & Beaumont, 1995) or Figure 1. End-member models for strain partitioning in oblique subduction zones. (a) "Oblique thrusting" above an oblique subduction zone. The relative motion of the subducting slab (1) and the continental plate (2) is combined as oblique thrust motion at the plate margin (3). (b) "Crustal strain partitioning" above an oblique subduction zone. The subducting slab (1) slides underneath the continental plate (2) with movement oblique to the plate margin. On the continental side, oblique movement is divided onto a (right lateral) shear zone (4) and margin-normal subduction (5). A volcanic arc (6) or other weak zone in the continental crust may facilitate margin-parallel movement of a continental sliver (7). three-dimensional models with initial subduction geometries that do not vary or insufficient spatial resolution to resolve discrete strike-slip shear zones in the upper plate (e.g., Haynie & Jadamec, 2017;Jadamec et al., 2013;Upton et al., 2003). Although some three-dimensional numerical modeling work has explored strain partitioning with a higher average spatial resolution (2.5 km; e.g., Whipp et al., 2014), factors such as variable convergence obliquity and subduction dip angle were not the focus.
In this work we systematically investigate the influence of various subduction zone geometric and mechanical factors on strain partitioning using a lithospheric-scale 3-D numerical geodynamic model. The model is designed to represent the general geometry of subduction of an oceanic plate beneath a continent with variable convergence obliquity, subduction dip, and with or without a weak zone in the continent. This provides insight into how each of these subduction zone parameters influences the formation of a continental sliver in a 3-D obliquely convergent subduction zone. The goal is not only to understand the main factors controlling strain partitioning in a finite-length oblique subduction zone with a curved subducting plate but also to compare the results to the Northern Volcanic Zone (NVZ) of the Andes (5°N-3°S), where the Nazca Plate subducts obliquely beneath South America. This region is host to a well-defined continental sliver that has developed since 15 Ma (e.g., Aizprua et al., 2019;Alvarado et al., 2016;Nocquet et al., 2014; note that herein we use "a" for years ago and "yr" for the duration of 1 year) and provides a study region where 3-D variations in several geological factors may contribute to the development and motion of the continental sliver. In addition, this is a region of incomplete strain partitioning, where questions about what controls the magnitude of partitioning of the lateral shearing component of the convergence velocity onto various strike-slip faults can be addressed.

Geological and Tectonic Setting
The Andes provide an excellent natural example of a subduction system with significant along-strike changes in subduction parameters such as the convergence obliquity and subduction dip angle, as well as variations in the occurrence of strain partitioning. The Andean orogen extends for about 7,000 km along the entire western coast of South America. This impressive orogen formed as a consequence of subduction of the Phoenix and Farallon/Nazca oceanic plates beneath the South American plate starting in the Triassic (e.g., Cobbold et al., 2007;Müller et al., 2019). The modern mountain system known as the Andes resulted from plate-tectonic activities during the Cenozoic (Oncken, Hindle, et al., 2006). At present, both subduction dip and obliquity angle vary along the margin. The strength of the continental crust is also not homogeneous along strike but varies due to the presence of fault and shear zones such as the Chingual Cosanga Pallatunga Puná (CCPP) shear zone in the Northern Andes (Figure 2; e.g., Alvarado et al., 2016;Yepes et al., 2016) or due to the presence or absence of an active volcanic arc. The remainder of this section presents an overview of the geological history of the Northern Andes that is relevant to the present-day tectonics of this region. A general overview of the geological history of the Andes can be found in Oncken, Chong, et al. (2006) and Ramos (2009), for example.
The Northern Andes (Figure 2) occupy the western parts of Colombia and Ecuador with a geological history that is complex and different from that of the Central Andes to the south. Here we highlight essential details of the geological history of the Northern Andes and how this region differs from the Central Andes, based largely on the detailed geological history of the Northern Andes presented in Kennan and Pindell (2009). Crustal compression or extension resulting from subduction in the Central Andes was fairly continuous since at least the Jurassic, while deformation over this time in the Northern Andes reflected back-arc basin and passive margin formation prior to accretion of a series of oceanic terranes (Gansser, 1973;). Accretion of these oceanic plateau and island-arc terranes produced a number of sutures that facilitated extensive dextral shearing during and after their formation (see terrane maps and tectonic reconstructions in Kennan and Pindell, 2009). The eastern boundary and composition of these terranes are apparent in Bouguer gravity data from this region (e.g., Mooney et al., 1979), placing it roughly in the location of the Dolores Guayaquil Mega (DGM) shear zone (Figure 2; Campbell, 1974;Gutscher et al., 1999;Jaillard et al., 2009;Moberly et al., 1982). Furthermore, paleomagnetic studies confirm and support the accreted terranes being allochthonous (Cosma et al., 1998;Roperch et al., 1987).
Deformation of the Northern Andes prior to the Miocene has likely influenced the formation of the well-defined continental sliver that is being transported parallel to the plate margin at present (Figure 2; e.g., Nocquet et al., 2014). As mentioned above and described in detail by Kennan and Pindell (2009),

10.1029/2019TC005886
Tectonics terrane accretion and associated faulting during the Cretaceous-Eocene resulted in widespread dextral shearing in the Northern Andes, displacing some terranes hundreds of kilometers parallel to the present-day Nazca-South American plate margin. Shearing during this period occurred along a fairly linear set of strike-slip faults, extending from present-day Colombia south to northwestern Peru (see Figure 5 of Kennan and Pindell, 2009). The connection of strike-slip faults in this region to the Gulf of Guayaquil occurred only in the Eocene and fault slip related to the present-day continental sliver in the Northern Andes appears to have only initiated since 15 Ma (Figure 2; e.g., Alvarado et al., 2016;). Furthermore, the depositional record in basins such as the Tumbes basin in the Gulf of Guayaquil, which opened as a result of northeastward migration of the continental sliver, suggests rapid subsidence associated with sliver movement starting only in the past~1.8 Myr (Witt et al., 2006;Witt & Bourgois, 2010). This is similar to the timing of the Pliocene arrival of the Carnegie Ridge at the subduction trench, suggesting a possible role of increased interplate coupling in driving sliver movement (Figure 2; Egbue & Kellogg, 2010;Witt et al., 2006).
Although the pre-Neogene deformation history has likely affected sliver formation and displacement, the emerging picture is one of a continental sliver involving motion on several fault segments, both new and  Gutscher et al., 1999, andYepes et al., 2016). The Nazca Plate (1) and the South American Plate (2) converge with a velocity relative to South America of 5-7 cm/year and with an obliquity angle (3) of 31-45°relative to the trench normal. A shear zone proximal to the volcanic arc is evidence of strain partitioning (4). Coast-parallel sliver movement (about 1 cm/year) of the North Andean Sliver (green and blue shaded region) is bounded by the Chingual Cosanga Pallatunga Puná (CCPP) shear zone , a region analogous to the earlier North Andes Block (green shaded area) bounded by the Dolores Guayaquil Mega Shear (DGM; Gutscher et al., 1999). The 100-and 200-km isolines of the subducting slab as well as the Carnegie, the Cocos, and minor Coba and Malpelo Ridges are indicated. Ages of the oceanic crust are shaded in gray (after Cobbold & Rossello, 2003). The black triangles are the approximate locations of Holocene volcanoes. The position of nine MW ≥7.0 earthquakes of the twentieth century in the NVZ is marked with white (MW ≥7.0) and black (MW ≥8.0) stars (year; MW): (1906; 8.8), (1942; 7.8), (1953; 7.6), (1958; 7.6), (1970; 7.2), (1979; 8.1), (1987; 7.1), (1995; 7.0), and (1998; 7.1).

10.1029/2019TC005886
Tectonics SCHÜTT AND WHIPP old. Earlier descriptions of the sliver had it bounded by the DGM (Campbell, 1974;Gutscher et al., 1999;Moberly et al., 1982), which reaches from the Gulf of Guayaquil northeastward through the Andean Cordillera of Ecuador and Colombia ( Figure 2). This right-lateral fault zone marked the eastern boundary of the sliver, referred to by those authors as the North Andean Block (NAB; green shaded region in Figure 2; Bourgois, 2013;Gutscher et al., 1999Gutscher et al., , 2000, a continental sliver being displaced at~1 cm/year northeastward parallel to the continental margin. However, the definition of the NAB has recently been revised with the name North Andean Sliver (NAS; green and blue shaded regions in Figure 2; Alvarado et al., 2016;Yepes et al., 2016), reflecting both internal deformation of the sliver that is inconsistent with the nondeforming behavior expected in the "block" model and a different tectonic boundary to the east. The tectonic map of Ecuador by Alvarado et al. (2016) now links sliver motion to the CCPP shear zone, a right-lateral strike-slip fault sharply bounding the eastern edge of the NAS, and suggests that the modern DGM is inactive at least in central-northern Ecuador and southern Colombia. Regardless, both right-lateral strike-slip faults belong to an extensive system of (abandoned, reactivated, and new) strike-slip and thrust faults located in the Andean Cordillera of Ecuador and Colombia Veloza et al., 2012).
Another factor that may have been important in the history of the NAS is the NVZ (Figure 2), the region of active arc volcanism in the Northern Andes that lies within the region of right-lateral faulting in Ecuador and Columbia mentioned above. In this area, the oceanic Nazca Plate subducts with a moderate subduction dip angle of about 25-35°and with a convergence obliquity angle of 31-45°relative to the trench normal. About 40 volcanoes (active from~100 ka to today) are located in the NVZ (Stern, 2004, and references therein). In addition to preexisting features such as suture zones, thermal weakening of the crust as a result of arc volcanism is often thought to influence the development of strike-slip faults and partitioning of strain at obliquely convergent plate margins (e.g., Beck, 1983;Fitch, 1972). In the case of the NVZ, it is possible that the volcanic arc may have influenced the locations of strike-slip faults accommodating sliver movement, though the majority of NVZ volcanoes in Ecuador lie west of the sliver-bounding CCPP fault system and are distributed over a broad region . This suggests that the volcanic arc may play a more limited role in sliver formation than in other areas with arc volcanoes more proximal to the sliver-bounding faults, such as Central America (e.g., LaFemina et al., 2009).
In addition to the fault systems mentioned above, some other tectonic and bathymetric features have also played important roles in the evolution of the Northern Andes. The most notable is the Carnegie Ridge, an aseismic oceanic plateau that marks a transition in the geometry of the subducting Nazca slab as well as having elevated topography that may increase friction in the subduction zone. The dip angle of subduction is moderate north of the Carnegie Ridge (Gutscher et al., 1999), while the ridge itself is undergoing flat subduction due to its buoyancy. In addition the ridge's higher elevation relative to the surrounding sea floor results in increased coastal relief where it pushes into the continental plate. Because of the timing of its arrival at the trench and the subsequent subsidence of basins trailing the displaced NAS, it has been suggested that increased interplate coupling above the ridge may have helped in the increased movement of the NAS since~1.8 Ma (e.g., Egbue & Kellogg, 2010;Witt et al., 2006).

Anticipated Model Behavior
Strain partitioning at obliquely convergent plate boundaries is thought to be driven by margin-parallel shear forces acting along the plate-bounding fault system and resisted by the mechanical strength of the upper plate material (e.g., McCaffrey, 1992;Platt, 1993). An additional important factor is the interplate coupling, which can also be related to the steepness of the slab subduction (e.g., Jarrard, 1986;Wortel et al., 1991). For an infinite length orogenic wedge or continental sliver, McCaffrey (1992) shows that the driving shear stress for margin-parallel motion increases in proportion with the obliquity of the convergence angle with respect to the normal to the margin and the strength of the plate-bounding shear zone (interplate coupling), and inversely in proportion to the dip angle of the basal detachment or subducting plate. In other words, the driving forces increase when the obliquity is large, the plates are strongly coupled, or the subduction dip angle is low. In addition, upper plates containing zones of crustal weakness will generally be more favorable to strain partitioning, whereas oblique thrusting is more likely when the crust in the upper plate is strong (e.g., McCaffrey, 1992;Platt, 1993). Finally, in continental slivers or orogenic wedges of finite along-strike length, the along-strike length of the sliver/wedge and the strength of the crust in the upper plate will also affect the occurrence of strain partitioning (e.g., Chemenda et al., 2000;Whipp et al., 2014). In this case the geometry of the sliver/wedge requires shortening and extensional structures at its lateral ends. These structures develop in the crust of the upper plate and thus resist strain partitioning, particularly when the crust is mechanically strong (Whipp et al., 2014). The influence of these "end" structures in determining strain partitioning behavior is also affected by the along-strike length of the sliver/wedge. They will play a larger role in resisting strain partitioning when the sliver/wedge has a length comparable to its across-strike width, and their contribution to resisting strain partitioning will decrease as the lateral length of the sliver/wedge increases.
In the Northern Andes in the NVZ, oblique subduction occurs with an obliquity angle of~40°, at a moderate subduction angle of~30°, and an active volcanic arc is present, which may act as a weakness in the crust of the upper plate. In the absence of a weak region in the upper plate, related to the volcanic arc, suture zones, or other features, we expect that any continental sliver that forms would be narrow and moving at a low velocity along strike. This is based on the theory presented above and past results from strain partitioning experiments with homogeneous upper plates (e.g., Braun & Beaumont, 1995;Burbidge & Braun, 1998;McClay et al., 2004). If the volcanic arc or another similar zone of weakness is present in the upper plate, strain partitioning is expected to be well developed with a coherent sliver moving at a relatively high velocity along strike. Changes to the obliquity angle or subduction dip angle will also affect the strain partitioning behavior, but the relative roles of these factors in partitioning strain in a subduction system with variable plate interface dip are less clear. We explore each of these parameters in series of numerical experiments described below.

Model Description
To investigate the conditions under which strain is partitioned in an oblique subduction zone and its implications for sliver movement in the Northern Andes, we focus on the influence of three factors: the strength of the continental crust (uniform or with a preexisting weak zone), the plate convergence obliquity angle, and the subduction dip angle. The models run for 1 Myr to find a stable velocity solution for the given geometry and material properties and are then compared to geological observations. Details of the model design are below.

Numerical Model
The influences of the plate convergence obliquity angle, subduction dip angle, and strength of the continental crust on sliver formation in the overriding plate are investigated using 3-D lithospheric-scale numerical experiments using the program DOUAR Thieulot et al., 2008). DOUAR solves the Stokes equations to determine the velocity fields in model materials, and the resulting uplift and strain rates. Details about the relevant equations solved in DOUAR can be found in Appendix 1. A comprehensive description of the program is in Braun et al. (2008).

Model Design
The reference model is designed to include the essential features of an oblique and moderately steep subduction zone based on the northern part of South America, the NVZ ( Figure 2). The model measures 1,600 km × 1,600 km in lateral extent with a 162.5-km thickness ( Figure 3a) with a horizontal spatial resolution of 12.5 km and vertical resolution of 3.125 km (128 × 128 × 52 elements). The topography is defined by a surface within the finite-element mesh to resemble the first-order topographic features in the Northern Andes. Sea level is at 155-km height in the model domain, from which the Andes rise at a gradient of 6% to form a 96-km-wide, 5-km-high elevated region parallel to the plate margin at a distance of about 200 km inland (Figure 3a). This upper surface does not move during the experiments, equivalent to assuming that any uplifted rock is instantaneously eroded and sediment fills any region that subsides. The subducting plate geometry is parabolic along dip, with a dip angle that varies from 0°at the plate margin to~40°at the base of the model. Both the topography and the geometry of the subducting plate are offset along strike to be parallel to the plate margin (Figures 3a and 3b) producing a central zone of oblique subduction along the plate margin with an obliquity angle γ that varies between 10 and 50° (Figure 3b). A zone of mechanical weakness in the continental crust is also included in some experiments (e.g., Model 2). It is located parallel to the plate margin and at the model equivalent location of the volcanic arc or existing suture zones in the Northern Andes (Figure 3a).

Boundary Conditions
The velocity boundary conditions that drive plate convergence are designed to simulate advancing subduction in which the subduction zone advances into the upper continental plate with time. The oceanic plate subducts beneath the continent at the subduction velocity v s = 5 cm/year. In addition, the subduction zone advances into the model continent at a rate of v adv = 2 cm/year to account for subduction zone advance into South America as shown in Figure 3a in a reference frame with the continent fixed. Combined, this produces a total plate convergence velocity of v conv = 7 cm/year (Figure 3a), similar to the observed Nazca-South America plate convergence velocity (e.g., DeMets et al., 2010;Kendrick et al., 2003). In the model, the reference frame is shifted by subtracting the subduction advance velocity v adv from the defined boundary condition velocities, which ensures that the margin of the two plates remains in the center of the model domain during the experiments. The subducting plate influx velocity along the left-hand side of the model is constant over the entire area and balanced by the out flux along the bottom of the model. Thus, the velocity field on the model base is defined for the subduction geometry and kinematically drives subduction. Assuming crustal shortening only in the upper continental plate, its mantle lithosphere is also removed at the base of the model within the subduction zone. Without removal of the continental mantle lithosphere, the mantle shortens below the continental crust and leading to excess uplift of the surface. The upper boundary of the model is open.

Material Parameters and Properties
For simplicity all model materials are either viscous or frictional plastic, and the model material properties do not vary as a function of temperature (i.e., purely mechanical models). Relevant material properties are Generalized topography of the study area including a continental crustal weakness (e.g., volcanic arc or suture zone) indicated with the gray box, which corresponds to the gray zone in the side view. Basal surface: Oblique (angle γ) zone of subduction between Nazca Plate and South American Plate. All gray zones have the same material properties as the host material except a friction angle φ = 2°instead of φ = 4°(slab) and φ =15°(continent), respectively. The shades of green represent the vertical velocity field of the oceanic plate with light green for zero vertical velocity, and the darker shades indicate increasing downward velocity in the subduction zone following the shape of the slab. White along the base of the continental side is for zero vertical velocity. largely adopted from Currie et al. (2015) and Beaumont et al. (2006;Figure 3c). Other than the asthenosphere beneath the oceanic plate (Figures 3a and 3c), a high nominal viscosity ensures frictional plastic deformation in all other model materials (see Appendix 1 for details on how plasticity is implemented in DOUAR). These materials are thus plastic without diffusion or dislocation creep. Variations in the internal angle of friction of these materials thus allow exploration of the effects of mechanical strength on deformation in the model.
Two main sets of experiments are considered, one with and one without a weak zone in the continental crust (Figures 3a and 3c). The friction angle for the normal (unweakened) continental crust is φ = 15°, reduced from the typical friction angle of~30°for dry sand by hydrostatic pore fluid pressure. This value is further reduced to φ = 2°within the crustal weak zone in models containing an upper-plate crustal weakness. In addition, the interface between the two plates is weak in all cases in order to facilitate subduction of the oceanic plate. Because of the resolution of our experiments, the subduction interface strength is specified as that of the oceanic crust. Thus, the model oceanic crust has a low initial friction angle φ = 4°, which is further reduced to φ = 2°along the plate interface at the start of the experiments to simulate the effects of strain softening, the reduction of rock strength as it undergoes deformation (e.g., Read & Hegemier, 1984), along the plate interface. During the experiments, the strength of the oceanic crust is continually reduced along the plate interface from φ = 4°for total strain values at or below 50% to a minimum value of φ = 2°for strain above 150%. Other than this small strength reduction, the plate interface (oceanic crust) strength does not vary in our experiments. Finally, we note that the absolute values of the friction angles in our experiments may differ from natural systems where parts of the crust deform in a viscous manner (crystal plasticity). However, we anticipate the same general behavior observed in our experiments when the subduction interface and regions of the continental crust are weakened relative to the strength of typical continental crust.

Model Assumptions and Restrictions
The uppermost surface in the models is fixed for simplicity. Several experiments with a free surface were also run, but unreasonable uplift or subsidence occurred in some cases leading to numerical problems. For example, subsidence often occurred near the plate margin as a result of the subduction fault system being distributed across several elements in the model. Changes in the geometry of the free surface would affect the isostasy and gravity of the whole system. However, as these models are not thermomechanical, where variations in temperature might affect the density of model materials, the isostatic effects of a changing free surface are not included. In addition, the anomalous uplift and/or subsidence in the early model designs with a free surface that led to the numerical problems did not affect the first-order tectonic behavior of the models and the development of strain partitioning (see Schütt, 2018, for additional discussion). Thus, we have chosen to use a fixed upper model surface.
The slab and its subduction path are kinematically defined using the velocity boundary conditions instead of driving subduction by slab buoyancy. This may influence the evolution of the subduction zone geometry over long time scales. However, because the model is only 160 km thick, there is insufficient space to include the deeper parts of the subducting plate that would gravitationally induce subduction. In addition, a deeper (i.e., upper mantle scale) model design would reduce the spatial resolution of the continental crust, which is the focus region for understanding strain partitioning behavior in this work. Lastly, kinematic subduction of the slab in all runs allows for more direct comparison between the different parameters for slab subduction angle and obliquity.
The weak zone in the continental crust, which could represent the volcanic arc or an existing structural weakness such as a suture zone, is a continuous and linear weakness, and the continental crust does not strain soften. The choice of continental weakness geometry is an end-member but is a simple geometry that allows the frictional strength of the crust to be varied along the length of the oblique segment. We also do not consider further weakening of the continental crust in the upper plate via strain softening. The strength of natural rock is often reduced when deformed (e.g., Read & Hegemier, 1984), which can act to localize deformation and introduce strength heterogeneity to otherwise homogenous rock. We do not include this effect because our focus is on deformation of the upper plate of an obliquely convergent subduction margin either completely lacking a weak zone in the crust or with a preexisting weak zone (volcanic arc, suture zone, etc.) present from the start of our experiments. While strain softening may be important in the development of strain partitioning in these settings, we have chosen to focus on the cases with and without a well-defined weak zone to clearly identify the effect of the presence or absence of such a feature.
Finally, temperature and its effects on the rheology of the model materials are not included (i.e., the models are purely mechanical). This means no weakening of the crust, lithosphere, or subducting slab due to diffusion/dislocation creep or (partial) melting is included. Experiments including temperature-dependent material properties are beyond the scope of this work, but we have conducted some experiments including temperature dependence for the NVZ reference model (Schütt, 2018) and preliminary results do not differ significantly from the purely mechanical results of sliver formation presented here. Because of this choice, these models are most appropriate for comparison with to relatively cold subduction margins where deformation is dominated by faulting and frictional sliding. Comparisons should be done more carefully in locations where the upper plate has been thermally weakened due to rapid surface erosion (e.g., New Zealand).

Results
A suite of numerical experiments (Table 1) have been conducted in order to investigate the influence of the strength of the continental crust (Models 1 and 2), the subduction dip angle (Experiment Series 1), and the convergence obliquity angle (Experiment Series 2) on strain partitioning in obliquely convergent subduction zones.

Model 1: Reference Model With Homogeneous, Strong Continental Crust
The reference model includes a homogeneously strong continental crust, an obliquity angle of 40°, and an average dip angle of about 24° (Figure 3 and Appendix 2). After 1.0 Myr of model time, subduction of the oceanic plate has led to shortening in the upper plate and formation of a narrow zone of material being transported along strike ( Figure 4). The average northward velocity of the narrow sliver in this model is 0.4 cm/ year (pink region along upper surface of model in Figure 4). A slightly higher velocity (approximately 0.8 cm/year; cyan region along upper surface of model in Figure 4) appears in close proximity to the margin,  in the model equivalent of the accretionary prism. The sliver that forms here is only present close to the plate margin and develops only above the slab where the subduction angle is shallow (<10°). The total velocity vector field (black arrows in Figure 4) shows a strong downward movement following the curvature of the subducting slab. This movement stems from the basal velocity boundary conditions. Except for the minor shortening in the overriding plate close to the margin, the velocity field is almost zero in the continental (upper) plate.

Experiment Series 1: Varying the Subduction Dip Angle
In this set of experiments the subduction angle was varied to investigate its effect on the behavior of the continental crust proximal to the margin. The dip angle was varied by changing the location at which the slab exited the model relative to the plate margin as well as the location of the weak zone between the slab and the overriding plate. In this experiment series we considered average dip angles from 21 to 42°. Details on how the average dip was calculated can be found in Appendix 2.
As expected, decreasing the average subduction dip angle leads to strain partitioning across a wider part of the upper plate in a homogeneously strong continental crust ( Figure 5). In the case of pure oblique thrusting along the subduction interrace, the margin-parallel velocity would be zero in the upper plate (i.e., no partitioning). Thus, the increase in margin-parallel velocity for different average subduction dip angles in Figure 5 reflects an increase in the amount of strain partitioning. The velocity of the sliver peaks at the plate margin and decreases exponentially inland. In addition, a more rapid exponential decrease is observed for higher subduction slab dip angles. For smaller dip angles, the velocities are stable with around 0.5 cm/year over a distance of around 200 km roughly 100 km inland from the margin. Although the margin-parallel velocity component increases for experiments with smaller dip angles, strain partitioning is still only very weak and without a clearly defined boundary for the partitioned region.

Model 2: Including a Crustal Weakness in the Upper Plate
For this experiment a weakness is included in the crust in the overriding plate ( Figure 3). The addition of the linear zone of weakness parallel to the margin results in a strong northward moving sliver as a result of a decreased continental strength resisting strain partitioning. It can clearly be seen that the sliver stretches from the plate margin to the continental crustal weakness and that it occupies the entire area from the surface to the slab interface (blue-green region in Figure 6). The sliver is transported northward parallel to the margin with a velocity of up to 1.3 cm/year, corresponding to roughly 40% of the margin-parallel velocity for complete partitioning.

Experiment Series 2: Varying the Convergence Obliquity Angle
This set of experiments explores the effect of the convergence obliquity angle on strain partitioning. The rate at which a continental sliver can migrate along strike will ultimately be limited by the plate convergence velocity and obliquity angle. However, a higher obliquity angle alone does not ensure more rapid margin-parallel sliver movement because of the influences of the strength of the continent and friction along the subduction interface, for example. In the following, experiments with obliquity angles varying from 10 to 50°are conducted in two sets: one with a homogeneously strong continental crust and a second with the continental crustal weakness included.
As would be expected for increasing components of convergence parallel to the margin, velocities near the margin and in the continental sliver increase with increasing obliquity angle (Figure 7). The northward velocity component maxima in all cases of these experiments are at or near the plate margin. In the cases where no continental weakness is present, the northward velocity drops almost immediately across the plate margin. In contrast, the velocity in cases with a continental weakness present remains high over several hundred kilometers and across the entire area between the plate margin and the continental weakness (gray bar in Figure 7) for obliquity angles larger than 20°. It is thus clear that sliver formation is enhanced in the cases that include a continental weakness. In these experiments the sliver formation begins at an obliquity angle of 20°. The velocity is rather small (between 0.3 and 0.6 cm/year), but for an obliquity angle of 30°, the entire area between the plate margin and the continental weakness is part of the sliver. In contrast, the model without the continental weakness does not show any sliver movement at 30°obliquity angle. In the experiments without the continental weakness sliver movement is only observed for an obliquity angle of 50°, similar to past numerical and analog modeling studies with a homogenous continental crustal strength (e.g., Braun & Beaumont, 1995;Burbidge & Braun, 1998;McClay et al., 2004). In cases with smaller angles, the northward  directed velocity component is only high in the trench region. This is due to very shallow subduction angle proximal to the margin as was observed in section 5.2 (variable subduction dip angle experiments).

Comparison of Model Results to Observations in Nature
Comparing the obtained results to geodetic observations in nature, it stands out that the sliver observed in Model 2, which includes a continental crustal weakness, is in very good agreement with the extent and velocity of the sliver observed in the NVZ (Figure 8). The weak zone seems not only to provide a predetermined breaking point for the movement of a coherent sliver but also to reduce the shear stress needed for strain to partition. The obliquity angle in Model 2 is sufficiently large, the subduction dip angle is below the high values that might reduce the driving shear stress, and the sliver is translated along strike at around 1.3 cm/year, similar to the appearance of the NVZ sliver described by Nocquet et al. (2014), for example. The sliver translation velocities observed in nature (1.4 cm/year; Nocquet et al., 2014) are nearly equal to those in the model, and sliver movement is well defined when the continental crust has a weakness in the vicinity of the model volcanic arc or suture zone. We do note, however, that the models utilize friction angles that are not based on observations from the study region. This topic is discussed further in section 6.3.

Which Subduction Zone Parameter Has the Largest Influence on Strain Partitioning?
The reference model (Model 1) with a homogeneously strong crust does not show significant strain partitioning or sliver movement in contrast to Model 2 with a weak zone present in the continental crust (compare Figures 4 and 6). The two series of experiments varying the subduction dip angle and convergence obliquity, respectively, show that strain partitioning is favored when the subduction dip angle is shallow or the obliquity angle is high. This raises the question: which of the geometric or mechanical factors considered in this work has the largest role in affecting strain partitioning in obliquely convergent subduction zones?
As we have shown, the lack of a weak zone in the crust of the upper plate does not produce any significant sliver movement (i.e., Figure 4), other than in the regions of very shallow subduction dip. In Model 1, a slightly higher margin-parallel velocity (approximately 0.8 cm/year) appears in close vicinity of the margin (e.g., Figure 4). This zone of more rapid along-strike movement is restricted to the immediate margin, where the overriding crust is very thin and more easily detached from the rest of the continent. This area can be compared to the area of an accretionary prism. Although this would qualify as weak partitioning in the proximal margin, this behavior is clearly not formation of a continental sliver. In contrast, Model 2 includes a continental crustal weakness and shows a coherent, rapidly moving sliver ( Figure 6). This means that the presence of a volcanic arc or other similar zone of weakness is cutting the continental plate such that the wedge above the subducting slab between the margin and the weak zone is treated as an independent block under the strong influence of the obliquity of the subduction. Further, this means that such a continental weakness is necessary for a proper sliver to develop and that obliquity alone is necessary, but not sufficient for strain partitioning in a moderately steep subduction zone.
Experiment Series 1 shows how strain partitioning is dependent on the subduction dip angle ( Figure 5). The areal extent of the strain partitioned region is increasing with decreasing subduction dip. With a moderate dip angle the interface between the sinking slab and the overriding plate in the frictional domain of the upper lithosphere is small, and thus, the area of friction reduced. As the plate subducts into the viscously deforming regime of the lithosphere closer to the margin, the frictional driving shear stresses at the interface are significantly reduced and no sliver movement likely. This is consistent with the highest velocities being observed in direct vicinity of the plate margin.
In Experiment Series 2, focusing on the effect changes in obliquity angle have on sliver formation and movement, it is clear that the obliquity angle is a first-order factor for strain partitioning within the range of parameters we consider, though average rates of sliver movement reflect partitioning of only 15-45% of the subduction velocity. The results of this experiment series show that the sliver width and rate of translation decrease with decreasing obliquity angle (Figure 7). We can also see that the sliver in the models is moving at a rate below that expected for complete partitioning of the subduction velocity (i.e., less than the full margin-parallel velocity component). For example, the expectation for complete partitioning of a 5 cm/year subduction velocity for a 40°convergence obliquity would have the sliver moving parallel to the margin at around 3.2 cm/year. In contrast, we find a maximum margin-parallel rate of sliver movement around 1.3 cm/year, and an average rate of northward motion of roughly 0.9 cm/year, corresponding to 37-40% partitioning. A velocity below that for full partitioning is an expected result based on predictions from analytical solutions for strain partitioning (e.g., McCaffrey, 1992;Platt, 1993) and reflects the resistance of the continental weak zone to sliver motion. In other words, the sliver motion will only begin when the shear stress from oblique subduction reaches the required yield stress in the upper plate, and motion will increase with increased obliquity from that point. This is observed in Experiment Series 2 (Figure 7).
Thus, based on the range of parameters explored in our experiments we can conclude that the first-order factors for strain partitioning are the obliquity angle and the presence of a weakness in the continental crust, while the subduction dip angle has only a secondary influence. In spite of the lower sensitivity to subduction dip angle, we note that a well-defined sliver can be obtained in experiments that combine a high obliquity angle or continental crustal weakness with a steeper subduction dip angle (see Schütt, 2018, for additional experiments). In addition, despite the simplicity of these experiments (purely mechanical), results that compare well to geodetic observations have been obtained. Furthermore, results with thermo-mechanical models show similar behavior for the NVZ region (Schütt, 2018).

How Strong Is the Continental Crust and How Strong Is It Proximal to a Volcanic Arc?
The continental weak zone corresponding to a volcanic arc or preexisting shear/suture zone in the continental crust has been implemented in our experiments by reducing the friction angle in this region, but is the reduction used in this study reasonable in nature? When the crust is thermally weakened by a volcanic arc, the friction angle does not likely change significantly, but rather the temperature increases and reduces the thickness of the crust where frictional sliding occurs (similar to the effects of erosion; e.g., Koons et al., 2002). To estimate the magnitude of this effect, we calculated the integrated strength of the a 35-km-thick continental crust in numerical experiments that include temperature-sensitive rheologies using the implementation of visco-plasticity in DOUAR  and the wet quartzite properties from Currie et al. (2015). We consider three different scenarios (Figure 9 and see Appendix 3 for details): (1) a typical continental crust (friction angle 15°; blue dotted line in Figure 9), (2) a weak (low friction angle) crust (friction angle 2°; red dotted line in Figure 9), and (3) a thermally weakened continental crust (friction angle 15°and temperature fixed at 900°C in a "magma chamber"; green line in Figure 9). The magma chamber in this case is meant to represent the thermal weakening effects over some thickness of the crust, as a first approximation. Additional detail about the material properties and thermal model design can be found in Schütt (2018).
In this rough estimation one sees that with typical continental crust as a reference, reducing the friction angle to 2°results in approximately the same integrated strength as introducing a magma chamber affecting a 15-km-thick section of the upper crust in the weak zone. Further consideration of the thermal and mechanical effects of arc magmatism in the crust is beyond the scope of this paper.
At the scale of a crustal shear zone, estimation of the effective friction angle for rocks in the shear zone is difficult (e.g., Burov, 2011). Thus, our approach for selection of suitable frictional angles is based on a combination of theoretical predictions and testing via numerical experiments. The theory for predicting conditions for strain partitioning developed in Whipp et al. (2014) allows for estimation of the crustal-scale friction angles in regions of oblique convergence with observed strain partitioning as a function the strength of the subduction interface and the convergence obliquity angle. Preliminary experiments using our model design indicated that the subduction interface must be weak even for kinematic (i.e., not buoyancy-driven) subduction to occur. This is consistent suggestions of subduction interface strength (e.g., Duarte et al., 2015) and limits the possible range of friction angles for the continental crustal and weak zone. Considering the weak subduction interface and theory of Whipp et al. (2014), it was estimated that the rear shear zone must be quite weak and similar to that of the subduction interface for strain partitioning to occur. Although the friction angle used here for the weak zone (2°) is quite low, preliminary experiments indicated that higher friction angles prevented strain partitioning and decreasing it below 2°led to shortening along the rear shear zone, consistent with observations from Whipp et al. (2014). We note that these numerical experiments were run at an intermediate resolution, so it is possible that the friction angles for the subduction interface and continental weakness are slightly lower than their natural equivalents. This is because shear zones in the numerical model are typically distributed over three to four elements, which could result in shearing both within the weak zone and into the nearby (stronger) continental crust.

Where Does the Sliver Go?
Once the sliver has formed and the continental crust is moving parallel to the margin, the sliver will be pushed into crust positioned further up the margin, requiring either local structures to accommodate shortening or space into which the sliver can be displaced. In addition, the sliver material either needs to be replaced or displacement will cause local extension near one lateral end of the sliver. In other words, the geodynamics of sliver movement include not only motion of the sliver itself but also the more complex tectonics around the sliver.
As seen in structural maps of the region north of the NAB (see Figure 2; e.g., Symithe et al., 2015), there are several structures within the South American continent and at its northern boundary that accommodate northeastward translation of the NAS. The two most prominent structures in the area of the boundaries of the Caribbean and South American Plates are a shear zone running SW-NE with a turn to W-E through Venezuela and a subduction zone further north, where the Caribbean Plate is subducting underneath the South American Plate. Global Positioning System measurements show that the crust west of the shear zone is moving north-eastward as an extension of the sliver movement observed further south (Symithe et al., 2015). This is overriding the strongly eastward moving Caribbean plate and accommodating shortening owing to sliver movement.
At the opposite end of the NAS near the Gulf of Guayaquil, sliver movement would be expected to cause local extension within the South American continent and its proximal margin. Strike-slip motion on the DGM and/or CCPP faults clearly accommodates escaping along faults of northern Colombia and Venezuela toward the Caribbean Plate (Gutscher et al., 1999;Symithe et al., 2015), but which faults may be involved in sliver motion near the Gulf of Guayaquil is less clear. However, as a result of the sliver movement, local extension has occurred in the region of the Gulf of Guayaquil, as expected for the trailing edge of a continental sliver. Offshore basin records suggest an onset of local subsidence of the proximal margin in this region at least since the earliest Quaternary (e.g., Witt et al., 2006;Witt & Bourgois, 2010), which appears to be related to the partitioning of strain along strike and motion of the NAS. Although our focus in this work is on the conditions under which strain partitioning occurs, we observe similar extensional and shortening structures at the leading and trailing edges in the numerical experiments where a continental sliver develops (e.g., Model 2).

Application to Global Subduction Zones
The experimental results presented above are also applicable to other obliquely convergent subduction zones globally. The Sumatra-Java subduction zone is characterized by a typical subduction slab dip angle of around 30°and an approximate obliquity angle of 30°(e.g., Hayes et al., 2012). Based on the presented results, the overriding plate would be expected to develop sliver movement at a slightly smaller rate compared to the NVZ due to the lower obliquity angle. In accordance with this prediction, sliver movement of~0.6 cm/year is observed along the Sumatra Fault, which runs approximately parallel to the subduction margin (e.g., McCaffrey, 1992).
The Kuril-Kamchatka Japan subduction zone is geometrically more complicated due to the triple point of the Eurasian, Philippine Sea, and Pacific Plates but can still be related to the results of our experiments. From northern to southern Japan, the subduction zone characteristics change from shallow (20°) to moderate (40°) subduction dip angles and from a very small obliquity angle (0-10°) to higher obliquity angle (30-40°;Fitch, 1972;Hayes et al., 2012). Our experiments would predict little to no strain partitioning in the northern segment of the Japan subduction zone (north of 35°N) and intermediate sliver movement in the southern segment (south of 35°N). These predictions are all excluding the possible local variations in lithospheric strength and subduction geometry. Again, similar to the predictions from the numerical experiments, the sliver movement observed in the northern segment of the Kuril-Kamchatka Japan subduction zone is near zero, while the southern segment has developed a sliver moving at approximately 0.6 cm/year along strike, the Nankai Forearc sliver (Fitch, 1972).
A third example of sliver movement owing to strain partitioning is the Cascadia subduction zone along the northwestern margin of North America. The subduction slab dip angle in Cascadia is low, and the setting can be classified as flat slab subduction (Gutscher et al., 2000;Hayes et al., 2012). Our numerical experiments suggest no sliver movement for obliquity angles less than 30°and only moderate sliver translation velocities for obliquity angles larger than 30°. Cascadia is characterized by an obliquity angle of around 20-30°and a convergence velocity of~3.4 cm/year, yet sliver movement of~0.7 cm/year is observed (Flesch, 2007). The Cascadia subduction zone is one example where the generalization of the model results from this work may not hold true. Here the large-scale influence of the adjacent plate tectonics also needs to be considered. Sliver movement in Cascadia can be an escape window for the margin-parallel running San Andreas Fault, which is to the south adjacent to the Cascadia region (Flesch, 2007). Alternatively, the plate coupling between the North America and Juan de Fuca Plates may be stronger than that used in our experiments based on the Northern Andes, though there appears to be little consensus on the strength of coupling of the subducting and overriding plates in the Cascadia region (e.g., Chapman & Melbourne, 2009;Wang et al., 1995). Further comparison of our results to Cascadia would require additional experiments varying the interplate frictional strength, which is beyond the scope of this work.

Conclusions
The general results of our 3-D numerical experiments of oblique plate convergence and subduction suggest the following: 1. The plate convergence obliquity angle has a larger influence on the strain partitioning in the overriding plate than the dip angle of the subducting plate or the presence of a weak zone in the crust of the upper plate. 2. A continental weakness in the form of a volcanic arc or preexisting structure strongly enhances the formation of a coherent continental sliver. 3. The magnitude of partitioning of the margin-parallel velocity component onto strike-slip shear zones in the upper plate is incomplete, with only up to 50% of that velocity component observed in sliver motion. 4. The subduction dip angle and the width of a continental sliver that forms as a result of strain partitioning are inversely proportional. In other words, shallow subduction favors formation of a wider continental sliver. This phenomenon is observed in all of the experiments, where a small part of the continental crust closest to the margin will always partition due to the very shallow subduction dip and limited shear resistance force in the thin continental crust.
In the context of the Northern Andes, the model with a continental crustal weakness represents the conditions observed in the NVZ well: sliver formation and its rate of movement are in good agreement with observations in nature despite of the simplicity of our (purely mechanical) models.
where i, j, and k are symbolic coordinate indices and σ′ is the deviator, defined as where Tr[σ] describes the trace of σ and I is the unit matrix. The Mohr-Coulomb criterion used in the code for plastic deformation can now be expressed using stress invariants where ϕ is the material angle of internal friction and θ L is the Lode angle (A7) DOUAR calculates stresses from the strain rate tensor _ ε and its invariants, instead of directly from the stress tensor. Accordingly, the first, second, and third strain rate tensor invariants are I 1 ¼_ ε ii ; (A8a) where the deviatoric strain rate tensor _ ε ′ is Tr I 1 ½ I: In each model iteration, the elemental effective stress The elemental effective viscosity η eff is recalculated in order to place the element on yield, whenever the elemental effective stress exceeds the Mohr-Coulomb yield strength σ y : Because the effective viscosity depends on the strain rate and the strain rates depend on the viscosity (i.e., the system is nonlinear), velocities are calculated iteratively. The velocity iterations are considered converged when the L 2 -norm of the difference in velocity for each vector i between iterations k and k + 1, normalized by the largest velocity in the model, is smaller than a stated threshold value ϵ, In the results presented in this article ϵ = 0.01. Note that this convergence criterion is slightly different than the version used in Braun et al. (2008).
If the elemental effective stress below the plastic yield strength, the initial material viscosity (Figure 3c) is used.
Appendix 2: Approximation of the Subduction Dip Angle For a slab of a parabolic geometry, the subduction angle changes along the length of the slab. For example, the elevation of the upper surface of the parabolic slab in our experiments is defined as where x is the distance along the x axis, z is the distance along the z axis, a = (z min − z max )/(x end − x start ) 2 , c = z max , and d = x start . In the context of our models, z min is the elevation of the base of the model ("base" in Figure A2.1), z max is the elevation of the free surface at sea level ("vertex" in Figure A2.1), x end is the distance along the x axis where the parabola intersects the base of the model ("base" in Figure A2.1), and x start is the point along the x axis where the parabola begins ("vertex" in Figure A2.1). In order to reference different subduction geometries with a single dip value, subduction angles were approximated using a triangular equivalent of the subducting slab's geometry. The angle is that between the horizontal base and the line connecting the vertex (the highest point) of the parabola at the plate margin and the point where the slab exits at the base of the model is defining the subduction dip angle ( Figure A2.1). The effective angle is in reality smaller at the margin and larger at the bottom of the model domain, but for a parabolic geometry it is difficult to determine and to generalize for easy comparison of various experiments.

Appendix 3: Strength Approximation of the Continental Crust
As described in the main text, the continental crust weak zone has been approximated with a material friction angle of 2°instead of 15°applied in the homogeneously strong continental crust. This weakness could reflect lower strength in a suture zone, using a similar friction angle to past studies including a suture (e.g., Whipp et al., 2014), or may result from thermal weakening of the crust as a result of arc magmatism. Here we roughly estimate the magnitude of the effect of thermal weakening on the crust by comparing the integrated strength of the crust in three different scenarios ( Figure A3.1): typical continental crust with no thermal weakening (15°friction angle, blue line in Figure A3.1), weak continental crust with a low friction angle and no thermal weakening (2°friction angle, red line in Figure A3.1), and thermally weakened continental  Strength profile for uppermost 35 km of a strong (15°friction angle, blue), a weak (2°friction angle, red), and a mechanically strong but thermally weakened (15°friction angle and 15-km-thick "magma chamber," green). The strength ratio shown in Figure 9 is calculated across this depth interval. crust (15°friction angle with a 900°C magma chamber; green line in Figure A3.1). A magma chamber of certain size is realized by setting the temperature to 900°C in the respective depth interval. The size of the magma chamber refers to its depth extent ( Figure A3.1); the hot region extends along the length of the oblique segment and has a 20-km width perpendicular to the trend of the oblique segment. Within the region strongly affected by the magma chamber, the viscous strength of the rock is below that for frictional sliding even for a friction angle of 2°(red and green lines in Figure A3.1), reducing the integrated strength of the crust (area left of curves in Figure A3.1) within that region. Figure A3.2 shows the temperature implementation for the thermally weakened case. The undisturbed geotherm has been calculated analytically following Stüwe (2007).