Analog and Numerical Experiments of Double Subduction Systems With Opposite Polarity in Adjacent Segments

In this work we study the dynamics of double subduction systems with opposite polarity in adjacent segments. A combined approach of numerical and analog experiments allows us to compare results and exploit the strengths of both methodologies. High‐resolution numerical experiments complement laboratory results by providing quantities difficult to measure in the laboratory such as stress state, flow patterns, and energy dissipation. Results show strong asymmetries in the mantle flow that produce in turn asymmetries in the trench and in the downgoing slab deformation. The mantle flow pattern varies with time; the toroidal cells between the plates evolve until merging into one unique cell when the trenches align. In that moment the maximum upward flow is observed close to the trenches. The interaction between the mantle flow produced by each subducting plate makes the rollback processes slower than in a single subduction case. This is consistent with the observed energy dissipation rate that is smaller in the double subduction system than in two single subductions. Moreover, we provide a detailed analysis on the setup and boundary conditions required to numerically reproduce the analog experiments. Boundary conditions at the bottom of the domain are crucial to reproduce their analog counterparts. Numerical results are compared to natural examples of multi‐slab subduction systems in terms of upper mantle seismic anisotropy, relative trench retreat velocities, and composition of subduction‐related magmatism.


Introduction
The study of subduction zones is of prime importance since they play an essential role as main driving mechanism for plate tectonics and mantle dynamics (Coltice et al., 2019). A large part of the current subduction zones, where a single oceanic plate descends beneath another continental or oceanic plate, are associated with nearly linear trenches that extend over distances of some thousands of kilometers (e.g., Aleutians, Peru-Chile, Kermadec-Tonga, Japan-Kurile, Java-Sunda). Mantle dynamics in such subduction zones is dominated by poloidal cells induced by the entrained flow beneath the slab and the corner flow produced by the down dip of the slab. Toroidal cells around the slab edges can be developed in response to the return mantle flow associated with the trench migration (Gable et al., 1991).
This mantle flow pattern can be largely modified by the presence of vertical and horizontal slab tears and slab fragmentation (e.g., Liu & Stegman, 2012;Long, 2016;Magni et al., 2017). Furthermore, in settings where oceanic plates are highly segmented, subduction is commonly characterized by very arcuate trenches of limited extent showing different slab dip orientations. Indeed, we can observe double subduction systems with parallel trenches with inward-dipping polarity as in Luzon (Bautista et al., 2001), outward-dipping polarity as in Molucca Sea (Zhang et al., 2017), or same-dipping polarity as in Philippine Sea . In addition, there are subduction systems characterized by two adjacent slabs retreating in opposite directions as in Taiwan (Lallemand et al., 2001), New Zealand (Lamb, 2011), the Western Mediterranean (Vergés & Fernàndez, 2012), and the Alps-Apennine junction in Italy (Vignaroli et al., 2008). All these complex subduction systems can modify the sub-lithospheric mantle flow generating seismic anisotropy and the geochemical signature of subduction-related magmatism (e.g., Faccenda & Capitanio, 2013;Ma et al., 2019;Magni, 2019).
The dynamics of double subduction systems have been investigated by a wealth of numerical and analog experiments involving different configurations (e.g., Di Leo et al., 2014;Holt et al., 2017;Mishin et al., 2008;Pusok & Stegman, 2019). In particular, subduction systems characterized by two adjacent slabs retreating in opposite directions have been tested with 3-D numerical (Király et al., 2016) and analog experiments (Peral et al., 2018). Whereas Király et al. (2016) focused on the interactions between the return flows generated by the respective retreating slabs, Peral et al. (2018) investigated plate deformation and variations in trench retreat velocities. Although both models use similar geometries and share main outcomes, differences in the model setup make it difficult to compare them.
Combining computational and laboratory models of the same geodynamic process may help understanding its dynamic evolution by complementing each method's weaknesses and strengths (e.g., Mériaux et al., 2018;Panien et al., 2006). While laboratory experiments provide the physical realism and high temporal and geometrical resolution, numerical models allow for fully controlling and quantifying physical parameters, as velocity and stress, which cannot be directly obtained from laboratory experiments. This said, the numerical reproduction of laboratory results in terms of temporal evolution of the subduction process, trench curvature, and slab geometry helps identifying and understanding boundary conditions, rheological behavior, and other effects in analog experiments.
Only a couple of studies compare analog and numerical experiments of single plate subduction and elaborate differences and similarities obtained by the different methods (Mériaux et al., 2018;Schmeling et al., 2008). Schmeling et al. (2008) pointed out the importance of a zero-density weak top layer to simulate the free surface behavior when trying to reproduce laboratory experiments numerically. More recently, Mériaux et al. (2018) presented 3-D numerical models with the objective of reproducing single plate subduction laboratory experiments. They proposed that surface tension effects of the syrup representing the mantle material in analog experiments may affect the resulting trench retreat velocity and flank stability of subducting plates.
Here, we perform numerical models aiming to reproduce analog experiments of complex double subduction systems as recently published in Peral et al. (2018). In this regard, the geometrical setup and material parameters have been chosen to best represent values applied in the laboratory. The main objective of our work is to better understand the evolution of small-scale subduction systems with opposite polarity in adjacent segments by determining the relevant physical parameters that characterize these processes. The study is divided into three main objectives: (1) Testing the effects of boundary conditions, rheology, and plate thickness on the evolution of single and double plate subduction systems to evaluate how the assumed simplifications in numerical models can affect the results. Furthermore, those results also help to identify the nature of analog boundary conditions and rheological uncertainties. (2) Establishing a reference subduction system with opposite polarity to analyze plate deformation and mantle flow interaction by quantifying velocities, stresses, and forces within such systems. (3) Comparing numerical results to natural systems and interpreting the geological consequences during the evolution of these subduction systems in terms of trench retreat velocities, plate deformation, magmatism, and seismic anisotropy.

Analog Model
Laboratory experiments presented by Peral et al. (2018) were performed in a plexiglass tank with dimensions of 150 × 150 × 50 cm 3 . Materials consist of both linear viscous glucose syrup and silicone putty representing the upper mantle and subducting plates, respectively. Subduction is driven by a density difference between the initially floating plates and the mantle material. All material parameters and respective scaling are described in Table 1. Among the models included in Peral et al. (2018) we have chosen the double subduction configuration consisting of 10 cm wide plates and the single subduction configuration with a 30 cm wide plate for comparison to numerical experiments. The length of the plates (l) measures 30 cm, and the distance between the trailing edges (s) is 39 cm (Figure 1). The plates are fixed at their trailing edge, and the upper/lower mantle boundary is simulated by placing a fixed base at 11 cm depth. Subduction is initiated by manually pushing down the leading edge of the plates into the syrup to about 3 cm depth. For more details about the analog experiments, we refer to Peral et al. (2018).

Numerical Model
The numerical counterparts of single plate and double plate subduction models with opposite polarity have been performed by a three-dimensional code based on the finite difference method with a fully staggered Eulerian grid and freely advecting Lagrangian markers (I3ELVIS; Gerya, 2009;Gerya & Yuen, 2007). Similar to the analog experiments we applied a linear (Newtonian) rheology. The numerical model solves the equation for conservation of mass assuming incompressibility and conservation of momentum (Stokes equation) where v i are the velocities in x-, y-, and z-direction, x i are the spatial coordinates, P is the dynamic pressure, ρ is the density, and g i is the gravitational acceleration. Deviatoric stresses are defined by  (Table 2).
where η denotes the viscosity.
The governing equations are solved with an OpenMP-parallelized multigrid solver running on 16 threads. The duration of single time steps is capped to ensure that no marker moves further than 1/10 of a nodal cell size, with a maximal value of 1 min (average time steps defined by marker movement are around 4 s). Numerical experiments needed around 2 weeks to complete. Respectively, wall time of conducted experiments average out at around 5,400 h (runtime multiplied by amount of threads). The nodal resolution of the models depends on the size of the computational domain (Table 2). All experiments initially contain eight Lagrangian markers per nodal cell. Marker properties are interpolated to nodes by an arithmetic averaging scheme, whereas the velocity field calculated on the Eularian grid is back-interpolated onto the markers by applying the fourth-order Runge-Kutta method.

Model Setup
Geometrical and physical parameters of numerical models have been chosen to best reproduce analog experiments from which we have selected those with narrow plates in the double subduction system and wide plates for the single plate system. All experiments with only one plate exhibit a vertical plane of  (Peral et al., 2018). Due to the relatively large cell size in model N1, plates are initially spaced 2 cm while in the rest of the models the initial lateral distance is of 1 cm, corresponding to the separation that is observed in the laboratory experiments immediately after subduction initiation. All models have a height of 12 cm. The initial distribution of markers is characterized-from bottom to top-by 11 cm of mantle and 1 cm of "sticky-air," which is a low-density, low-viscosity phase imposing low shear stresses along its interface with the mantle/plate and allowing the system to develop a surface topography. The top boundary of the sticky-air layer is no-slip. We have also tested models with free-slip top boundary conditions showing very similar results though consuming a much longer computation time.
The density of the mantle material is 1,445 kg/m 3 , the plates have a density of 1,505 kg/m 3 , and the "stickyair" layer has a density of 1 kg/m 3 ( Table 2). All models but one exhibit a linear viscous rheology for both plates and mantle (Table 2). Plates have an initial thickness of 1, 1.2, or 1.35 cm and are located at the top of the mantle, in contact with the "sticky-air" layer ( Figure 1). Plate lengths are 30 cm with a width of 10 cm for double plate and a width of 10 or 30 cm for single plate models (Table 2). Plates are fixed at their trailing end by predefined null nodal velocity in all directions at a distance of 39 cm to each other parallel to the plate extent ( Figure 1). To initiate density-driven subduction, a small slab perturbation is initially imposed at their tips where the slab penetrates 3 cm into the mantle at 45°dip angle. One additional model has been conducted with a non-linear, strain rate dependent viscosity for the subducting plate (model N15). The non-linear viscosity is calculated by where η 0 denotes the reference viscosity, n is the power law exponent, and ε II the second invariant of the strain rate tensor ( Table 2).
Subduction of the plate(s) in the numerical models is dynamically self-consistent in the sense that it is driven by density contrast only, and no material flux is allowed into and out of the model domain. Lateral and bottom boundary conditions are variably prescribed as no-slip or free-slip (zero shear stresses along boundary) to test their effect on mantle flow and subduction retreat (Table 2).

Comparing Analog-Numerical Experiments
Preliminary numerical models have been run reproducing the laboratory conditions to better understand how the numerical/analog modeling constraints affect the final results in the natural prototype rather than to the laboratory experiments itself. A three-step parametric study has been performed allowing to design the optimum numerical experiment to study the double subduction systems: First, numerical models of double subduction systems with 10 cm wide plates were designed to test the effects of variable model domain size. Second, the importance of applied numerical boundary conditions was investigated. Finally, numerical models with a single 30 cm wide plate were conducted to reveal the effects of plate stiffness by varying plate viscosity and thickness. Geometric and physical parameters of all presented models are listed in Table 2.
Results are described taking into account the different phases of the evolution of a subduction system with opposite polarity in adjacent segments described in Király et al. (2016): (i) initial stage corresponding to the evolution of the system until plates reach the base of the upper mantle (Phase 1); (ii) approaching trenches, starting with the acceleration of slab rollback after the slabs interact with the lower mantle and finishing when trenches intersect, that is, are aligned to each other (Phase 2); and (iii) diverging trenches, spanning from trench intersection until subduction completes (Phase 3). We define trench intersection as the transition from Phase 2 to Phase 3.

Influence of Model Domain Size
Numerical models with different box sizes show similar plate geometries during the different phases of subduction ( Figure 2). Phase 1 is not shown in this section as the trenches are too far from each other to produce any interaction between plates. Phase 2 shows that plates tend to approach each other, this effect being more intense for the medium and small box experiments. In Phase 3, subduction continues, and the slabs show a flattened asymmetric shape lying on the bottom of the model. This asymmetric deformation is less intense in model N1, where plates are initially more separated. The flow pattern at 6 cm depth from the top of the model domain is similar in all models (Figure 2), though the radius of the toroidal flow is smaller as the box size decreases. The maximum mantle velocities (~1.4 mm/min) at this depth are registered in front of the trench and behind the slab during the entire subduction process. Figure 3 illustrates the amount of trench retreat versus time for numerical experiments with different box sizes in comparison to the reference analog experiment. Trenches of the analog experiment retreat with a roughly linear trend (orange and black points in Figure 3 with their linear regression as black line). In

Influence of Applied Boundary Conditions
We first analyze the role of the conditions applied to the lateral boundaries of the model domain. Figure 4 illustrates the temporal evolution of trench retreat of models with free-slip (N1, N2, and N3) and no-slip (N4, N5, and N6) conditions. Applying free-slip conditions allows the mantle material to move along the lateral walls of the box without resistance. Contrariwise, no-slip lateral boundary conditions prevent lateral flow on the surface of the walls. There are no noticeable differences in the resulting trench retreat related

10.1029/2020GC009035
Geochemistry, Geophysics, Geosystems to these two end-member boundary conditions for models with box size of 150 × 150 cm 2 and 80 × 80 cm 2 . However, models N3 and N6 (box size 40 × 30 cm 2 ) show a measurable offset with free-slip boundaries resulting in a slightly faster trench retreat ( Figure 4).
Second, we explore the effects of applying a free-slip or no-slip boundary condition at the bottom of the model domain ( Figure 5). Laboratory experiments show that subducted plates when reaching the basal plate, at least to a certain extent, deform and move horizontally. This observation requires testing different bottom boundary conditions for the numerical counterpart (Peral et al., 2018). A double subduction model with no-slip boundary condition at the bottom implies zero horizontal velocity for the mantle and the slab material at the base of the model domain. Therefore, lateral movement and stretching of the plate lying on the floor of the model is restricted (Figure 5a). On the other hand, free-slip boundary conditions at the bottom allow the slabs to stretch in horizontal directions and slabs can move laterally along the bottom boundary depending on the resulting velocities ( Figure 5b). Trench retreat velocities are up to 15% higher for free-slip than for no-slip conditions and the toroidal mantle flow produces a lateral movement of both slabs that is not observed for the no-slip models.

Influence of Plate Rheology and Plate Thickness
Previous laboratory experiments of single and double subduction have shown that the trenches of subducting plates exhibit a more intense curvature than numerical models when applying the measured spatial and rheological parameters (Figurea 6a and 6b). This effect, which is more pronounced for wider plates, raises the question whether the rheology of the plates and the mantle differ for laboratory and numerical models. Therefore, a series of numerical models with a 30 cm wide single plate was conducted to test the effects of plate rheology and thickness on the subduction process and particularly on the trench curvature (see Table 2 and Figure 6).
Results show that the trench curvature, defined as the ratio between the chord and the sagitta of the circular segment delineated by the trench (Peral et al., 2018), decreases and the deformation and stretching of the slab increases with a decreasing viscosity ratio between the lithosphere and the mantle (γ = η l /η m ). Reducing the viscosity ratio from a reference value γ = 195 to γ = 100 is not sufficient to reproduce the slab deformation observed in the laboratory experiment ( Figure 6c). However, reducing at the same time the viscosity ratio to γ = 100 and the plate thickness to 1.2 cm maintains the trench curvature showing a more pronounced slab deformation (Figure 6d). A similar effect has been observed for a non-Newtonian plate

10.1029/2020GC009035
Geochemistry, Geophysics, Geosystems viscosity ( Figure 6e). The numerical model that best represents the laboratory experiment in terms of plate deformation (slab deformation and trench curvature) is obtained by reducing the plate thickness to 1.0 cm ( Figure 6f). However, the time evolution of this model (model N11) differs notably from that performed in the laboratory.

Reference Model for the Double Plate Numerical Experiment
Reducing the computational domain allows to increase the model resolution without increasing the computational cost and therefore to simulate the laboratory experiment with more spatial accuracy. Regarding the size of the computational domain, negligible differences in terms of plate geometry and trench retreat velocity are obtained when reducing the numerical domain from 150 × 150 cm 2 to 80 × 80 cm 2 and 40 × 30 cm 2 for 10 cm wide plates (Figures 2 and 3). In terms of trench retreat the numerical model that best matches the laboratory experiment is model N3 (40 × 30 cm 2 ), which allows the highest resolution (Table 2). However, the toroidal flow for this model is narrower than for the medium and large models, due to the proximity of the lateral walls to the slabs. Indeed, the effect of applied free-slip or no-slip lateral boundary conditions on trench retreat velocity is only noticeable for the small-box numerical experiments, indicating undesired boundary effects (Figure 4). Furthermore, velocity magnitudes at the lateral walls of free-slip experiments of the small box model N3 are significantly larger than for models N1 and N2 ( Figure S1).  Table 2. "c" indicates the trench curvature defined as the ratio between the chord and the sagitta of the circular segment delineated by the trench (see Peral et al., 2018).

Geochemistry, Geophysics, Geosystems
The bottom boundary condition strongly affects the slab geometry ( Figure 5). The slab deformation and temporal evolution of the subduction system observed in the laboratory experiment are better reproduced by numerical models having no-slip boundary conditions at the bottom of the model domain.
The important effect of rheology and plate thickness on plate deformation is evident from the results obtained from numerical models of single subduction with 30 cm wide plates ( Figure 6). The model with non-linear viscosity plate (model N15; Figure 6e) shows the most acceptable similarity with the laboratory experiment in terms of plate deformation and time evolution, indicating that the materials used in the laboratory may not be perfectly linear viscous. On the other hand, changing the plate thickness or the linear viscosity contrast between plate and mantle strongly affects the temporal evolution of the subduction system ( Figure 6).
Following the above discussed results, to study the evolution of subduction processes with opposite polarity in adjacent segments, we choose a reference numerical setup exhibiting a medium box size (80 × 80 cm 2 ) with boundary conditions of free-slip at the lateral walls and no-slip at the bottom of the model domain (model N2; see Table 2). We chose this model because it shows the best trade-off between numerical resolution, trench curvatures, deformation of plates, and trench retreat velocities. The differences observed in the trench retreat versus time between model N2 and the analog experiment ( Figure 3) can be explained by the larger trench curvature and the slab friction with the bottom boundary exhibit by the analog model, which consumes more energy slowing down the process.

Numerical Model of Subduction With Opposite Polarity in Adjacent Segments
In the following, we present the results obtained from the numerical model N2 consisting of two 10 cm wide plates with an initial separation of 1 cm ( Table 2). The interaction between both plates is investigated by analyzing (1) mantle flow, (2) stress and energy dissipation, (3) plate deformation, and (4) trench retreat velocity.
Results are compared with those from a numerical model of single plate subduction (model N8) and the respective double subduction laboratory experiment (model L1, Table 2).

Mantle Flow
The mantle velocity field induced by the double plate subduction process is calculated at different times and depth levels and entirely presented in Figures S2-S5 in the supporting information. Figure 7 shows the three-component velocity field at 6 cm depth of the model domain corresponding to intermediate mantle depths with the magnitude-less horizontal mantle flow direction in the background. At 65 min, mantle flow exhibits a rotational symmetry of second order with the rotation point at the center of the model (Figure 7; top row). This symmetry pattern is also observed through the whole evolution of the system, which is characterized by four large toroidal cells with flows converging toward the frontside of the trenches and diverging outward from the backside of the trenches. The symmetry axes of the cells are roughly orthogonal and rotate counterclockwise through the different phases as a result of the progressive trench retreating (Figure 7). The orientation of these axes, as well as the size of the cells and its symmetry, depends on the initial geometry of the system (plate width, plate separation, and box size). In our experiment, during Phase 1 (from minute 0 to 53), Phase 2 (from minute 53 to 102), and early Phase 3 (from minute 102 to 140), the induced toroidal mantle flow is asymmetrical with respect to the longitudinal axis (x-direction) of each plate, particularly in the backside of the trenches (see also Figure 9). The two toroidal cells around the adjacent lateral slab edges push the plates toward each other and merge into a single cell during trench intersection. During late Phase 3 (>140 min), the interaction between the adjacent plates vanishes and the toroidal flow cells become nearly symmetrical in the backside of the trenches but strongly asymmetrical in the frontside (Figure 7).
The velocity component v z , parallel to the initial subduction trenches, is the most affected by the interaction between the two plates (

10.1029/2020GC009035
Geochemistry, Geophysics, Geosystems slab that is vanishing along Phase 3. In the inter-plate region, the mantle material flows upward associated with the front side of the slabs ( Figure S5). These upwelling flows approach each other as the subduction progresses changing during trench intersection when upwelling is associated with the backside of the slabs. During intersection, the mantle flow in the inter-plate region changes direction generating a downward vertical flow that is not observed in other regions of the mantle or during other phases (Figures 7c and S2-S5).

Stress and Energy Dissipation
A different way to figure out how plates interact with each other through the generated mantle flow is analyzing the stress distribution and the total energy dissipated during the subduction process. Figure 9 shows the second invariant of the stress tensor σ eff ¼ 1 Á 1=2 calculated at intermediate mantle depth (6 cm from the top of the model domain) and the direction of the horizontal velocity

10.1029/2020GC009035
Geochemistry, Geophysics, Geosystems field. At this depth, stress values of 0.6 Pa are registered in the mantle regions adjacent to the slabs during the entire subduction process. However, as the subduction progresses, the stress increases from 0.05 Pa to more than 0.2 Pa in the vicinity of the internal sides of the slabs reaching maximum values (>0.45 Pa) in the inter-plate region during the trench intersection. In the course of Phase 3, the mechanical coupling vanishes progressively.
Additionally, we have calculated and compared the energy rates dissipated dW ¼V · dt · _ ε ij · τ ij by the mantle of both single and double plate subduction models as shown in Figure 10. In the double subduction model, the maximum rate of energy dissipation occurs during Phase 1 reaching~1.2 · 10 −6 W, before the slabs reach the base of the box when energy dissipation rate decays rapidly to~0.3 · 10 −6 W ( Figure 10; red line). During Phase 2, the energy dissipation rate increases again until 0.9 · 10 −6 W before it decreases to a stable value of 0.7 · 10 −6 W. From minute 150 on to the end of subduction (minute 200) there is a steady increase in energy dissipation rate up to 0.9 · 10 −6 W ( Figure 10; red line). The dissipated energy rate for the single plate model (model N8) shows a similar pattern with values that are slightly larger than half of those calculated for the double plate model (Figure 10; gray line). Indeed, the total dissipated energy calculated for a single plate (0.0044 J) is half of that corresponding to a double plate system (0.009 J) considering 200 min of evolution in both cases. This released energy represents~16% of the total potential energy of the system, the rest being partly held as potential energy as plates are still subjected to their trailing edges, and partly being dissipated by plate deformation.

Plate Deformation
To visualize plate deformation during the double subduction process corresponding to model N2, we illustrate the mantle velocity field at 1.2 cm depth from the top of the model domain, that is, 2 mm into the plates from the sticky-air interface (Figure 11), and the second invariant of the strain rate within the two plates

10.1029/2020GC009035
Geochemistry, Geophysics, Geosystems (Figure 12). The velocity field observed at the plate level indicates that, from the beginning of the experiment until early Phase 3, the mantle flow pushes the plates toward each other in a direction perpendicular to their respective longitudinal axes producing a lateral movement of both plates (Figure 11). This is also observed in Figure 12 in which strain rate increases from the center of the plates to their borders as a result of lateral movement. Trenches show a symmetric curvature that slightly decreases with time accompanied of a counter clockwise rotation due to the lateral movement of plates (Figures 11 and 12). As expected, the highest strain rates occur in the regions where plates bend vertically, that is, near the trenches and at the bottom of the model where plates get horizontal. Figure 13 shows the time evolution of trench retreat velocity for single and double plate subduction systems. In the double subduction model, trench retreat velocities of both plates are roughly equal due to the symmetry of the initial setup (minor differences may occur related to the technique of trench localization which is also dependent on random marker distribution). Phase 1 is characterized by a fast trench retreat until 40 min, followed by a velocity decrease until the tips of the plates reach the base of the model box at around 53 min (Figure 13). During Phase 2, trench retreat velocity increases again reaching a maximum of~1.6 mm/ min followed by a short period of velocity decrease before trenches intersect at around 102 min. Phase 3 is characterized by a progressive increase in retreat velocity reaching maximum values of~1.9 mm/min at the late stage of evolution. The single plate model (model N8) shows similar trench retreat velocity variations during all phases of the subduction process although is slightly faster reaching a maximum velocity of 2.1 mm/min during Phase 1 (Figure 13). The higher speed computed in the single plate model relative to the double plate model is consistent with the calculated energy dissipation rate, which is slightly higher than half of that for the double plate model ( Figure 10). As the total energy is exactly half, the time over which the energy rate is integrated must be lower and, therefore, the velocity faster.

Discussion
In the following, numerical results are discussed with respect to the three main objectives introduced earlier, which include (1) testing the effects of rheology and boundary conditions and compare them to laboratory experiments in terms of temporal and structural evolution, (2) investigating and quantifying mantle flow and plate deformation in a double subduction experiment, and (3) comparing the obtained numerical data to natural cases of double subduction.

Comparison of Numerical and Laboratory Experiments
Numerically reproducing laboratory experiments is not straightforward. Indeed, other authors already found differences in the sinking and retreating rates and plate morphology when attempting to numerically reproduce laboratory experiments of single plate subduction (Mériaux et al., 2018;Schmeling et al., 2008). Performing both models simultaneously helps us to control all the parameters affecting the evolution of the system. The agreement between results from numerical and laboratory models is fairly good when considering an intermediate domain size (80 × 80 cm 2 ) with free-slip boundary conditions for the lateral sides and no-slip at the bottom of the numerical model, although the effect of applying different lateral boundary conditions is weak. First-order observations like temporal evolution of subduction, plate and trench geometries, plate deformation, and mantle flow are similar in both models (Figures 2 and 12). However, the two methodologies produce also slight differences. For example, the experimental uncertainties related to the handling process of the laboratory experiments can lead to asymmetries between plates in the initial stages of the evolution of the system or to slight modifications in the rheology of materials. Nevertheless, we found that these factors, which are complex to quantify, produce only second-order differences.

Geochemistry, Geophysics, Geosystems
Aside these differences related to inaccuracies from the experimental handling, there are others related to simplifications made in the numerical approach. For example, the overall slab deformation as well as the temporal evolution of the subduction process is more accurately reproduced when applying no-slip boundary conditions at the bottom of the model domain ( Figure 5). However, the laboratory experiments show that the plates stretch and extend horizontally once lying on the model box floor (Figure 12), indicating that the analog boundary condition ranges in between those end members (free-slip/no-slip). An intermediate boundary condition is also expected for the horizon separating the upper from the lower mantle, where an increase in viscosity with a factor of 10 to 100 indicates resistance but does not prohibit viscous drag (e.g., Goes et al., 2017;Király et al., 2017).
Furthermore, numerical experiments cannot reproduce simultaneously the large trench curvature and the retreat velocity observed in laboratory experiments, particularly for wide plates, even when varying the rheological parameters (Figure 6). This points toward potential effects that were not taken into account as the formation of a thin "crystalized" layer at the surface of the syrup, which might generate surface forces not considered in the numerical model (e.g., Mériaux et al., 2018).

Geochemistry, Geophysics, Geosystems
Concerning double subduction systems with opposite polarity, numerical and analog experiments are in general agreement in terms of overall duration of the subduction process. Furthermore, stresses generated in the inter-plate region in both numerical and analog models, tend to separate the plates from each other (Figures 9 and 11; and Figure 7 in Peral et al., 2018). However, some discrepancies between the results of the two modeling approaches should be discussed in more detail: Analog experiments show a divergent displacement of the two plates, whereas in numerical experiments the resulting deformation of plates shows a lateral movement that brings the central parts of the plates closer together (Figures 11 and 12). This discrepancy can be explained by the different radii of the toroidal cells generated in both experiments, which seems to be larger in the numerical than in the analog models independently on the size of the modeling domain ( Figure 2). Another discrepancy is observed related to trench retreat velocities. Laboratory experiments show retreat velocities increasing until plate intersection and decreasing thereafter (Peral et al., 2018), whereas velocities remain roughly constant or even increase toward the end of the subduction process for numerical experiments (Figures 3 and 13). A possible explanation is the different mantle flow patterns and associated evolution of energy dissipation (Figure 10) that are affected by implemented boundary conditions in numerical models and experimental uncertainties in the analog setup, as well as different rheological properties between analog and numerical materials. Concerning the shapes of subducting plates, asymmetric trench curvatures are observed in both numerical and laboratory experiments suggesting that these asymmetries are related to the mantle flow interactions associated with double subduction systems ( Figure 12). A more detailed study of the mantle flow in analog experiments is needed to better understand the observed differences.

Dynamics of Double Subduction Systems
Mantle flow in an ideal single slab subduction system is characterized by a toroidal component that is symmetric relative to the longitudinal axis of the slab. In the case of two adjacent slabs with opposite polarity, the toroidal flows induced by both slabs interact with each other generating different flow patterns in the inter-plate and outer mantle regions. Consequently, the resulting toroidal component of the mantle flow becomes asymmetric relative to the slabs axes (Figures 7 and 11). This asymmetry prevents for applying a half-space model domain as for example in Király et al. (2016), whose setup implies an indefinite repetition of slabs laterally by introducing free-slip side boundaries cutting the opposite polarity slabs. Such boundary conditions are strictly valid for plates exceeding 1,000 km width or for tectonic scenarios where the opposed

10.1029/2020GC009035
Geochemistry, Geophysics, Geosystems polarity repeats several times. In any case, the resulting flow pattern has important implications for subduction driven by rollback, where the toroidal flow makes up 95-100% of the entire mantle flux (Schellart et al., 2007). For example, assuming that seismic anisotropy in the upper mantle is related to the lattice preferred orientation of olivine as a result of large-scale mantle flow (e.g., Long & Becker, 2010;Silver, 1996) allows detecting complex subduction scenarios and understanding related mantle dynamics (Alpert et al., 2013;Wei et al., 2016). For double subduction systems similar to the ones presented in this study, the orientation of the symmetry axes of the toroidal cells may help interpreting inter-plate coupling and the spatial framework of adjacent trenches (Figure 7a).
Deformation of plates is coupled to the mantle flow and, therefore, affected by mantle flow asymmetry producing lateral movement of the plates and additional deformation of the slabs in double subduction systems ( Figure 12). The deformation experienced by the plates is related to the mantle drag exerted by the net outward flow separating one plate from the other. The strength of this drag decreases with the distance between the slabs, so its effect is maximum during the intersection of the two trenches. Similar behavior has been reported by Király et al. (2016). The interaction between the flows generated by the two slabs reduces the energy dissipation rate compared to two isolated plates. As a consequence, the trench retreat velocities in a double subduction system are slowed down with respect to a single subduction process. An exception is the end of Phase 3 (t >150 min), where plates are sufficiently far from one another and trench retreat is accelerating rather than keeping a constant velocity as in the case of a single plate subduction (Figures 10 and 13).
Despite the strong interaction between plates in a double subduction system, the evolution of the trench velocity through time shares some common features with the single subduction system. Obviously, these similarities are stronger during the initial and final stages of the evolution, when the two slabs are distant and the interaction of their induced mantle flows is weaker. For example, as reported elsewhere (Funiciello et al., 2003;Funiciello et al., 2006;Schellart, 2004;Strak & Schellart, 2014), trench velocity increases rapidly at the beginning of the model while the slab is sinking and the negative buoyancy increases (Figure 13; 0 to~40 min). When the slab reaches the bottom of the domain (40 to 55 min), the trench velocity slows down and accelerates again (from 55 to 75 min) until the steady-state subduction is achieved (75 to 180 min). During the steady-state period, the single subduction system shows some minor periodic accelerations and decelerations probably due to changes in the slab angle ( Figure 13). These changes are not so evident in the double subduction system as the interaction between slabs and mantle flow is probably overprinting the slab dip changes.

Relevance for Natural Prototypes
Complex subduction systems formed by several slabs, with reduced dimensions and with advancing or retreating trench migration velocities, can produce changes in the stress, strain rate, and pressure and temperature conditions in the mantle modifying the lattice preferred orientation of olivine crystals and seismic anisotropy as manifested by shear-wave splitting (Faccenda, 2014;Faccenda & Capitanio, 2013). Similarly, the incorporation of sediments and hydrous fluids into the sub-lithospheric mantle during subduction and the rise of the asthenosphere through slab tears or in the back-arc extensional basins can result in volcanic activity with a variety of geochemical signatures (e.g., Lustrino et al., 2011;Melchiorre et al., 2017).
As mentioned in section 1, natural prototypes where double subduction systems with opposite polarity in adjacent segments have been proposed to occur have been identified in Taiwan, New Zealand, the Western Mediterranean, and the Alps-Apennine junction, among others. However, extracting conclusions applicable to these natural scenarios is not straightforward as the presented study is based on oversimplified models of double subduction systems. Main limitations are the lack of overriding plates, the assumption of a viscous rheology for the lithosphere and upper mantle, and the dimensions of plates, which at nature scale are 600 km wide and 1,800 km long separated by a 60-120 km wide transform zone. Despite these limitations, we can infer some expected effects derived from the interaction between the plates and the upper mantle in terms of seismic anisotropy and magmatism.
One of the regions where these effects are supported by observations is the Western Mediterranean. There, the Alboran-Tethys and Algerian-Tethys segments retreated in opposite directions inferred from the present distribution and vergence of the metamorphic complexes in the Betic-Rif and Tell-Kabylies chains (e.g., Casciello et al., 2015;Fernàndez et al., 2019;Vergés & Fernàndez, 2012). It must be noted that the 10.1029/2020GC009035 Geochemistry, Geophysics, Geosystems dimensions of the Ligurian-Tethys domain, whose closure gave rise to the present Western Mediterranean, were narrower than the slabs modeled in this study, and the trench intersection occurred before the slabs reached the lower mantle (i.e., during Phase 1). Fast polarity directions (FPD) inferred from SKS shear wave splitting onshore show a consistent anisotropy pattern oriented parallel to the Betic-Rif orogen (e.g., Díaz et al., 2015;Miller et al., 2013). Despite the poor coverage of anisotropy data offshore, a regional 3-D azimuthally anisotropic model of Europe shows a WSW-ENE alignment of FPD in the Alboran Basin changing to NW-SE in the Algerian Basin at depths of 70-200 km (Zhu & Tromp, 2013). At shallower depths however, Pn and Sn tomography shows an almost perpendicular anisotropy pattern for the uppermost mantle with NW-SE oriented FPD in the southern margin of the Alboran Basin varying to NE-SW in the Algerian Basin (Díaz et al., 2013). These variations in the anisotropy pattern are related to the interaction between the mantle return flows produced by the slabs retreating in opposite directions, the northwest to west displacement of the Alboran-Tethys slab and its tightening, and the absolute motion of the African plate (e.g., Spakman et al., 2018).
The vertical component of mantle flow shows maximum velocities in the back-arc and inter-plate regions being more active during Phases 1 and 2 (Figures 7 and S2-S5). This mantle upwelling from the deeper parts of the upper mantle can produce widespread melting by adiabatic decompression with variable composition depending on the water content and the depletion degree of the mantle source, which changes spatially with the retreating of slabs as proposed from recent 3-D numerical modeling (e.g., Magni, 2019). This variable volcanic geochemical signature is observed in the Western Mediterranean magmatic manifestations changing from tholeiitic in the early Oligocene Malaga dikes to calc-alkaline/HK calc-alkaline and alkaline magmas during Late Miocene to Quaternary. The combination of subduction and rifting events operating in the region produced recycling and depletion of mantle rocks making difficult to associate the derived volcanic products with the present-day geodynamic setting (e.g., Carminati et al., 2012;Lustrino et al., 2011;Lustrino & Wilson, 2007;Melchiorre et al., 2017). The limited observations of anisotropy in the offshore regions of the Western Mediterranean, together with the large variety of magmatic compositions, make difficult to favor a unique geodynamic interpretation among those proposed for this region (e.g., Duggen et al., 2005;Faccenna et al., 2004;Platt et al., 2013;Spakman & Wortel, 2004;Van Hinsbergen et al., 2014;Vergés & Fernàndez, 2012).
Some of the regions where opposed subduction polarity in adjacent plate segments is observed show a strong trench asymmetry with a tight curvature in one of the ends of the trench (e.g., western Alps, western Betic-Rif, southwest Ryukyu trench). Király et al. (2016) proposed that the stress propagation through the mantle produced by the adjacent retreating slabs contributed to the strong curvature of the Western Alps and the SW-Ryukyu trench. However, according to our results, reproducing the arc tightening observed in the above-mentioned regions would require additional "tectonic" conditions to slow-down or even preventing the slab retreat in one of the plate edges to force the trench curvature.
The case of New Zealand where the Pacific plate dips to the west beneath the North Island and the Australian plate dips to the east beneath the South Island is more elusive. The interaction between both slabs is doubtful as they are presently separated by the 500-600 km long right-lateral transform fault, which is the critical distance for stress propagation. The curvature of the Hikurangi Trough could be caused by the fast clockwise rotation of the North Island and the collision of the Chatham Rise microplate with the northern part of the South Island (Wallace et al., 2004;Wallace et al., 2009). In addition, the comparison with the presented experiments is hindered because subduction of both plates is not synchronous, initiating 25 Ma in the North Island and 10-12 Ma in the South Island (King, 2000;Schellart et al., 2006). Summarizing, the present work allows for identifying some distinctive features related to subduction systems with opposed polarity in adjacent segments. However, its application to natural prototypes requires a detailed analysis of the tectonic evolution of the study region and the incorporation of a more sophisticated model setup.

Concluding Remarks
This study presents 3-D numerical experiments of subduction systems with opposite polarity of adjacent segments with geometrical and rheological parameters taken from published analog experiments. First, the effects of input parameters such as Eulerian domain size, boundary conditions, and plate geometry and 10.1029/2020GC009035 Geochemistry, Geophysics, Geosystems strength were tested and compared to analog results. Second, numerical data were used to better quantify physical processes characterizing the evolution of such double subduction systems and finally for comparison to natural examples.
Comparing numerical to analog experiments allows to conclude the following: 1. The Eulerian domain size has a second-order effect on the overall geometry of plates and the trench retreat velocities during the subduction process. However, elevated velocities along the boundaries of the smallest box setup (30 × 40 cm 2 ) indicate that there is a critical small domain size beyond which mantle flow is affected. 2. The choice of lateral boundary conditions (free-vs. no-slip) has a minor effect on the subduction process for the intermediate size experiments. For the bottom boundary, a no-slip condition results in plate geometries better comparable to analog experiments, while a free-slip condition moves the bottom-laying plates away from each other. 3. The trench retreat evolution of subduction systems (single and double) is reproduced very well with rheological and geometrical input data retrieved from analog experiments. However, the trench curvature values obtained from numerical experiments are smaller than the corresponding analogs. This might be due to uncertainties in the measured viscosities from analog experiments or spurious deformation of the plates during handling. 4. Numerical experiments with thinner and weaker plates can reproduce the trench curvature of analog experiments but fail in reproducing the trench retreat velocity. This might be due to the crystallization of a fine film on the surface of the syrup.
The physical characterization of double subduction systems and its comparison to natural examples gives the following insights: 1. The mantle flow induced by the subduction of adjacent plates pushes the plates against each other. Simultaneously both downgoing slabs and trenches deform asymmetrically, and the subduction process evolves slower than the single plate model. During trenches intersection, maximum stresses contribute to further plate deformation. 2. In the horizontal plane mantle flow forms four toroidal cells with symmetry axes rotating during trench retreat. The flow lines converge toward the frontside of the trenches and diverge from the backside of the trenches. 3. The upward component of mantle flow is maximum around the sinking slabs before the intersection of trenches and decreases as the slabs retreat concentrating around the edges of the slabs after trenches intersection. 4. The energy dissipation rate of a plate in a double subduction system is smaller than that of a single subduction. Accordingly, double plate systems exhibit slower trench retreat velocities and a longer duration of the subduction process than a single plate. 5. The subduction models with opposite polarity in adjacent segments predict complex patterns of seismic anisotropy, magmatic composition, and plate deformation associated with the interaction between the sinking slabs and the surrounding mantle. However, its application to natural scenarios needs of a more sophisticated model setup in terms of plates geometry, rheology, and incorporation of overriding plates.