Dynamics of Oceanic Slab Tearing During Transform‐Fault Horizontally‐Oblique Subduction: Insights From 3D Numerical Modeling

Oceanic‐plates vertical tearing is seismically identified in the present‐day Earth. This type of plate tearing is frequently reported in horizontally‐oblique subduction zones where transform‐faulted oceanic plates are subducting (or subducted). However, the mechanisms behind vertical slab tearing are still poorly understood, thus we utilize 3D time‐dependent Stokes' flow thermo‐mechanical numerical models to further study this problem. We find that (a) the age offset of transform fault and (b) the horizontal obliqueness of subduction fundamentally control the tearing behavior of two generic, materially homogeneous oceanic slabs separated by a low‐viscosity zone. The two slabs sequentially bend, which combined with the age‐thickness difference between slabs, causes the differential sinking of them. Based on the modeling results, well‐developed slabs vertical tearing would happen when the oblique angle of subduction is ≥30° or the age ratio of the secondly bent to firstly bent slab being ∼<0.6. Quantifying the horizontal distance‐vector between sinking slabs, we find that subduction at medium‐low horizontal‐obliqueness angles (≤40°) of young lithosphere (slabs‐average ∼15 Myr) tends to produce fault‐perpendicular tearing. Contrastingly, old‐age slabs (average ≥ 30 Myr) with medium‐large obliqueness angles (∼>20°) tend to produce fault‐parallel tearing, related to differential slab‐hinge retreat or rollback. Correlations between slabs' (a) computed tearing horizontal‐width and (b) scaling‐theory forms of their subduction‐velocity differences, are reasonable (0.76–0.97). Our numerically predicted scenarios are reasonably consistent with plate‐tear imaging results from at least four natural subduction zones. Our modeling also suggests that continual along‐trench variation in subduction dip angle may be related to a special case of oblique subduction.

Plate tearing can be categorized into horizontal and vertical tearing patterns (Rosenbaum et al., 2008). The horizontal tearing normally occurs along a sub-horizontal plane, in which the lower segment of the slab detaches from the upper and then sinks into the mantle. In this case, the negative buoyancy of the subducted slab contributes to the slab break-off (tearing) process (Duretz et al., 2011;Gerya, 2022;Hildebrand et al., 2018;Kufner et al., 2017;L. Liu et al., 2021). On the other hand, it is generally considered difficult for a coherent subducting slab to tear along sub-vertical planes perpendicular to the trench (Cui & Li, 2022) due to the lack of major force gradients in the along-strike direction. However, according to studies of seismic tomography, vertical slab tearing (sub-vertically-oriented anelastic (brittle-visco-plastic) tears, that may be propagating in the subducted-subducting portions of oceanic slabs) seems to have developed in many subduction zones of present-day Earth (Figure 1), including the Mariana subduction zone (Miller et al., 2004;Miller et al., 2006), the Aegean and Anatolia regions (Govers & Fichtner, 2016;Jolivet et al., 2015), the Central America subduction zone (Carciumaru et al., 2020;Dougherty & Clayton, 2014;Stubailo et al., 2012;Thorkelson, 1996), the South American subduction zone (Pesicek et al., 2012;Rosenbaum et al., 2008;Vargas & Mann, 2013), the Caribbean subduction zone (Meighan et al., 2013), and the Kamchatka subduction zone (Levin et al., 2002).
Slab quasi-vertically-oriented tears (and frequently sub-vertically propagating) within the mantle generally require along-strike physical variations (such as gradients in morphology, forces, material properties, etc.), with several mechanisms being proposed: non-uniform rollback due to lateral temperature and age differences in subducting plates or step faulting at the edges of active margins (Burkett & Billen, 2010;Govers & Wortel, 2005), opposite rotation of slab segments (Gianni et al., 2019), oceanic/continental plate transition along-trench (Li et al., 2013;Magni et al., 2014;Pusok et al., 2018), subduction of spreading ridges (Burkett & Billen, 2010), and perhaps finally, variable orientations of passive margins (Fernández-García et al., 2019). Therefore, various factors can affect the specific subduction mode, with slabs undergoing translation, rotation, and deformation during the sinking, with tearing being an extreme case of slab deformation involving topological changes. Due to this complexity, in this study, we would narrow our focus to the scenarios where the oblique subduction of transform-faulted slabs develops.
As proposed in previous studies, due to the contrasting seafloor ages (and mechanical properties) across a transform fault, the vertical tearing of transform-faulted oceanic plate likely happens (Burkett & Billen, 2010). In addition, we consider the horizontal obliqueness of the subduction direction with respect to the trench as playing a key role in combination with the transform-fault age contrasts.
In this study, by systematically conducting 3D time-dependent numerical modeling to investigate the effects of the two factors for vertical oceanic-plate tearing, we will summarize the key conditions for different observed patterns for vertical slab (Figure 1 and Table 1) tearing and use a scaling method to reveal the controlling variables behind them. We hope this work will improve our understanding of these "abstract" systems, but also our understanding of real-Earth systems that comprise geological/physical ingredients of this kind; and to help us for example, when attempting to infer the pre-subduction obliqueness angle of a subducted plate that may lack hotspot-track constraints, and also to improve the accuracy of plate reconstructions.

Model Setup
We use the finite element software ASPECT (version 2.3.0) to solve the incompressible Stokes' viscous flow Equations 1 and 2 coupled with thermal advection-diffusion Equation 3 (Bangerth et al., 2021;Glerum et al., 2018). Equation 4 describes the evolution of additional (compositional) fields transported along with the velocity field u, without diffusion. For the Stokes' flow equation, we employ the Taylor-Hood spatial element.  (2013); Alaska (ala): X. Yang and Gao (2020). The identifications of slab-1 and slab-2 are defined as follows: slab-1 is the slab that would arrive and contact the trench earlier, whereas slab-2 would do it later. The vertical color bar represents the age of the oceanic lithosphere, and the horizontal color bar is the magnitude of the angle between horizontal plate velocity and the direction perpendicular to the trench. The angle sign is coded-indicated by the numbers after the second hyphen, the first number corresponds to the left-side slab when standing above the subduction zone and looking toward the continent, and the second indicates the right-side slab. The arrows indicate the plate velocity in the seafloor spreading-aligned reference frame (Becker et al., 2015). The data on the subduction zone, seafloor age, and color bar are referred to Hayes et al. (2018), Müller et al. (2008), and Crameri (2021), respectively. (b) Cartoons illustrate the two typical slab vertical tearing patterns (Rosenbaum et al., 2008 where is the ith compositional field that defines and tracks subsets of distinct chemical compositions and properties is the linear strain rate; η, ρ, C p , and k are the fluid's shear viscosity, mass-density, specific isobaric heat capacity, and thermal conductivity, respectively; and g is the gravitational acceleration, p is the pressure, T is the temperature. A 3D Cartesian model box with spatial dimensions of 2,500 km × 3,000 km × 660 km is used. Figure 2 shows the initial model configuration demonstrated in composition, including a partially subducted oceanic lithosphere: The slabs convergence angle and lithospheric ages use the averaged values within 300 km of both sides of the corresponding symbol sites on the map ( Figure 1). b Other types of subduction tectonic environments are distinguished by green stars (Figure 1), for example, ridge subduction.  The oceanic lithosphere's weak central region (transform fault zone) is 30 km wide and separates the plate in two distinct slabs. The two fundamental coordinate systems are shown, with trench -perpendicular (X) and -parallel (Y) axes, and transform-fault -parallel (TF//) and -perpendicular (TF⊥) axes. (c) The initial composition of the oceanic lithosphere (blue) and vertical slices of density (red). (d) The vertical cross-sections on the lateral edges of subducting lithosphere and adjacent plate slice lithosphere. The initial dip angle of the subducting slab is α = 30° (vertical angle, same value for all the models). The horizontal angle θ parameterizes the horizontally-oblique subduction. The length of the shortest subducted slab segment is 100 km (in slab-2) and increases gradually to the other side (slab-1).
XIN ET AL.
10.1029/2023JB027041 5 of 23 two slabs separated by a weak transform fault-zone in the middle, a continental lithosphere as the overriding plate, and two side plates ( Figure 2a). Although these denominations may conflict with previous studies, this will be, for the sake of cleanliness, our formal language-terminology from this point onward in this article: our model has one single oceanic plate or lithosphere, separated in two slabs by the transform-fault zone.
The width of the weak transform fault-zone is 30 km, and its thickness is consistent with that of the thicker lithosphere on both sides. A low-viscosity oceanic crust defines the uppermost part of the subducting lithosphere (slabs), and acts as a weak contact zone with the overriding continental plate. The initial vertical dip angle of the subducting slabs is fixed at 30°. The subducting oceanic plate is decoupled from the side plates by applying two side faults, as shown in Figure 2b. In contrast to the central plate, which is 10 km from the left boundary and undergoes free subduction, no subducting slab is prescribed for the side plates ( Figure 2d).
The half-space cooling model is used to calculate the oceanic plate's initial-condition thermal-age thickness as ≈ 2.32 √ (h max = 100 km) (Turcotte & Schubert, 2002), where κ is thermal diffusivity, and t is seafloor age. In this sense, the thickness ratio between slabs 1 and 2, that is, h 1 /h 2 , is dynamically similar to the square root of the ratio of lithospheric ages √ 1∕ 2 . The ocean plate's initial (unsubducted) maximum thickness is limited to 100 km, consistent with thermo-rheological observations (Schubert et al., 2001). The model-box's top and bottom boundary temperatures are fixed at 293 and 1775 K (calculated using the mantle adiabatic thermal gradient), respectively.
Our representation assumes the oceanic-plate as a numerical compositional field, with two distinct variable-thickness domains, layers or slabs whose initial-condition upper and lower surfaces (and thicknesses) are defined by two isothermal surfaces of the half-space cooling model (whose maximum depth depends on the horizontal transform-fault-parallel coordinate as previously indicated). The slabs' density and viscosity are assumed to remain constant (this is a major simplification: these model-variables do not depend on pressure, temperature, or strain rate) with well-established values (e.g., S. Yang et al., 2021) (Table 2); thermal expansion coefficient is assumed null and thermal conductivity is fixed. This simplification is only a first-order approximation which allows faster computations. In the mantle domain, an adiabatic temperature profile is applied as initial condition.
Our model slabs can change their thickness as they evolve in time, but solely by virtue of divergence-free flow advection and deformation (they can thicken/thin if their perpendicular area decreases/increases), a dynamical consequence of 3D-time-dependent velocity field. The thermal field spatio-temporally evolves and its coldest regions fairly approximate the compositional slabs (especially at early times given slow-enough thermal diffusion), and this provides additional supporting information on the system (later on used in two of our figures).
The velocity field boundary conditions (at all times) are of tangential flow boundary, that is, free-slip (zero shear stress) and no-flow-through (zero normal velocity), on all the model-box surface boundaries except for that at the bottom: the lower boundary is fixed (3D zero velocity vector) and intended to simulate the effect of the 660 km phase-transition interface.
Finally, the computational grid is adaptively refined five times until its smaller cell reaches a minimum size of 6.51 km × 6.25 km × 6.88 km, and the coarsest cell size is 8 times larger. The maximum number of cells is ∼3 × 10 5 and the number of unknowns is ∼8 × 10 7 . Usually, it takes 24-48 hr for each model, running with 220 computing cores, to simulate a geological time of about 30 Myr.

Results
A total of 69 experiments/simulations are conducted with different relevant combinations of oceanic lithosphere's (a) horizontal-obliqueness angles and (b) slabs age offset. To quantify the size of the tear width, we visualized and post-processed the output data in ParaView, specifically by selecting horizontal slices of the oceanic lithosphere at various times and depths and measuring the horizontal vector distance between the two slabs ( Figure 3).

Tested Conditions for Vertical Slab Tearing
In our simulations, the slabs' tearing process can be separated, before and after, by the event of slab 1 reaching the 660 km depth (as mapped by the oceanic-lithosphere compositional field (in ASPECT software) value >0.5 within 15 km from such depth). In particular, for our reference model (t 1 = 33 Myr, t 2 = 27 Myr, θ = 20°; , with the variable "d" being the time-dependent horizontal vector distance (of magnitude "width") of slabs' tearing window. The vector-components ratio "dx/dy" and the oblique-convergence angle θ are used to characterize the slabs' tearing patterns as explained in Section 3.2 and 4.3.
7 of 23 abbreviated as 3327_20), we show snapshots of its evolution in Figure 3a by isolating temperatures ≤1380 K, so the sinking slabs appear bounded by isothermal surfaces (with the hotter ambient mantle made invisible).
Using ASPECT-software compositional fields, for all the simulations, we isolate the slabs from the embedding ambient mantle (as for the reference simulation in Figure 3b), and then we measure the tear horizontal-widths at different depths in both time-stages (as in Figure 3b).
As can be seen from the white contour lines (depths: 100-200-300-400 km) in Figure 3a, vertical tearing can occur at deep locations (>200 km depth). As the subduction progresses, the finite-width tear gradually (after 12 Myr) expands and vertically propagates to the shallow part (∼100 km). After slab 1's tip arrives at the box bottom, its sinking velocity decreases, contributing to the increase in tearing width ( Figure 3) at intermediate depths.
After identifying the 660 km-depth contact event of slab 1, we select two sampling times: 2 Myr before and 2 Myr after the event. At each one of those times we evaluate the geometry, we take the maximum value of tearing width d max over the different depths for the subsequent determination of the tearing conditions (as in Figure 3b). Figure 4 displays the tearing judgments/decision array of the full set of numerical experiments.
Each one of the three major rows on the diagram in Figure 4 summarizes the results of simulations with the same arithmetic-mean age 1 + 2 2 of the plates separated by the transform fault. To classify each simulation, we use relative age ratio rather than age difference as a way to avoid discussing parameters that deviate far from nature, that , For the same age ratio, a younger mean age suggests a smaller absolute age difference.
For any given mean age, significant oceanic-plate tearing is favored by a large obliqueness angle and/or by older slab-1 age relative to slab-2. Most models fail to tear without horizontal obliqueness of the subduction (angle θ = 0°) except for medium-aged plates with substantial age differences.
In particular, due to thin slabs, Model 1 (1318_20) in the low mean absolute age group (15 Myr) exhibits lateral break-off, increasing tearing width. At 60 Myr-30° and 60 Myr-40°, the tearing width does not increase beyond a certain extent before slab 1 reaches the bottom, so its tearing width is small during stage 1. In cross-sectional comparison, we find that the tearing area is larger in stage 2 than in stage 1, especially the 60 Myr group changes from transition to dark green color code. In the vertical comparison, we found that the tear area increases and then decreases with increasing mean age, which will be further analyzed in the subsequent discussion.
In conclusion, our simulations suggest that for any values of obliquity, the subducting plates tear for age ratios less than 0.6, meaning that an older-thicker slab-1 entering the trench first, then sinking first and faster, progressively tears off from slab-2. Furthermore, regardless of the transform-fault age ratios, all the plates tear at horizontal obliqueness greater than 30°. This latter behavior is certainly an effect of the successive bending of an obliquely subducting slab-1 first and then later on slab-2. The larger the horizontal-obliqueness angle, the larger the difference between the slabs' subducted areas (especially during the early stages of the process) and the more intensive and extensive the strain of the weak fault zone between the slabs. Finally, as crucial reference cases, we must note the case when the horizontal-obliqueness θ = 0° and simultaneously slabs' age-ratio = 1, meaning two identical slabs separated by a weak zone and converging perpendicular to the trench. Those simulations show tearing displacements that when existing and being measurable at all, have values d max < 50-75 km, meaning that tearing almost does not occur (symbolized as triangles, in dark red and transition in Figure 4).

Patterns of Vertical Slab Tearing
In this study, we consider the horizontal distance vector of tearing displacement between the slabs: d = d x e x + d y e y , measured at depth in their vertical portions inside the ambient mantle.
The maximum tearing width (displacement magnitude) over different depths of the subducting plate (depths: 100-200-300-400 km, from slice-contours as shown in Figure 3 and subsequent simulations' snapshot figures) is used to indicate whether tearing occurs (summarized in Figure 4). In addition, the ratio of the tearing vector components in the transform-fault-based system (given an obliqueness angle θ, see Figures 2 and 3 for definitions) is employed to characterize the tearing patterns. We project the tearing distance-vector components d x and d y (referred to trench perpendicular-and-parallel coordinates, respectively) onto the initial horizontally-obliqueconvergence direction and its normal direction (i.e., onto the coordinate system defined by the transform-fault 9 of 23 direction and its perpendicular), respectively, to obtain the fault-parallel fault∕∕ and fault-perpendicular-fault⟂ tearing components. We obtain: 3, the tearing pattern is defined as fault-parallel tearing (related to slabs' differential hinge retreat or rollback) ( Figure 5). In contrast with | | | Plate-motion horizontal obliquity influences the fault-parallel and also the transverse/fault-perpendicular direction of forces and resulting tearing displacements. However, it seems dominant only in young slabs, such as the region less than 40° in the 15 Myr group in Figure 5. On the other hand, large absolute age differences generally result in large sinking velocity differences that project in the convergence direction, such as the region θ ∼>20° and t 2 /t 1 ∼< 0.6 (bottom left corner) in the 60 Myr group in Figure 5.
As shown by Figure 6, considering older slabs at the initial condition, with mean ages from 15 to 30 and to 60 Myr, fault-perpendicular tearing patterns tend to disappear (occurring only in some cases), meaning that it is mostly a young-age tearing behavior. On the other hand, fault-parallel manifests commonly and extensively only for old slabs with a mean age of 60 Myr, with only one exception at a mean age of 30 Myr. Figure 7 shows the evolution of the typical models of the two tearing patterns, simulation 1713_30 for fault-perpendicular tearing and simulation 6753_30 for fault-parallel tearing. The larger mean age/thickness enhances the differential retreat of the two slabs in the convergence direction, causing the vertical tearing (while inside the mantle) pattern to be dominated by fault-parallel tearing (Figures 7a and 7d). In fact, at old slab ages (such as 60 Myr), the same age ratios as in the other age groups are attained with larger absolute-age differences, and these latter ones imply larger absolute sinking-and-subduction velocity differences (see Section 4.1). Therefore, old-age slabs tend to favor differential slabs' retreat and fault-parallel tearing ( Figure 5). Additionally, the distance from the trench to the 90° slab-dip angle position, called the bending length l b , varies between the two differentially retreating slabs and with time evolution, and this will be further analyzed in the subsequent discussion section (Figures 7b and 7c).

Analytical Study of the Slab Tearing Dynamics
In order to understand how the various controlling factors operate in our simulations and to enhance the verifiability of the model results, we utilize scaling similarity theory and available analytical-numerical studies.
We hypothesize that the horizontal width (d) of the tearing window is correlated to the difference in subduction velocity Δv between the two slabs separated by a transform fault. In a 2D subduction model, the oceanic slab is primarily subjected to (a) downward-buoyancy force δρ·g·h·l per unit length (with δρ density anomaly, g gravity . Check-box diagrams for the models' tested conditions. Angle means the horizontal obliqueness angle in the systems and defines the small columns inside each panel. The arithmetic-mean lithospheric age of slabs 1 and 2, ((t 1 + t 2 )/2), is (a) 15, (b) 30, and (c) 60 Myr, respectively, defining the big panel-rows of the array plot; whereas the ratio t 2 /t 1 defines the small rows inside panels. A slab reaching the 660 km-depth event defines the big panel-columns of the array plot: stage 1 (before) and stage 2 (after). Each element on one side panel and its corresponding on the other side panel represent one single simulation. The category of slab-tearing discrimination: Colors indicate the geometric state of slab tearing at each particular time-stage (at 1 (left) and at 2 (right)), as Red tones: maximum tearing horizontal width d max < 50 km, green tones: d max > 100 km, and transition d max : 50-100 km, where the maximum values are taken among the different sampling depths (100-200-300-400 km). Symbol shapes indicate the "overall" situation integrated through both time-stages: If the tearing width remains <70 km throughout all time stages and additionally becomes <50 km at least at one stage, the model is classified as "no tearing" (triangles). If the tearing width >100 km at any of the stages the model is classified as "tearing" (circles). Cases in between are transitional modes (rectangles). These symbol shapes are necessary for the subsequent figures of this article.
XIN ET AL.
10.1029/2023JB027041 10 of 23 acceleration, h the oceanic plate thickness, and l the slab tip-length along convergence direction), (b) surrounding mantle resistance F 1 along the tangential direction of the slab η mantle ·v (where v is the plate speed, η mantle the ambient-mantle viscosity), and (c) unbending resistance F 2 perpendicular to the slab η plate ·v·h 3 /l b 3 , where η plate is the oceanic-plate viscosity, l b is the bending length from the slab tip to the unbending place behind the hinge point, similar to the previous description in (Li & Ribe, 2012;Ribe, 2001). Then we take the one-slab tangential velocity 2 = ⋅ mantle sin ⋅ ℎ ⋅ (in an idealized, terminal steady-state scenario where the slab, driven by the downward-buoyancy force, is "gliding" inside the viscous mantle) along the slab subduction direction, as the major component that defines the tangential velocity difference Δv between slabs, yielding where α is the subduction dip (vertical) angle.
We assume that at the initial-condition of the subduction zone, in each one of the two slabs subducted portions (Figure 2), their thicknesses h 1 and h 2 remain constant. Then in this initial configuration and due to horizontal obliqueness, the initial slabs subducted-length l varies along the trench direction (y-direction) following ( ) = + ⋅tan cos , with min = cos , max = +2 ⋅tan cos , where c is fixed at 100 km ( Figure 2b). Then we get 3 = ∫ 0 ∕ ∫ 0 1 , as the spatial-average velocity of the entire slab tip, with w the horizontal span of the slab along the trench axis.
Then, the simple-theoretical spatial-average velocity difference between slabs is obtained as (the unit is s −1 ) is a functional factor composed of several dynamic parameters kept constant in each model simulation for simplicity, b = w/cosθ is the fixed slab width (Figure 2b), and the functional factor = 2 (√ sin (cos ) 2 considers the geometrical parameters reflecting the horizontally-oblique subduction effects, with the ratios of the square root of time representing ratios of thermal-age thicknesses, as ≈ 2.32 √ gives each oceanic slab thickness h versus its seafloor age t.
Following the logic of the unbending-resistance form η plate ·v·h 3 /l b 3 , in order to find the correlation between the model-output computed tearing distance d and the theoretical underlying controlling factors for the velocity difference Δv, we consider an additional multiplying variable factor f in the former simple-theoretical expression. Figure 5. Check-box diagrams for the tearing patterns on time-stages 1 (left) and 2 (right). Angle means the horizontal-obliqueness angle. The arithmetic-mean lithospheric age of slabs 1 and 2, ((t 1 + t 2 )/2), is (a) 15, (b) 30, and (c) 60 Myr, respectively. The category of slabs tearing pattern, depending on the ratio " | fault∕∕∕ fault⟂| ." Orange: If | fault∕∕∕ fault⟂| < 0.7, the pattern is defined as fault-perpendicular tearing. Blue: In contrast with | fault∕∕∕ fault⟂| > 1.3, the pattern is defined as fault-parallel tearing. Gray: A transition-hybrid mode is 0.7-1.3. For each simulation and at each time stage, the tearing components ratio is averaged over the different sampling depths (100-200-300-400 km). The definition of symbol shapes is the same as in Figure 4. Each age-group stage-2 simulations enclosed by dashed-line squares are used in Figures 6 and 7.   Figure 6. Example models illustrating the three different slab tearing patterns (determined by the ratio of tearing horizontal-displacement components). End-member cartoons are provided for illustration. These three particular example simulations (signaled by dashed squares in Figures 5a-5c), while differing only in their slabs' mean age, all have age ratios of about ∼0.8 and obliqueness-angle θ = 30°. In the three example-model plots, the compositional fields of slabs and transform-fault zone are colored in light gray-blue and in dark blue, respectively.
12 of 23 f is a dimensionless parameter representing additional effects on Δv, such as unbending resistance along the slab-normal direction, differences between the two slabs' bending-sinking lengths and thicknesses, viscosity and width of the weak zone (transform fault), etc.; all of which depend on time as the sinking process evolves.
Therefore the original dimensional expression k·C·h is expanded, upgraded and refined to k·C·h·f. In this way, the dimensionless parameter f should accommodate for the aforementioned additional effects, and we use trial expressions of the form ≈ ( ℎ 2 ℎ 1 ) 1 2 to represent and account for them.
We used a set of fixed times (10, 14, 20 Myr) and depth range (100-400 km) to compare and correlate the tearing horizontal-width d with the spatial-average velocity difference Δv between slabs (brought by analytical scaling theory), and the results are shown in Figures 8-11. The figures are plotted for respective absolute mean ages of plates at {15, 30, 60} Myr and all of them together.
Upon constructing the correlations shown in Figures 8-11 and Table 3, we have removed four simulations from the data set, specifically simulations 1318_30, 1318_40, 1515_20, and 1515_30. This was done because these particular young "ultra-thin" slabs broke off during the subduction process, so we cannot precisely measure the tearing width in these models. Simulation 4317_40 can greatly increase the correlation coefficient in a way we distrust, so we also discarded it.
After a detailed analysis of our numerical experiments (Figures 8-10), we find that the tearing horizontal width d is sensibly and roughly proportional to the slabs' velocity difference. Our hypothesis is reasonably confirmed. If we consider the "full" scaling form for the velocity difference Δv = k·C·h·(h 2 /h 1 ) n , the correlation coefficient R 2 is in the range of 0.76-0.97 (see Table 3) considering all depths and times, age groups and simulation cases. For comparison, if employing the much simpler "incomplete" scaling form Δv = k·C, the range of R 2 is 0.69-0.96 for the same simulations cases as in the previous trial (see Table 3). Therefore, regarding the similarity between the two correlation-values intervals, it is inferred that the first-order governing factor is the oblique convergence given by the characteristic parameter C = C(θ), as defined in Equation 6. Figure 5. Panels (a, d): View from above, thermal contour (yellow, the legend is the same as Figure 3) and viscosity (green) structure (Table 2) Table 2).

Figure 7. Evolution of the two case simulations that exemplify the two tearing patterns, singled out and indicated by dashed-line squares in
XIN ET AL.

10.1029/2023JB027041
13 of 23 In particular, if we divide each group of time-depth diagrams along the diagonal line (Figures 8-10), the parameter changes can roughly reflect the tearing width evolution in different segments of the ocean plate. Most parameters remain consistent within each one of the areas delineated by the red slash, reflecting the previous deformation history inherited by the slab during most of the subduction. The variation of parameters in various intervals from the lower left to the upper right depicts the variation of governing factors with subduction, which is dependent on (h 2 /h 1 ) n , like n = 0 → −1 in the 15 Myr-age group and n = 1 → 0 → −1 in the 60 Myr-age group (Figures 8b  and 10b), and n = 0 → −1 at 300-400 km depth in the 30 Myr-age group (Figure 9b).
This overall reduced evolution represents the effect of unbending resistance F 2 in various plate segments. The oceanic lithosphere subducted length l of slab-2 is obviously smaller than that of slab-1 during the early stages of subduction (well before the two slabs are dipping sub-vertically). The smaller l in slab-2, the stronger the unbending resistance and the apparently "stiffer" the slab. Accordingly, "stiffer" slab-2 affects the change of tearing width when h = h 2 , n = 1. The overall subduction dip angle (including that of slab-1) increases as slab-2 subducts, and the effects of F 2 grow weaker at this point with h = h mean and n = 0. Later on, with continued slabs sinking, there is the singular event of slab-1 reaching 660 km depth: slab-1 is subjected to a force at the bottom that resists bending. At this point, the velocity-difference Δv turns to be controlled by slab-1, with the exponent transitioning to n = −1. As for the variation of the specific value of n in different age groups (15, 30, 60 Myr), we believe it is f that is also correlated with l b1 /l b2 . The changes in bending length compensate for some of the changes in slabs' thickness ratio.
Finally, it can be seen that the controlling parameters change in the slabs deepest regions (400 km) on stage 2 (20 Myr) (Figures 9b and 10b). We think it is slab-1 here that is bent at the bottom of the model, that is associated with the folding pattern in the 30 Myr-age group and the trench retreating pattern in the 60 Myr-age group (Figures 9a and 10a) (Li & Ribe, 2012).
The comparison among the models with different mean slab ages indicates that the older slabs of 60 Myr generally have a narrower tearing window than the mid-age slabs with mean ages of 30 Myr (Figures 4 and 11).
Certainly, an additional multiplication factor (h 30 /h 60 ) p must be applied in the 60 Myr mean-age simulations to fit in with the other age groups at similar correlation levels (Figure 11). At old-plate ages (mean 60 Myr), the full scaling form Δv = k·C·h·f (with f = (h 2 /h 1 ) n ) yields values that are too big (because h is big). Therefore, an additional post-factor <1 is needed to fit the observed values of d through the range of Δv values of the different simulations and to get similar correlation levels to those of the other age groups.
Beyond this, we believe that the smaller values of tearing width in the 60 Myr than in the 30 Myr age group may be related to the difference in bending length (to measure the effect of bending length (l b1 /l b2 ), we use the average slab-thickness ratio (h 30 /h 60 ) p , which can be compared quantitatively) involving viscous dissipation, with older-thicker slabs offering higher viscous resistance to bending and retarding the development of the bending-length increase. The sequence of slabs-1, 2 changing from small to large bending length is different when plate age-thicknesses are large ( Figure 6). The apparently rheologically weaker slab-1 changes "earlier" in the plate-thickness sinking-velocity parameter space than slab-2, increasing the tearing width from 15 Myr-age groups to 30 Myr-age groups and then decreasing (just to some extent) from 30 Myr-age groups to 60 Myr-age groups. This effect is more apparent in the later phases of subduction and at larger depths, as shown by the change of p in Figure 11. This feature is also shown in the shift of the age ratio cut-off with absolute age in Figure 4.
In summary, vertical oceanic-plate tearing is controlled by two dominant factors, (a) the subduction horizontal obliqueness: the larger the obliqueness angle, the larger the difference in the subducted area of the two slabs, the higher the likelihood for slabs tearing; (b) the difference between the ages of the two slabs: tearing across/along the transform-fault zone is favored by older, heavier, thicker slab-1 relative to slab-2, as this configuration plays constructively with the obliqueness. The plate unbending resistance effect varies among the individual slabs thicknesses and also among the different subduction zones (in real Earth), but seems to only play a secondary role in tearing generation and development.

On Transform-Fault-Perpendicular Tearing
In Section 3.2 we showed that tearing horizontal-displacement-components ratio, when expressed in the transform-fault-based coordinate frame can take diverse values far-from yet about 1. If dominated by fault-parallel component, we recall that slabs differential roll-back and differential hinge retreat is a relatively normal phenomenon in geosciences. Contrastingly, if dominated by fault-perpendicular component the slabs may be slightly separating inside the asthenosphere, parallel to their transversal/hinge axes. We will shed light on this hitherto yet-to-be-understood phenomenon.
We already provided kinematical (visual displays of geometry vs. time) evidence that this transform-fault-perpendicular tearing displacement can dominate the tearing, and that it manifests only for horizontal-obliqueness θ ≤ 40° and more frequently in young slabs (Section 3.2, Figure 5). We shall now show some supporting dynamical evidence for this tearing behavior.
Let's consider X' as the fault-parallel axis and Y' the fault-perpendicular axis. Figure S1 in Supporting Information S1 displays the viscous horizontal shear stress on the slabs surface for simulation 1713_30 (lithospheric ages t 1 = 17 Myr, t 2 = 13 Myr; horizontal obliqueness θ = 30°) which exemplifies transform-fault-perpendicular tearing. The viscous horizontal shear stress   is shown at 6 and 12 Myr of computed geological time. Given the horizontal-obliqueness angle θ = 30°, the relative motion between the slabs and the overriding continental plate, particularly the tangential/trench-parallel component, influences the local stress field on the surface of these bodies: the surface viscous friction imprinted by the relative motion dictates the shear stress. The effects of this shearing may not be equal on the two slabs regarding their different ages and subducted volumes versus time, and the interaction with the overriding plate may even be unsteady in time, resulting in an increasing separation between the slabs. Figure S2 in Supporting Information S1 (in its left column) shows the viscous horizontal normal deviatoric stress on a 2D vertical slice along the transform-fault-zone central-axis for simulation 1713_30 exemplifying transform-fault-perpendicular tearing. The stress   is shown at {4, 6, 8, 10, 12, 14} Myr of computed geological time. The normal deviatoric stress   (left-side panel) inside the transform-fault-zone is extensive in-and-around the subduction hinge (bent), and slightly compressive where the slabs are close to each other far from the trench in the horizontal portions; and this stress pattern changes only slightly during this particular subduction process. With time, the fault-zone extensive stresses slightly oscillate in magnitude, meanwhile and incrementally, the slabs separate horizontally inside the ambient mantle (right-side panel).

On the Trench-Based Alternative Tearing Patterns
In this study we opted to express the tearing horizontal-displacement vector in reference frames defined by the transform-fault-based system in each simulation, case by case depending on the horizontal-obliqueness angle. An interesting alternative is to use a universal (for all the simulations) reference frame, the one that defines the modelbox: we shall now express the tearing in the "universal" trench-based system, with axes X (trench-perpendicular) and Y (trench-parallel) as in Figure 2. This is important as the overriding continental plate with its fixed trench, is a common single entity which all the different simulated systems have and interact with.
In this "universal" reference frame, the tearing displacement (d x ,d y ) is used to directly compute the magnitude of the tearing horizontal-components ratio |d x /d y | to characterize the tearing patterns. The values of this ratio tend to be <1 due to the horizontal obliqueness: when the angle θ > 0° the "d y " component tends to be non-zero in most cases, but "d x " can be about zero in several cases. When the ratio is <0.5 trench-parallel tearing dominates thoroughly and more than half of the simulations are represented, and when >1 the pattern is trench-perpendicular-dominated but only 1 (in stage 1) and 4 (in stage 2) out of 69 cases were found.
Particularly, when the subduction horizontal-obliqueness is zero, the two coordinate-frames coincide, and fault-perpendicular equals trench-parallel. Indeed, observing Figures 5 and 12, we can confirm that for angle θ = 0° all the simulations colored in orange in Figure 5 have a ratio <0.7 colored correspondingly in Figure 12, and the only simulation colored blue (ratio >1.3) in Figure 5 is correspondingly colored as >0.9 in Figure 12. Figure 11. The nine panels display the correlation between slabs tearing width (d) and the integrated controlling factor (Δv/k = C·h·f) for all models. The correlation coefficients (R 2 ) are shown in the top left corner of each panel. We find that an additional multiplication factor (h 30 /h 60 ) p must be applied in the 60 Myr mean-age simulations for the dv = k·C·h·f function-group to capture a certain influence of slabs thickness differences on tearing width, and in this way they fit in reasonably with the other age groups.
XIN ET AL.

10.1029/2023JB027041
17 of 23 For obliqueness θ > 0°, the relationship between the patterns from the two coordinate-frames-descriptions is non-trivial, specifically it is non-linear and given by trigonometric projections.
The overall tearing patterns in the trench-based system are the following: 1. Trench-parallel tearing is mostly a young-slabs at high horizontal-obliqueness phenomenon: There is a clear dominance of ratio values <0.5 (strong trench-parallel tearing displacement) for young-intermediate (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30) ages, and being more prevalent for high obliqueness angles but virtually independent of transform-fault age ratios. Aside from this, for young slabs this pattern accentuates with temporal evolution, more prevalent and intense in stage 2 after one slab reaches the bottom; and only slightly more prevalent in late stage for the 60 Myr age group; however, the opposite occurs for mid-age slabs (pattern is more prevalent at early stage). 2. The trench-perpendicular tearing-displacement component, being relatively smaller in general, becomes significant only for intermediate-old (30-60 Myr) ages and mostly low horizontal-obliqueness (for the transform-fault-parallel tearing to project effectively on that axis). The components-ratio becomes comparable or larger than unit (ratio ≳ 1) only in a few simulations ( Figure 12). Also, trench-perpendicular tearing component becomes (excepting for two simulations) more important (prevalent and intense) in late times in stage 2, as an effect of one slab reaching the bottom and blocking the motion for some time.
The evident dominance and prevalence of trench-parallel tearing (ratios < 0.5), being common for young slabs' mean ages and high obliqueness angles, is interpreted as: Young ages signify smaller sinking-velocity differences (constructively combined with older slab-2 preventing excessive vertical separation, in the 15 Myr group only) combined with high obliqueness angles that signify higher trench-parallel relative velocities between slabs and trench, thus higher viscous trench-parallel shear stress on the contact. These conditions, of low vertical velocity-displacement difference between slabs, with high tangential forces on them, may favor the dominance of trench-parallel tearing. On the other hand, the manifestation of trench-perpendicular tearing is more complex, as it not only depends on the time stage, but perhaps on an obvious trade-off: similar constructive-tearing results for low age-ratio (relatively older slab-1) and high obliqueness angle, which can result in non-uniform (differential) rollback and smaller (delayed) bending-length of "apparently stiffer" plate 2, respectively. Finally, as seen in Figure 12, the overall temporal change of patterns between stages 1 and 2 is due to one slab reaching the bottom, and the consequential difference in sinking-subducting style that ensues, which affects the slabs tearing displacement.

Geological Application
As shown in Figure 13, examples of slab tearing of natural subduction zones are projected onto our model phase diagrams. Some cases (such as "Central America subduction zone (cam-1)" and "Alaska subduction zone (ala)") occur simultaneously in two mean-age groups, given our sparse and imprecise sampling of this variable. Except for the "Alaska (ala)" point, all the natural cases fall in our model-predicted regime where tearing occurs. The misfit of the "Alaska (ala)" case may be due to a secondary (Yakutat) plate with different seafloor ages (X. Yang & Gao, 2020). Furthermore, because the horizontal-obliqueness angle of "Caribbean (car)" subduction is ≥60°, we have not carried out corresponding simulations considering the numerical stability problems they face. According to our analysis, the Aleutian slab window's low-velocity zone emplaced on a high oblique-convergence angle may undergo tearing due to slab-1 bending (Bai et al., 2020).
We measured the tear widths observed in the literature. We projected the tearing widths onto the convergence (transform-fault-parallel) direction and its perpendicular direction, which were used to determine the corresponding   (Meighan et al., 2013); South America (sam-1): (Pesicek et al., 2012); Central America (cam-1): (Stubailo et al., 2012); Central America (cam-2): (Carciumaru et al., 2020)) ( Figure  S3 in Supporting Information S1). Colored cartoons represent the tearing patterns in nature, whereas the tearing patterns predicted by the model are represented by colors in the phase diagram ( Figure 13).
We find that fault-perpendicular tearing manifests in "Central America (cam-2)," the transitional node manifests in "Central America (cam-1)," and fault-parallel tearing manifests in "Caribbean (car)" and "Mariana (mar)." Our simulations do not adequately explain the real "South America (sam-1 and sam-2)" tearing pattern, with these real cases being likely in transition modes that, if actually tearing off, are more inclined to fault-parallel tearing as suggested by observations (not appreciable in Figure 13g cartoon), most likely due to the real plate being thin and experiencing break-off during subduction (something that may be occurring in the Nazca plate under central Chile (Beck & Zandt, 2007;Beck et al., 2000;Cahill & Isacks, 1992;Portner et al., 2020), and Ecuador-Colombia (Wagner et al., 2017)), as shown in Figure 1 and Table 1, which increases the tearing width in the along-transform-fault direction. According to our understanding, the vertical tearing pattern of "Kamchatka (kam)" is most likely fault-parallel tearing. Figure 12. Check-box diagrams for trench-based tearing patterns on time-stages 1 (left) and 2 (right). Angle means the slabs' horizontal-obliqueness angle. The arithmetic-mean lithospheric age of plates 1 and 2, ((t 1 + t 2 )/2), is (a) 15, (b) 30, and (c) 60 Myr, respectively. The trench-based coordinate frame is as shown in Figure 2: X axis (trench-perpendicular), Y axis (trench-parallel). The definition of symbol shapes is the same as in Figure 4. Black triangles are the no-tearing simulations. The category of slab trench-based tearing pattern, depending on the ratio "|d x /d y |." Red-color tones: |d x /d y | < 0.5, contrasting with Blue-color tones: |d x / d y | ≥ 0.5. The smaller the ratio |d x /d y | the stronger dominance of trench-parallel displacement component and the more extreme the trench-parallel-dominated pattern, and most simulations are represented in this group. Throughout time, only a few simulation cases attain displacement ratios ≳1 qualifying as trench-perpendiculardominated tearing (TP P DT).

Figure 13.
Comparisons between our numerical-model predictions on time-stage 2 ((a-c) phase diagrams color-coded for tearing pattern and shape-coded for overall tearing occurrence) and (d-g) natural subduction-zones observations (small black circles) along with their simplified observation-based interpreted structures (cartoons). (d) The case with fault-perpendicular tearing pattern, (e) the case with fault-parallel tearing pattern, (f) the case with transitional mode pattern, (g) the case in transitional mode that, if actually tearing off, are more inclined to fault-parallel tearing according to the observations. The "Alaska (ala)" natural case falling in the predicted no-tearing zone may be affected by a possible secondary slab tear. In contrast, the "South America (sam-1 and sam-2) differential retreat" case inclined to a fault-parallel pattern may be affected by a lateral break-off in the younger slab but also lack more detailed local observations. Colored diagram: the same as in Figure  In order to investigate whether our model can be applied to other types of subduction, such as "Mariana (mar)" and "South America (sam-2)" involving ridge subduction. We performed comparative experiments of oblique subduction with/without a ridge ( Figure S4 and S5 in Supporting Information S1). In these cases, weak zones near the front and back of the ridge may exist in the form of rift or trough, for example, in the Central American region (Morell et al., 2012;Smith et al., 2011), and New Zealand's southeast area (Grobys et al., 2007). When the ridge undergoes subduction, as long as a weak zone exists in front of the ridge, the presence of this latter does not change the vertical tearing results/pattern of the forward slab, indicating that our model can be potentially applied to cases where a weak zone exists in front of the ridge ( Figure S6 in Supporting Information S1).
In addition, we find that when the oceanic ridge is at the front edge of the plate, and the weak zone is behind, its oblique subduction results in a pattern similar to the tearing pattern of the opposite rotation of slab segments (Gianni et al., 2019) ( Figure S6 in Supporting Information S1). Details of oblique subduction of a ridge need to be studied further in dedicated future studies, for example, using numerical models with non-Newtonian rheology, strain-rate-dependent viscosity.
We also carried out additional simulations altering the viscosity of the weak transform-fault zone and found that a high-enough fault-zone viscosity can (in some cases) preclude and inhibit tearing. However, we must note that when slab-2 is not entirely subducted during the early stages of subduction evolution, the variable slab length along the trench generates transverse non-uniform slab "stiffness" (Li & Ribe, 2012). In these cases, the deformation results in a continuous variation in subduction dip angle ( Figure S7 in Supporting Information S1), like the imaging on the right side of "South America (sam-1)" (Figure 1).

Model Limitations
Our model is simplified by having constant (yet adequate) density, viscosity, and thermal properties values. This was chosen to isolate the first-order controls on the slab tearing during transform-fault subduction. This simplification, also applied in order to facilitate the quantitative scaling analysis, may lead to some inaccuracies related to the following: 1. The age difference between the two sides of the transform fault should lead to additional lateral variations in the oceanic plate's physical properties (e.g., viscosities). 2. The viscosity and density should vary continuously with depth (pressure and temperature), particularly through and between lithospheric and asthenospheric mantle regions. 3. Additional complexities and effects on the sinking-tearing plate, that is, more specific and refined forms of the post-factor function f, were not tried. 4. The obliqueness angle and age-ratio in the natural cases may have changed through time. Therefore, their values during the early stages of subduction (Myr ago) may differ from those observed at the present day.
All the aforementioned neglected factors/processes can influence the slabs' tearing vector (magnitude and direction) at all stages during subduction. However, we believe that the first-order controlling mechanisms are adequately represented in our study.

Conclusions
Subduction horizontal obliqueness and slabs' transform-fault age-offset play dominant first-order roles in the development and evolution of slabs tearing during the subduction of such complex systems. The viscous unbending resistance varies among the two subducting slabs during the subduction process, and also varies among different subduction zones, yet it seems to only have a secondary role in slab tearing.
When the subduction horizontal-obliqueness angle is ≥30° or the age ratio of slab-2 to slab-1 is less than 0.6, fast and well-developed oceanic plate vertical tearing occurs and develops inside the mantle without any other requirements. In contrast, no tearing occurs for small obliqueness angles (<20°) if slab-2 is approximately older and thicker (age ratio >0.6-0.8).
The vector nature of slab tearing can occur/be decomposed in two patterns or modes when referred to the convergence direction or transform-fault orientation: fault-parallel and fault-perpendicular tearing. Fault-perpendicular tearing manifests mostly (not solely) in young-age slabs and for medium-low subduction horizontal-obliqueness angles (≤40°) and somewhat more favorably for low-age ratios (slab-1 being older). Fault-parallel tearing (related to differential slabs-hinge retreat or rollback) is relatively well-known and manifests almost solely for old-age slabs (60 Myr-mean-age group) and mostly for medium-high obliqueness (∼>20°) angles but almost independently of age ratios.
Transform-fault-perpendicular tearing of subducting slabs' is here kinematically identified fairly/slightly above numerical errors and model uncertainties, and classified in a parametric way using the tearing horizontal-distance vector. Dynamically speaking, some of our simulations suggest that this tearing mode is accompanied by viscous horizontal normal deviatoric stresses of extension across the weak fault-zone, in-and-around the hinge (subduction bent) region, with some slight variations in time; all of this concomitant to the slabs' increasing horizontal separation inside the ambient mantle. Future research is needed to better understand this particular tearing behavior, but hereby we have provided a first step in that direction.
Compared to the several transform-fault-based systems in our simulations, the universal trench-based reference frame is defined by a single common entity (overriding plate), and when the slabs' horizontal-obliqueness angle is 0° (normal subduction) the two reference frames coincide. Trench-parallel tearing is frequent and dominant for young slabs, and may be caused by low vertical-velocity difference between slabs, and moving at high horizontal-obliqueness angles (∼>20°) causing high tangential shear stress at the trench. Trench-perpendicular tearing is relatively infrequent and likely an old-age phenomenon (∼>30 Myr), favored by mid-low obliqueness angles, and closely related to non-uniform slabs rollback.
Our modeling suggests that the continual spatial variation in dip angle along the trench observed in real-Earth subduction zones could be potentially related to the initial evolution of subducting slabs under horizontally-oblique convergence. Further beyond, knowing the age combination across a transform fault (ratio and average) and the horizontal obliqueness of an oceanic subduction scenario would, with a framework like the one here presented, allow a prediction of the tearing occurrence, potential pattern, and the likely evolution of the subducting slabs, and perhaps of the associated locally-perturbed mantle flow.