Effect of Plate Length on Subduction Kinematics and Slab Geometry: Insights From Buoyancy‐Driven Analog Subduction Models

Subduction style is controlled by a variety of physical parameters. Here we investigate the effect of subducting plate length on subduction style using laboratory experiments of time‐evolving buoyancy‐driven subduction in 3‐D space. The investigation includes two experimental sets, one with a lower (~740) and one with a higher (~1,680) subducting plate‐to‐mantle viscosity ratio (ηSP/ηM). Each set involves five models with a free‐trailing‐edge subducting plate and variable plate length (20–60 cm, scaling to 1,600–4,800 km), and one model with a fixed‐trailing‐edge subducting plate representing an infinitely long plate. Through determining the contact area between subducting plate and underlying mantle, plate length affects the resistance to trenchward motion of the subducting plate and thus controls the partitioning of the subduction velocity (vS) into the subducting plate velocity (vSP) and trench velocity (vT). This subduction partitioning thereby determines the subduction style by controlling the dip angle of the slab tip once it first touches the 660 km discontinuity. The low ηSP/ηM models display two subduction styles. Short plates (≤40 cm) induce a higher subduction partitioning ratio (vSP/vS), promoting trench advance and slab rollover geometries, whereas longer plates (≥50 cm) lead to a lower vSP/vS, producing continuous trench retreat and backward slab draping geometries. In contrast, the high ηSP/ηM models exclusively show trench retreat with draping geometries, as the high ηSP/ηM enables less slab bending before its tip touches the 660 km discontinuity. Our study indicates that future modeling work should consider the effects of plate length on the style and evolution of subduction.


Introduction
Natural subduction systems come in many forms and can be grouped according to their specific subduction kinematics and slab geometries. Subduction kinematic studies focusing on trench motion commonly distinguish between "trench retreat" and "trench advance" modes in nature (Heuret & Lallemand, 2005;Jarrard, 1986;Kincaid & Olson, 1987;Schellart et al., 2008). Furthermore, a variety of slab geometries can be observed from seismic tomography images (e.g., Amaru, 2007;Chang et al., 2015;Goes et al., 2017;Li et al., 2008;Simmons et al., 2012;van der Meer et al., 2018;Van der Voo et al., 1999;Zhao, 2001). The two trench motion modes combined with the slab geometries describe the three main subduction styles of natural subduction systems: (1) trench retreat and slab rollback with a backward slab draping geometry at the 660-km discontinuity (lazy "S" or "Z" shape, e.g., Izu-Bonin subduction zone; Figure 1a), (2) trench advance and slab roll-forward with a rollover geometry (a U shape on its side) in the upper mantle (e.g., India-Eurasia collision zone; Figure 1b), and (3) alternating trench retreat and trench advance and a steep slab with a slab folding geometry at and below the 660-km discontinuity (e.g., Mariana subduction zone; Figure 1c). These different subduction styles depend on the plate and plate boundary velocities, the internal properties of the subducting slab, the ambient mantle rheology, and the interaction between the plates. Both the properties of the subduction zone itself as well as external parameters are important controlling factors for the variable subduction styles. Subduction zone width Schellart et al., 2007;Stegman et al., 2010;, slab-mantle viscosity ratio Funiciello et al., 2008;Li & Ribe, 2012;Ribe, 2010;Schellart, 2008), different trailing edge boundary conditions on the subducting plate (Funiciello et al., 2004;Schellart, 2005) and on the overriding plate (Čížková & Bina, 2013;Heuret et al., 2007), thickness of the subducting plate (Conrad & Hager, 1999;Funiciello et al., 2008;Schellart, 2004a), mantle viscosity stratification (Christensen, 1996;Funiciello et al., 2003;Kincaid & Olson, 1987), mantle density stratification (Christensen, 1996), and the mechanical coupling at the subduction zone interface (Čížková & Bina, 2013;Duarte et al., 2013) are a few examples of physical parameters that have been investigated in studies that aim to understand the formation of the different subduction styles.
Most of the studies that analyzed the mechanism behind the controlling factors of subduction styles use analog and/or numerical modeling techniques (e.g., Bellahsen et al., 2005;Billen, 2008;Boutelier & Cruden, 2008;Capitanio et al., 2007;Carluccio et al., 2019;Cerpa et al., 2018;Christensen, 1996;Čížková & Bina, 2013;Crameri & Lithgow-Bertelloni, 2018;Di Giuseppe et al., 2008;Duarte et al., 2013;Funiciello et al., 2003;Garel et al., 2014;Gerardi et al., 2019;Griffiths et al., 1995;Guillaume et al., 2009;Holt et al., 2015;Kincaid & Olson, 1987;Li & Ribe, 2012;Ribe, 2010;Schellart, 2004a;Schellart et al., 2007;Stegman et al., 2006;. These modeling studies often investigate more than one parameter and demonstrate that the controlling factors sometimes work together to produce a specific subduction style. For example, Schellart (2004a) found that an increase of slab thickness combined with an increased slab-mantle density contrast promotes trench retreat by increasing the negative buoyancy of the slab. In addition, an increased slab thickness also increases the bending radius to keep the bending resistance low, which promotes trench retreat even further (Capitanio et al., 2009;Irvine & Schellart, 2012;Schellart, 2004a). Another parameter that influences the kinematic characteristics of a subduction zone is the boundary condition at the subducting plate trailing edge. A fixed trailing edge exerts a resisting force on the trenchward motion of the subducting plate that produces a backward sinking slab and only results in trench retreat, while a free trailing edge or a pushed trailing edge promotes a more vertical or forward sinking slab (Funiciello et al., 2004;Schellart, 2004aSchellart, , 2005 and can produce all three main different subduction styles, depending on other physical parameters (Bellahsen et al., 2005;Funiciello et al., 2008;Schellart, 2004aSchellart, , 2008. Yet another example is the slab-mantle viscosity ratio (η SP /η M ), which can be inferred from models. With the use of a variety of methods, different viscosity ratio ranges (~100-500) Figure 1. Illustration of three types of slab geometry (first row) interpreted from cross sections of seismic tomographic images (second row) from (a, c) two subduction zones and (b) from one continental subduction zone, with red lines showing location in the map (third row). Seismic images are from Amaru (2007) and were obtained from the SubMachine website (Hosseini et al., 2018). (a) Backward slab draping geometry, (b) rollover geometry, and (c) steep slab with slab folding geometry at and below the 660-km discontinuity.
were identified empirically for nature Ribe, 2010;Schellart, 2008;Wu et al., 2008). The different η SP /η M produce different slab geometries by influencing for example the relative viscous drag around the slab and by controlling the bending angle of the subducting slab (Bellahsen et al., 2005;Funiciello et al., 2008;Schellart, 2008). Schellart (2008) and Ribe (2010) proposed that different values of η SP /η M generate different degrees of slab bending before the slab tip touches the upper-lower mantle boundary, thereby determining the subduction partitioning and slab geometry during subsequent subduction. In addition, several models also show that the initial slab dip angle is another important factor that determines the subduction kinematics and geometry, by controlling the slab dip angle when the slab arrives at the upper-lower mantle boundary (Crameri & Lithgow-Bertelloni, 2018). A lower slab dip angle at the upper-lower mantle boundary promotes the slab to drape backward and the slab to form a rollback geometry, whereas a nearly vertical slab dip favors slab roll forward with a rollover geometry or vertical penetration into the lower mantle.
In terms of energy loss, the potential energy of the negatively buoyant subducting slab is dissipated into four different mediums: (1) the sublithospheric mantle to drive mantle flow, (2) the slab during bending, (3) the subduction zone interface by shearing, and (4) the overriding plate during deformation (Capitanio et al., 2007(Capitanio et al., , 2009Chen et al., 2015;Conrad & Hager, 1999;Duarte et al., 2013;Irvine & Schellart, 2012;Krien & Fleitout, 2008). Some early studies supported the idea that slab bending consumes most of the potential energy (reaching 60-95%) (Bellahsen et al., 2005;Conrad & Hager, 1999;Funiciello et al., 2003). However, more recent studies argue that only a minor part of the potential energy (<20%) is needed for slab bending and a much higher percentage (70-97%) is required to drive the flow in the mantle (Capitanio et al., 2007(Capitanio et al., , 2009Gerardi et al., 2019;Irvine & Schellart, 2012;Krien & Fleitout, 2008;Leng & Zhong, 2010;Schellart, 2009Schellart, , 2010Stegman et al., 2006). Although a large component of energy is dissipated in driving mantle flow around the slab and bending the slab, a certain amount of energy must be present to drag the subducting plate toward the trench. Stegman et al. (2006) calculated that of the total 82% potential energy that dissipates through mantle flow approximately 69% contributes to the toroidal mantle flow and 13% contributes to the poloidal flow, and part of this is dissipated in the viscous drag experienced at the base of the subducting plate. Energy dissipation through slab bending is controlled by slab viscosity, thickness, bending radius and subduction velocity, whereas plate size plays an important role regarding energy dissipation in the mantle below the base of the subducting plate. Plate motion is responsible for part of the~13% of the energy dissipated in poloidal mantle flow generated by the subduction system due to the viscous shear drag that occurs at the base of the plate. However, the effect of plate width has been studied extensively (Bellahsen et al., 2005;Schellart, 2004a;Schellart et al., 2007;Stegman et al., 2006Stegman et al., , 2010, while plate length has seldom been systematically investigated, especially using analog modeling methods. Plate length has been studied in combination with plate velocity using observational data (Conrad & Hager, 1999) and the authors found that there is no correlation between plate velocity and plate length. Plate length has also been investigated by some numerical modeling studies (Holt & Becker, 2017;Li & Ribe, 2012;Ribe, 2010). Li and Ribe (2012) and Ribe (2010) found that plate length plays only a minor role on the slab geometry, while Holt and Becker (2017) found that the subducting plate length affects the subduction kinematics and the slab geometry. However, numerical models from Li and Ribe (2012) and Ribe (2010) investigated a relatively short range of plate lengths and Holt and Becker (2017) used a relatively small number of increments of different plate lengths. Moreover, many analog subduction experiments have been performed, but with different plate lengths used in different studies (Table 1), which prevents a direct comparison between the models. Thus, it remains unclear how exactly subducting plate length affects subduction evolution and the style of subduction. Since previous analog studies have not investigated plate length and its effects on the subduction evolution, and previous numerical studies have not investigated a large range of plate lengths, we conducted this study to investigate how plate length influences the kinematic and geometric style of subduction, as well as the viscous drag force at the base of the plate. The experiments were conducted in a fully dynamic environment, with negative buoyancy as the only driving force for subduction. We varied three parameters: (1) subducting plate length, (2) the subducting plate to upper mantle viscosity ratio η SP /η M , and (3) the trailing edge boundary condition. With this study we aim to quantify the effect of plate length on the style of subduction and to increase our understanding of subduction evolution and dynamics.

Model Setup, Materials, Boundary Conditions, and Assumptions
We present 12 laboratory experiments to model subduction in an upper mantle reservoir. Four experiments have been conducted twice to ensure reproducibility. Experimental materials similar to previous modeling studies were used (Chen et al., 2016;Funiciello et al., 2008;Meyer & Schellart, 2013;Schellart, 2004aSchellart, , 2004b. The model box consisted of a rectangular plexiglass tank of 80 × 60 × 8.25 cm 3 ( Figure 2). The tank was filled with glucose syrup representing the sublithospheric upper mantle. A highly viscous sheet of mixed silicone putty and fine iron powder was placed on the glucose syrup to represent the subducting lithospheric plate. The viscosity of glucose syrup is strain rate independent but depends on temperature (Schellart, 2011). The syrup used in this study has a low viscosity of~38 Pa s and a density of~1,408 kg/m 3 at 21°C. At the start of each experiment, we measured the temperature and viscosity of the glucose syrup. The viscosity of the glucose syrup was measured with a sinking sphere using Stokes' law. The subducting silicone plate has a density of 1,508 kg/m 3 . The density contrast between the subducting plate and the sublithospheric mantle is 100 kg/m 3 . Such a density contrast represents a subducting oceanic plate of more than 80 km thick in nature assuming pervasive crustal eclogitization (Cloos, 1993). We used silicone putties with two different viscosities to obtain two different η SP /η M . The rheological properties of the subducting oceanic plate material were tested with a rheometer (Anton Paar MCR 702). The silicone putty mix with the lower η SP /η M has a viscosity of~2.8 × 10 4 Pa s (η SP /η M =~740), while the mix with the higher η SP /η M has a viscosity of 6.4 × 10 4 Pa s (η SP /η M =~1,680). In all experiments, a plate thickness of 1.2 cm (representing 96 km in nature) and a width of 10 cm (representing 800 km in nature) were used. We varied the plate length in 10 cm increments between 20 and 60 cm for two sets of experiments with different η SP /η M (Table 2). Each set of experiments includes five models with a free trailing edge varying only the plate length  Bellahsen et al. (2005) Analog modeling 70~4,080 Funiciello et al. (2008) Analog modeling 40 2,400 Schellart (2008) Analog modeling 55 3,300 Ribe (2010) 2-D numerical 2000 and 3,600 Li and Ribe (2012) 3-D numerical 1,150-2,300 Holt and Becker (2017) 2-D and 3-D numerical 2000, 3,000 and 4,000 This study Analog modeling 20-60 and Fixed trailing edge 1,600-4,800 and infinite  (Schellart, 2004a), the free trailing edge of the subducting plate represents the mid-ocean ridge and the fixed trailing edge may represent a slab attached to a very large, slowly moving tectonic plate. The plate with the fixed trailing edge can also represent an extreme boundary condition, representing an infinitely long subducting plate.
The top of the model is a free surface and the sidewalls and the bottom of the box are rigid (no-slip) boundaries, with the bottom of the tank representing an impenetrable upper-lower mantle boundary. The subducting plate was placed laterally in the middle of the tank before subduction initiation in all experiments, leaving a space of 25 cm from each side of the plate to the front wall and back wall of the tank. In the free trailing edge models, a distance of 3 cm was kept between the trailing edge of the plate and its nearest side wall before initiation such that in all experiments the subduction zone at the leading edge of the plate was located sufficiently far from all the lateral sidewalls (~20 cm in Exps. 5 and 12 to~60 cm in Exps. 1 and 7) during an experiment. The space between the plate and the sidewalls of the tank reduced the boundary effects on the experimental subduction system.
This study does not consider the effects of other physical parameters within a subduction zone, for example, subduction zone width, plate thickness, and the friction at the subduction interface. In addition, some other simplifications exist. The overriding plate is not modeled in this study, and hence, the subduction interface is as weak as the ambient mantle material. Thus, more energy is available to deform the mantle and subducting plate instead of being used to deform the overriding plate (Chen et al., 2015;Holt et al., 2015;Sharples et al., 2014) or to dissipate at the subduction interface (Duarte et al., 2013). This simplification may lead to an overestimation of the forces deforming the subducting plate. However, previous studies imply that the subduction zone interface is very weak (e.g., Duarte et al., 2013Duarte et al., , 2015Seno, 2009;Wang et al., 1995), indicating that the frictional resistance at the subduction interface is much smaller than the slab pull force. In addition, an earlier study has also shown that only a small portion of the potential energy of the subducting slab is used to deform the overriding plate (Chen et al., 2015).
In our experimental work we do not consider any vertical or lateral variation in temperature and petrological properties. In nature, the increase of temperature and pressure with depth triggers metamorphism that increases the negative buoyancy of the subducting slab at that depth (Cloos, 1993). However, metamorphic reactions cannot be modeled in our experiments and a constant slab density has to be applied here. The increased negative buoyancy is compensated by a higher initial slab-mantle density contrast used in this study. Even though lateral and vertical variation of temperature may affect the mantle flow pattern, we ignore this aspect since implementing realistic temperature gradients in our buoyancy-driven laboratory models is beyond the scope of the current work.
In nature, subduction styles are affected by the upper-lower mantle boundary, with a viscosity increase by a factor of 30-100 (Hager, 1984;Kaufmann & Lambeck, 2000;Steinberger & Calderwood, 2006). Our study applies a rigid boundary as the upper-lower mantle discontinuity at 660-km depth, following earlier subduction modeling studies (e.g., Bellahsen et al., 2005;Capitanio et al., 2007;Chen et al., 2016;Duarte et al., 2013;Funiciello et al., 2003Funiciello et al., , 2004Schellart, 2004aSchellart, , 2005. Such a rigid bottom boundary could affect the subduction styles differently compared to nature with a relatively soft and penetrable upper-lower mantle boundary. Funiciello et al. (2003) and Schellart (2008) discussed the differences in subduction kinematics in models with only an upper mantle reservoir involved and in models with a deep mantle. Schellart's (2008) study shows that subduction in the deep mantle is characterized exclusively by trench retreat, whereas subduction within the upper mantle shows both trench retreat and trench advance, depending on different slabto-mantle viscosity ratios. Schellart and Moresi (2013) have shown that the implementation of a rigid 660-km discontinuity rather than a more realistic viscosity step of~100 does not make significant difference on the experimental results for models with a narrow slab (width = 800 km). In the current study we present subduction models with an equally narrow slab (scaling to 800 km), and thus, the simplification of a rigid bottom boundary is justified.
In the lab experiments with silicone putty overlying syrup, surface tension exists as another resistance, which forms at the contact boundaries between the plate and the upper mantle. Its existence can be observed from the troughs along the subduction zone plate boundary (Schellart, 2008), although it is not easily quantified. The surface tension was estimated to be very small (Jacoby, 1976) and is generally compensated by increasing the slab/mantle density contrast in the analog experiments (Duarte et al., 2013;Meyer & Schellart, 2013).

Approach and Scaling
Our experimental approach is fully dynamic and subduction is accomplished by the force balance between the negative buoyancy force (gravitational force) of the subducting slab and the resisting forces within the experimental system (Kincaid & Olson, 1987;Meyer & Schellart, 2013). This is referred to as an "internal approach" because all forces driving and resisting motion are internal to the experimental system (no external velocities or forces are applied) .
A length scale ratio of l m /l n = 1.25 × 10 −7 (superscript m represents model, superscript n represents natural prototype) was adopted, such that 1 cm in the experiment represents 80 km in nature. This means that our shortest tectonic plate, with the absolute dimensions of 10 × 20 × 1.2 cm 3 , represents 800 × 1,600 × 96 km 3 , and the longest plate of 60 cm represents 4,800 km in nature. To keep dynamic similarity between our experiments and nature, the experiments should have a similar flow regime as natural subduction systems, which is determined by estimating the Reynolds number (Re). Re needs to be significantly smaller (more than an order of magnitude) than 1 (Schellart, 2008). Re is defined as such where l represents a characteristic length scale of the object in the fluid reservoir in the experiments such as the length, width or thickness of the slab, v is the characteristic velocity in the experiment such as the sinking velocity of the slab or the slab rollback velocity (v =~0.05 mm/s), and ρ and η are the density and viscosity of the ambient upper mantle material. In this work Re is calculated to be~3.7 × 10 −4 ≪ 1, which means that our subduction experiments are conducted in a symmetrical laminar flow regime (with upstream-downstream symmetry), where inertial forces are negligible and viscous forces dominate.
The models were scaled using the Stokes' settling law equation (Duarte et al., 2013;Jacoby, 1973;Schmeling et al., 2008;: where v represents the velocity of the subducting slab, Δρ is the density contrast between the subducting slab and the ambient upper mantle, l is a characteristic length, g is the gravitational acceleration, η is the viscosity of the upper mantle, and C is a constant. To achieve kinematic and dynamic similarity between our models and the natural equivalent, the force balance between the driving body forces and the resisting viscous forces has to be identical in both our models and the natural case (Duarte et al., 2013;Jacoby, 1973;Schmeling et al., 2008;: Since v in Equation 3 can be replaced by l/t, we can express the scaling ratio of viscosity between model and natural prototype as follows: We used a time scale ratio of t m /t n = 3.8 × 10 −12 (1 s represents 8,300 years in nature), which follows from Equation 4 and our choice of viscous analog materials to represent the upper mantle in nature. Therefore, the plate velocity (v = l/t) of 0.01 mm/s in the experiment represents 0.96 cm/year in nature. The density contrast ratio is Δρ m /Δρ n = 100 kg/m 3 /80 kg/m 3 = 1.25. The experiments are conducted in the Earth's gravity field, and hence for the gravitational acceleration g m = g n . Given the above scaling factors, we obtain a scaling ratio for viscosity η m /η n = 5.94 × 10 −19 . The glucose syrup used in the experiments has a viscosity of 38 Pa s that thus corresponds to a scaled viscosity of~6.4 × 10 19 Pa s for the sublithospheric upper mantle in nature, which is within the estimated range of 10 19 -10 21 Pa s (Artyushkov, 1983;Ranalli, 1995). The lithospheric plates used in the experiments have average viscosities that scale to natural values of~4.7 × 10 22 and 1.1 × 10 23 Pa s for the low and high η SP /η M , respectively.

Initiation and Recording of Experiments
The experiments were initiated by bending the leading edge of the subducting plate downward by pouring a controlled volume (~8 ml) of glucose syrup on top of the first~3 cm of the leading tip. This approach guaranteed to obtain a controlled and comparable initial slab perturbation with a slab length of~3 cm and añ 35-40°dip angle. By initiating subduction this way, we drastically reduced the (unwanted) variability in slab dip angle and slab length of the initial slab perturbation between the models.
Two cameras were used to record the subduction process from the top and from one side with a 20 s time increment. We used passive white tracers placed on top of the subducting plate to investigate the subducting plate velocity (v SP ), trench velocity (v T ), and subduction velocity (v S = v SP + v T ). The tracers were placed along the midline of the top surface as well as on both sides of the subducting plate at 2-cm intervals. We quantified v T by measuring the displacement of the trench in the center of the subduction zone for successive photographs, and quantified v SP by measuring the displacement of the subducting plate trailing edge. The trenchward motion of the subducting plate and trench retreat are defined as positive here. We also spread black tracers randomly around the subducting plate to visualize the mantle flow at the surface during subduction.

Results
Our 12 analog models show two main different subduction styles. To qualitatively and quantitatively evaluate the subduction behavior in our models we focused on (1) the subduction evolution and kinematics; (2) the slab bending radius, slab bending force, and basal drag force; and (3) the trench curvature.

General Subduction Evolution and Kinematics
In all experiments we can distinguish three phases during the subduction process, similar to previous studies (e.g., Funiciello et al., 2003;Schellart, 2008;Strak & Schellart, 2014). First, there is the "free sinking" phase, second the "transitional" phase and lastly the "steady-state" phase. During the free sinking phase, all experiments experience fast subduction with a drastic increase of all velocities (v SP , v T , and v S ) until the slab tip reaches the bottom of the tank. This first phase is characterized by trenchward motion of the subducting plate and retreat motion of the trench. A short transitional period follows during which the slab tip interacts with the 660-km discontinuity and all velocities decrease. This also results in geometrical changes of the slab as it rolls over and drapes forward forming a "U" geometry rotated 90° (Figures 3a-3c) or folds and drapes backward forming a lazy "S" shape (Figures 3d-3l). During the subsequent steady-state phase, two different styles of trench motion are observed. They are (1) the "advancing" mode with a slab moving toward the mantle wedge and a rollover slab geometry (Figures 3a-3c) and (2) the "retreating" mode with a slab rolling backward to the subducting plate side and a slab draping geometry lying flat below the mantle wedge (Figures 3d-3l). The "advancing" mode is observed in the lower viscosity ratio models with a shorter plate (plates of 20, 30, and 40 cm long). The "retreating" mode is observed in the lower η SP /η M models with a longer plate (50, 60 cm long plates and fixed trailing edge model) and in all higher η SP /η M models. 3.1.1. Low η SP /η M Models (Experiments 1-6) The velocities (v SP , v T , and v S ) in the experiments vary with plate length (Figures 4a-4d). The low η SP /η M Experiments 1, 2, and 3 are characterized by an "advancing" mode and have an overall progressive increase of the subducting plate velocity v SP with a temporary drop during the second phase (Figure 4a). A local peak of v SP of 0.13, 0.065, and 0.04 mm/s is reached at the end of the free sinking phase for Experiments 1, 2, and 3, respectively. During the transitional phase, v SP decreases temporarily and it increases again continuously during the final phase with peak velocities of 0.17, 0.15, and 0.12 mm/s at the end of subduction. For the "advancing" models the fastest v SP with the sharpest increase is obtained for the shortest plate (Exp. 1, Figure 4a). The "retreating" models Exps. 4 and 5 experience a rapid increase of v SP only during the first phase. The maximum v SP is 0.038 and 0.034 mm/s and occurs at the end of the free sinking phase. There is a decrease in v SP during the transitional phase, but the experiments keep a relatively steady v SP ranging between 0.012 and 0.018 mm/s during the steady-state subduction phase. The difference in v SP is negligible between the "retreating" models. For Exp. 6, we calculate that v SP due to the trench-normal surface extension of the subducting plate during subduction is only 4.7% of the average v T . Thus, considering that the contribution of v SP to v S is negligible for Exp. 6 with the fixed trailing edge, we assume that v SP = 0 during the whole experiment. The trench velocity v T (Figure 4b) evolves similarly for all lower η SP /η M models during the free sinking phase, with a local maximum between 0.05 mm/s (Exp. 1) and 0.09 mm/s (Exp. 6). After the transitional phase, v T for the "advancing" models with relatively short plates (Exps. 1-3) decreases rapidly and attains negative values, whereas v T for the "retreating" models with relatively long plates (Exps. 4-6) remains relatively constant during the steady-state subduction phase.
The relationship between v SP , v T and plate length during the free sinking phase is investigated statistically in this study ( Figure 5). We take the maximum value for both v SP and v T during this phase to calculate the least squares linear best fit trend lines and correlation coefficients (Pearson) and plot them as a function of plate length. The low η SP /η M experiments ( Figure 5a) show a strong negative linear correlation between plate length and v SP and a strong positive linear correlation between plate length and v T , with linear regression correlation coefficients R = 0.90 and 0.86, respectively.
The subduction velocity (Figure 4c) shows minor variation between models with different plate lengths during the free sinking and transitional phases. The peak v S at the end of the free sinking phase (Figure 5a) shows an approximately negative linear relationship with plate length, with R = 0.76. However, the difference between the highest (v S = 0.13 mm/s from Exp. 1) and lowest value (v S = 0.096 mm/s from Exp. 4) is moderate (Δv S = 0.034 mm/s). During the steady-state subduction phase (Figure 4c), the "advancing" mode models (Exps. 1-3) display a steady increase of v S and the shorter plates (20-30 cm) show a steeper increase in v S than the longer plate (40 cm). In the "retreating" mode models (Exps. 4-6) v S remains relatively steady overall during the steady-state phase, with some temporary deviations.
The ratio v SP /v S represents how v S is partitioned into v SP and v T , where a value of 0.5 represents equal partitioning between v SP and v T , and values >1 indicates trench advance. Values for models with plate lengths of 20-40 cm increase quickly to a peak v SP /v S ranging between 1.6 and 4 until the end of the free sinking phase, while afterward it remains stable at about 1.6 during the steady-state subduction phase (Figure 4d). For the models with a plate length of 50-60 cm, the v SP /v S also increases to a peak value at the end of the free sinking phase, but these peak values are significantly lower (~0.4-1.4), and the values drop significantly after the transitional phase fluctuating between 0 and 0.5.
The average values of v SP , v T , v S , and v SP /v S during the steady-state subduction phase for all low η SP /η M models are illustrated in Figure 6a. The average v SP shows a decrease with increasing subducting plate length, keeping a positive value, but there is a notable drop in velocity from the model with a 40-cm-long subducting plate producing a rollover subduction style to the model with a 50-cm-long plate that leads to a rollback subduction style. In contrast, there is a strong increase in v T between these two models when plate length is increased. In general, there is a positive correlation between v T and subducting plate length, with negative trench velocities (trench advance) for those models with a rollover subduction style. For the "retreating" mode models, v T is positive and increases with plate length. v S is positive for all models and slightly decreases with plate length for "advancing" mode models, and slightly increases with plate length for "retreating" mode models. v SP /v S shows no general trend as a function of plate length, but a nearly constant higher value for all the "advancing" mode models and nearly constant lower value for all the "retreating" mode models.

Journal of Geophysical Research: Solid Earth
For all these average values, there is a clear difference between the "advancing" mode models (Exps. 1-3) with shorter plates up to 40 cm, and "retreating" mode models (Exps. 4-6) with a minimum plate length of 50 cm. 3.1.2. High η SP /η M Models (Experiments 7-12) For the high η SP /η M models, only the "retreating" mode with a slab draping geometry is observed. Experiments 7-11 with a free trailing edge experience a rapid increase of v SP during the free sinking phase, reaching a peak v SP at the end of the free sinking phase (Figure 4e). The peak value of v SP decreases with increasing plate length, with Exp. 7 experiencing the highest (0.08 mm/s) and Exp. 11 the lowest v SP (0.04 mm/s). A clear negative linear correlation between the peak v SP and subducting plate length is found during the free sinking phase from this set of experiments with relatively high η SP /η M with R = 0.98 (Figure 5b). During the transitional phase, v SP temporarily decreases to negative values in Experiments 7-11, indicating that the subducting plate is moving backward (away from the trench). The subducting plate velocity increases again after the transitional phase before reaching a constant value during the steady-state phase. There is no obvious difference in steady-state v SP between the models as very similar values are obtained for the average v SP , ranging between 0.01 and 0.02 mm/s (Figure 6b).
The trench velocity v T for all high η SP /η M models (Figure 4f) increases to a local maximum at the end of the free sinking phase. There is a general increasing trend of the peak value of v T with increasing plate length. The trench velocity has a value of~0.05 mm/s for Exps. 7 and 8, whereas for Exp. 12 v T is~0.11 mm/s. A roughly linear positive relationship between v T and plate length is found (Figure 5b), with R = 0.89. During the first phase, v T is lower than v SP for models with shorter plates (Exps. 7 and 8), and higher than v SP for models with longer plates (Exps. 9-12). The transitional phase is characterized by a rapid decrease of v T until the subduction process reaches the steady-state phase. During the steady-state phase, models with a plate length of 20-50 cm (Exps. 7-10) show a minor increase of v T , whereas the models with a plate length of 60 cm with a fixed and free trailing edge show a minor increase followed by a minor decrease in v T . Importantly, v T is higher than v SP for all models during this phase.
The peak v S at the end of the free sinking phase was plotted against plate length (Figure 5b), and it shows almost no correlation (R = 0.01) and minor differences between different models (the difference between two extreme values Δv S = 0.01 mm/s). The peak v S shows some variability during the steady-state subduction phase, but, except for the fixed trailing edge experiment, those with a free trailing edge all show comparable values (Figures 4g and 6b).
During the free sinking phase, v SP /v S is between 0.1 and 0.62 for all high η SP /η M models (Figure 4h). Except for Exp. 8, the v SP /v S decreases rapidly during the transitional phase to a minimum of~−0.9−0.3 before reaching the steady-state subduction phase. Exp. 8 shows a quick increase of v SP /v S to 1.17 and then a decrease to negative values before reaching the steady-state subduction phase, which means subduction experienced a brief period of slight trench advance (between 7 and 8 min after subduction initiation) Figure 6. Diagrams showing the average velocities and average v SP /v S during the steady-state subduction phase for low (a) and high (b) η SP /η M models. The bars show the range of the increasing v SP and v S , and decreasing v T and v SP /v S of the shorter and advancing mode plates, corresponding to the low η SP /η M . Average values of the constant v SP , v T , v S , and v SP /v S were plotted for the longer, retreating mode plates and higher η SP /η M models. Note that v T = v S for models with fixed trailing edge; hence, v S overlaps with v T in the graphs for models with fixed trailing edge.
shortly after the slab tip reached the bottom boundary. During the steady-state subduction phase, v SP /v S remains relatively stable between 0 and 0.5 for all high η SP /η M models and shows little variation between the different models ( Figure 6b).
The v SP , v T , v S , and v SP /v S all show comparable averaged values during the steady-state subduction phase for all high η SP /η M models with a free trailing edge (Figure 6b), and there is no apparent dependence on subducting plate length. The model with fixed trailing edge shows higher v T and v S , while v SP and v SP /v S both equal to 0 because v SP due to trench-normal extension in the subducting plate is negligible (it is onlỹ 1.8% of the average v T ).

Slab Bending Radius, Bending Forces, and Basal Drag Forces
We measured the slab bending radius at the subduction zone hinge (r B ) and slab tip dip angle when it first touches the bottom (θ B ) for all models (Figure 7). The slab bending radius was measured when the slab tip reached a depth of 6.25 cm, 2 cm above the bottom of the model box. We choose this advanced stage of the free sinking phase, because the slab has enough time to evolve self-consistently, canceling the initial bending radius imposed at the start of the experimental run, while the slab is not much affected yet by the bottom boundary. The measurement was performed by manually fitting a circle to the outer surface of the bending part of the slab from the trench to a depth of 4.25 cm (Figure 7c), which is the depth range over which the slab experiences its greatest curvature, and herein the radius of the circle represents the slab bending radius, which is similar to the "depth" method described in Petersen et al. (2017). The slab bending radius varies between the different models ( Figure 7a). A higher bending radius was observed in the higher η SP /η M models compared to the lower η SP /η M models. Both lower and higher η SP /η M models with the free trailing edge (Exps. 1-5 and 7-11, respectively) show a general increase in bending radius with increasing plate length. For the 60 cm plate models with a fixed trailing edge (Exps. 6 and 12), the bending radius is slightly lower compared to that of the 60 cm plate models with a free trailing edge (Exps. 5 and 11). Both sets of experiments show a relatively strong positive correlation between the slab bending radius and the plate length with R = 0.90 for the low η SP /η M models and R = 0.80 for the high η SP /η M models. The slab tip dip angle was measured when the slab tip reached the bottom boundary ( Figure 7b). A tangent line was drawn from the slab tip along the surface of subducting plate from the bottom boundary up to the second tracer, and the angle between the tangent line and the bottom boundary was measured (Figure 7d). The θ B shows a negative correlation with plate length for both sets of experiments, which is quite strong for both with R = 0.91 and R = 0.94 (Figure 7b).
Here, we also calculated the slab bending resistance (F Be ) for all models when the slab tip reached the position in Figure 7c. We use the following equation to calculate F Be (Buffett, 2006;Conrad & Hager, 1999;Irvine & Schellart, 2012;Schellart, 2009): where A C is the ratio of slab curvature length to bending radius. We measured A C and found that it ranges between 1.2 and 1.7, which is comparable to values of A C = 1.4-1.9 measured by Schellart (2009) for a similar subduction stage, but somewhat smaller than the value of A C = 2 proposed by Buffett (2006). T SP represents the thickness of the subducting plate, r B is the slab bending radius, and W represents the width of the subducting plate.
We also calculated the buoyancy force of the slab (F Bu ): where L slab is the length of the slab, Δρ represents the density contrast between the slab and the upper mantle, and g is the gravitational acceleration.
We have calculated the F Be /F Bu ratio for all the experiments (Figure 8), which shows that it falls between 2.8% and 10% for the low η SP /η M models, and between 4.3% and 7.3% for the high η SP /η M models. The value F Be /F Bu is the highest and second highest in Exps. 1 and 2 with shorter plates and low η SP /η M , respectively, reaching 10% and 8%.
In addition, we have calculated the ratio (F D /F Bu ) of viscous basal drag force (F D ) to slab negative buoyancy force (F Bu ) for Exps. 1-5 and 7-12. We calculate F D using the following equation: (7) where L SP represent the length of the subducting plate and γ is the horizontal trench-normal shear rate in the sublithospheric mantle located between the base of the plate and the bottom of the tank. Since the upper mantle is a linear viscous material, we can derive γ from the difference between the subducting plate velocity (v SP ) and the velocity at the bottom of the tank (v bottom ): where d is the vertical distance between the base of the subducting plate and the bottom of the tank, and v bottom is assumed to be 0.
We calculate F D /F Bu in two different ways: (1) During the free sinking phase when the slab tip is at 6.25-cm depth, that is, 2 cm above the base of the tank (the stage shown in Figure 7c).
(2) For the free sinking phase using an average v SP and average L slab . Figure 9 shows the results of our calculations, showing a positive correlation between the effective plate length (L SP -L slab ) and F D /F Bu for both experimental sets and for both approaches, indicating that the basal drag force increases with increasing subducting plate length. The correlation for the high η SP /η M set of models is high (R = 0.99) for both approaches. For the low η SP /η M set of models, the correlation is also high (R = 0.97) for the second approach (Figure 9b) but is lower (R = 0.8) and has a less significant slope for the first approach (Figure 9a).

. Trench Curvature
We measured the trench curvature in all models (Figures 10a-10d) by fitting the trench with a circular segment and define the ratio (h/s) of segment height (h) and the chord of the segment (s) as the trench curvature (Figure 10e), following the method in Peral et al. (2018). The trench attains a concave geometry toward the mantle wedge in all experiments (Figures 10a and 10b). The lower η SP /η M models display a higher range of trench curvature (Figure 10c) ranging from 0.07 to 0.25 compared with the higher η SP /η M models (trench curvature ranges from 0.08 to 0.14) (Figure 10d).
The lower η SP /η M models show that Exps. 4-6, which are characterized by trench retreat, have a more curved trench geometry than Exps. 1-3, which are characterized by trench advance. Experiments 2 and 3 both show a decrease of the trench curvature after~20 and 30 min of the subduction. Experiments 4-6 show a somewhat variable behavior, but all show that, over the duration of the measurement period, the trench curvature has increased. The higher η SP /η M models generally display a progressive increase in trench curvature ( Figure 10d).

Effects of Subducting Plate Length on Subduction Style
In the experimental set with the low η SP /η M of~740, two different subduction styles were observed from two groups of models with shorter and longer plates. A reversed trench motion from retreat to advance has been documented in the models with shorter plates of 20, 30, and 40 cm, whereas continuous trench retreat was observed for models with plate lengths of 50-, 60-, and the 60-cm-long plate with fixed trailing edge. In the experimental set with a higher η SP /η M of~1,680, only one subduction style was found, namely, trench retreat and a draping slab geometry.
The experimental results can be explained by the forces involved in the subduction system. In nature and our experiments, the negative buoyancy force of the slab F Bu is the main driving force of subduction, and the main resistive forces are the viscous drag forces in the mantle, the bending resistance at the trench and the shear resistance at the subduction zone interface. Our results ( Figure 8) and earlier studies (Capitanio et al., 2007(Capitanio et al., , 2009Irvine & Schellart, 2012;Schellart, 2009) show that the bending resistance only accounts for a small part of F Bu . The shear resistance at the subduction zone interface (F SI ) is also thought to be low (e.g., Duarte et al., 2013). Thus, most of F Bu (>~70%) is dissipated in the sublithospheric mantle, as also concluded in other studies (Capitanio et al., 2009;Chen et al., 2015;Schellart, 2004b;Stegman et al., 2006), both to drive toroidal flow and to drive poloidal flow. Of the poloidal resistance, one component includes the viscous basal drag force of the subducting plate (F D ). Resistance to trenchward subducting plate motion increases linearly with effective subducting plate length (L SP − L slab ), because F D scales linearly with (L SP − L slab ) as shown in Equation 7. Furthermore, the experimental results in Figure 9 show a positive Figure 9. Ratio of viscous drag force (F D ) to slab negative buoyancy force calculated for Exps. 1-5 and 7-12 using two approaches, (a) F D /F Bu calculated during the free sinking phase when the slab tip is at 6.25 cm depth, that is, 2 cm above the base of the tank (the stage shown in Figure 7c). (b) F D /F Bu calculated for the free sinking phase using an average v SP and average L slab . For each data set a least squares linear best fit line has been plotted and a correlation coefficient (R) has been given.

10.1029/2020JB020514
Journal of Geophysical Research: Solid Earth linear correlation between (L SP − L slab ) and F D /F Bu . However, the slopes of the best fit lines are lower than one might expect from Equation 7. Indeed, a doubling of (L SP − L slab ) does not correspond to a doubling of F D /F Bu . This is because v SP (and thus shear rate) is not constant and an increase in (L SP − L slab ) comes with a decrease in v SP .
Thus, during the free sinking phase, lower resistance to trenchward subducting plate motion experienced by shorter plates promotes faster v SP (Figure 4a), and vice versa for models with longer plates. However, subducting slabs try to sink into the mantle at their Stokes' sinking velocity, which means slabs with long and short trailing plates want to sink at the same rate. This can be seen from the relatively comparable values for v S during the free sinking phase for models with different plate lengths (Figures 4c, 4g, and 5). Moreover, slab subduction and sinking can occur either through subducting plate motion (v SP ) or trench motion (v T ), with v S = v SP + v T . Thus, for long subducting plates, F D is relatively high, promoting a relatively low v SP and high v T , while for short plates, F D is relatively low, promoting a relatively high v SP and low v T . Earlier studies (Griffiths et al., 1995;Schellart, 2004a) have shown that the upper mantle slab dip angle has a strong correlation with v T , where a high v T promotes a lower dip angle and a low v T promotes a higher dip angle. Thus, a long plate, with a relatively low v SP and high v T , will attain a relatively low slab dip angle and a relatively high bending radius. In such a scenario, the slab tip will touch the bottom boundary with a relatively low dip (in our experiments between 83 and~97°, Figure 7b), which promotes the slab to subsequently fall

10.1029/2020JB020514
Journal of Geophysical Research: Solid Earth backward and roll back. In contrast, a short plate will attain a relatively steep slab dip, and, in case the slab does not have enough time to unbend in the middle-upper mantle, a slab with a steep dip will touch the bottom boundary with a relatively high (overturned) tip (in our experiments between~108 and 126°, Figure 7b), which promotes the slab to subsequently fall forward and roll over.
On the other hand, all models with a high η SP /η M show only trench retreat and backward slab draping. This phenomenon can be explained by the high η SP /η M , which causes a higher resistance to bend at the trench, giving rise to a larger bending radius ( Figure 7a) and a lower dip angle at the slab tip once the slab reaches the bottom boundary (Figure 7b). Even though v SP is higher than v T for shorter plates in the high-viscosity ratio models during the free sinking subduction phase, the smaller slab dip angle (~82-95°; Figure 7b) that forms once the slab reaches the bottom boundary is not sufficient to bend the slab in a U shape on its side and form a rollover slab geometry (Figure 7b).
All the models in the high η SP /η M set of experiments have a negative value of v SP during the transitional phase (Figure 4e), which means the subducting plate moves away from the trench during the interaction between the slab tip and the bottom boundary. This can be explained by the relatively high angle between the slab tip (~90°) and bottom boundary and the high resistance to buckle and fold at the bottom boundary, which is caused by the high η SP /η M . The high resistance to fold causes the slab to retain a relatively planar geometry, and, to continue sinking, requires it to fall backward (toward the subducting plate side), thereby applying a push force to the trailing subducting plate at the surface, which causes the subducting plate to move away from the trench.
We observe that experiments showing trench advance during the steady-state phase (Exps. 1-3) all show a continuously increasing v SP during this phase, whereas experiments showing trench retreat (Exps. 4-5 and 7-11) all show an almost constant v SP during the steady-state phase. We explain it as follows. Once the slab tip has reached the bottom boundary, it is facing backward (toward the subducting plate side) for the rollover mode and facing forward (toward the mantle wedge side) for the rollback mode. The bottom boundary is a no-slip boundary, so for the slab to slide/shear horizontally over it requires a significant shear force to be overcome. For the rollover mode, however, such horizontal sliding is not required during trenchward subducting plate motion, because v SP is mostly accommodated at depth by the forward rolling of the slab on top of the base (see, e.g., Figures 9d and 10d in Schellart, 2008). For the rollback mode, such horizontal sliding is very much required to accommodate trenchward subducting plate motion (or else the slab has to thicken or buckle at the bottom boundary to accommodate v SP ). Thus, during rollover subduction v SP is not hindered by shear stresses at the slab-bottom boundary contact area but mostly by shear stresses at the plate-syrup contact area at the base of the subducting plate, which explains why v SP progressively increases with decreasing plate length during a subduction experiment. In contrast, for the rollback mode v SP is affected both by the slab-bottom boundary contact area and the plate-syrup contact area. Considering that a decrease in plate-syrup contact area is accommodated by an equal increase in slab-bottom boundary contact area, the average resistance to v SP remains approximately constant during progressive subduction, and so a decrease in plate length does not result in an increase in v SP . Thus, in the high η SP /η M models, which all show a rollback slab geometry, v SP is independent of plate length and relatively low due to the high shear resistance at the bottom of the tank. This causes subduction to be accommodated mostly by slab rollback and v T to be relatively high for all the high η SP /η M models and causes v T and the subduction partitioning ratio to be independent of plate length as well, as shown in Figure 6b.
The evolution of trench curvature can be explained by the quasi-toroidal mantle flow, which is driven by the lateral migration of the slab through the upper mantle. Here we describe this mantle flow as quasi-toroidal, because it is not strictly toroidal, in that it has a vertical velocity component as well (see also Schellart, 2004aSchellart, , 2008Strak & Schellart, 2014). During the free sinking phase, the slab rollback activates the quasi-toroidal mantle flow and produces a concave trench geometry toward the mantle wedge. If the trench continues to retreat, the trench geometry becomes more curved by the continued quasi-toroidal mantle flow, in agreement with earlier works (Loiselet et al., 2009;Morra et al., 2006;Schellart, 2004aSchellart, , 2010. On the other hand, for the models with trench motion reversed to advance, the quasi-toroidal mantle flow reverses with mantle material moving from the mantle wedge side, around the lateral slab edges, toward the subslab side. This opposite quasi-toroidal mantle flow explains the progressively reduced original trench curvature in models with trench advance, as also documented in Schellart (2010). In addition, our models with a high η SP /η M show a less curved trench compared to the models in the same mode but with a low η SP /η M . This can be explained by the higher stiffness of the plate that reduces the deformation caused by the quasi-toroidal mantle flow, in agreement with earlier works (Loiselet et al., 2009;Morra et al., 2006;Schellart, 2010).

Comparison With Previous Work
Our results show that shorter plates have a higher v SP during the free sinking phase, which promotes trench advance with a rollover slab geometry. Previous studies have also shown that long-term trench advance is associated with a slab rollover geometry (Bellahsen et al., 2005;Di Giuseppe et al., 2008;Funiciello et al., 2008;Schellart, 2005Schellart, , 2008. Schellart (2005) varied v SP by exerting an external velocity boundary condition at the plate trailing edge and found that a lower v SP promotes trench retreat and a higher v SP promotes trench advance. Bellahsen et al. (2005) studied several factors that influence the subduction style, including the thickness, viscosity and width of the subducting plate, as well as thickness of the mantle. They found that these parameters affect subduction style by controlling the partitioning of v S into v SP and v T , and the slab bending radius. Funiciello et al. (2008) also suggested that v SP contributes to produce models showing trench advance and that v T is higher for models showing trench retreat. Faccenna et al. (2007) showed that trench advance is associated with a higher slab dip and higher v SP and that trench retreat is related to lower slab dip and lower v SP , and they also found that there is an inverse correlation between v T and v SP . Our results show the same correlation between velocities and the subduction styles, namely, that trench advance correlates with a higher v SP and lower v T and trench retreat shows a lower v SP relative to v T . Schellart (2011) presented a regime diagram showing how the partitioning of v S influences the subduction mode, with v SP /v S ≤ 0.5 resulting in a slab draping geometry and v SP /v S ≥ 1.5 resulting in a rollover slab geometry. Our experimental results agree with this earlier work, where the "advance" mode models (Exps. 1-3) show an overall v SP / v S > 1.5 during the steady-state subduction phase, while the models with a retreat mode (Exps. 4-12) all have a v SP /v S ≤ 0.5 (Figures 4d, 4h, and 6).
The calculations of F Be /F Bu show that F Be accounts for a small part of F Bu , which is between 2.8% and 10%. This agrees well with earlier studies (Capitanio et al., 2007(Capitanio et al., , 2009Irvine & Schellart, 2012;Schellart, 2009). For models showing exclusively trench retreat, the results show that a higher F Be /F Bu is coupled with a higher η SP /η M , in agreement with earlier works (Schellart, 2009). The calculations of F D /F Bu show that F D accounts for a small part of F Bu and ranges between 0.5% and 1.6%. This is even less than an estimated 8% from Schellart (2004b) for a similar subduction stage. The low values are generally consistent with the work of Stegman et al. (2006), who calculated a value of 13% for all the poloidal flow in the mantle, which includes the poloidal flow component related to F D . Schellart (2008) proposed a subduction regime diagram in which two parameters control the subduction kinematics and the subduction style, namely the slab-to-mantle viscosity ratio (η SP /η M ) and mantle-to-plate thickness ratio (T M /T SP ). Both parameters work together to control the slab bending angle before the slab tip reaches the transition boundary, which controls subsequent development of the trench kinematics and slab geometry. This slab bending angle can be represented by the ratio of slab bending radius to mantle thickness (r B /T M ). Based on r B /T M , four regimes of subduction styles were distinguished (Figure 11a). According to this regime diagram, the results from all of our models should be positioned in either Regime I or IV (Figure 11a) that both are characterized by trench retreat. However, our Exps. 1-3 with low η SP /η M show an "advance" mode, with the slab geometry falling into Regime III (rollover) in Schellart's (2008) diagram.
Other models  from this study all show exclusively a "retreat" mode, which agrees well with the regime diagram with results falling into Regime I or IV. This discrepancy could possibly be (partly) explained by the fact that Schellart's (2008) diagram was based on a limited amount of experiments, so that the actual locations of the boundaries separating the different regimes could be shifted. More importantly, the experiments chosen for Schellart's (2008) diagram are from models with a different plate length, the influence of which was not investigated in this diagram. The different lengths that were used in earlier studies are listed in Table 1. What is important to note here is that the plate lengths used in the regime diagram of Schellart (2008) (2008) using 2-D numerical models with the boundaries between regimes somewhat shifted, but otherwise looking very comparable. According to Ribe's diagram (Figure 11b), all of our experiments should display the FR (Folding-Retreat) mode. However, our models show either a retreat mode or an advance mode. The difference between our results and Ribe's 2-D numerical experiments can be ascribed to the following differences in model conditions and setup: (1) Ribe (2010) used 2-D numerical models, implying an infinitely wide plate and slab compared to a narrow plate and slab (10 cm scaling to 800 km in nature) in our experiments.
(2) Ribe (2010) used an infinitely deep mantle reservoir, while we used an upper mantle reservoir with a rigid 660-km discontinuity.
(3) Surface tension only exists in the analog experiments, which reduces v T , but not in the numerical models. In addition, Ribe (2010) concluded that plate length does not impose a significant difference on the slab geometry. From inspecting Figures 2b and 2c in Ribe (2010), a small difference in slab tip dip angle between models with a different plate length can be observed. Indeed, a longer plate has a nearly vertical slab tip dip angle and the shorter plate has a slightly overturned slab geometry. Although the difference in slab tip dip angle is minor, it can still determine the subsequent subduction kinematics and evolution. In addition, only two plate lengths were compared in the work from Ribe (2010), which correspond to 2,000 and 3,600 km respectively in nature. In our current study, we investigated a larger range of plate lengths (1,600-4,800 km in nature), and the subduction mode changed in the low η SP /η M models when the plate length was changed from 3,200 to 4,000 km. We suggest that if a larger range of plate lengths was applied, a different subduction mode may have occurred in the numerical models of Ribe (2010).
Li and Ribe (2012) performed a large number of 3-D numerical simulations varying different physical parameters, a number of which have comparable parameters as our models in terms of tank depth to plate thickness ratio (6.8 compared to 6.7 in our experiments), η SP /η M values of 600 and 1,000 (compared to~740 in our Exps.1-6), and a plate length to plate thickness ratio of 24 (compared to 25 in our Exp. 2). These models from Li and Ribe (2012) all show a retreat/slab rollback subduction mode (their Figure 11b). Our Exp. 2, however, shows a rollover slab geometry, which we ascribe to the very short subducting plate length. The cause for the difference is not immediately clear as Li and Ribe (2012) use a comparably short subducting plate, but could be due to the different initial slab dip angle (30°vs. 35-40°in our models), the free surface in our models vs. a lubrication top layer in Li and Ribe (2012) and the difference in slab width (a slab width to plate thickness ratio of 12 vs. 6.7 in our models).
The influence of subducting plate length on subduction kinematics and slab dip angle was studied by Holt and Becker (2017) using 2-D and 3-D numerical models with different plate lengths of 2,000, 3,000 and 4,000 km. They found that longer plates experienced higher resistance to subducting plate motion, producing higher v T /v SP ratios and lower slab dip angles and promoting trench retreat, and vice versa. Our results are consistent with these results. However, none of the short plate models from Holt and Becker (2017) has Figure 11. Regime diagrams from Schellart (2008) (his Figure 13) (a) and from Ribe (2010) (his Figure 11) (b). Naming of subduction styles use regime numbers from Schellart (2008) and the naming from Ribe (2010) for (a) and (b). "WR" represents weak retreating mode. "FR" represents folding retreating mode. "A" represents advancing mode. "SR" represents strong retreating mode. The two red dots show the position on the diagrams of the two sets of experiments reported in this study.
shown long-term trench advance and a rollover slab geometry. This can be caused by the difference in model parameters used in Holt and Becker (2017) and this study, including plate width (unlimited in 2-D and 2,000 km in their 3-D models compared to 800 km in this study), box boundary conditions (free slip boundaries in their models vs. no-slip and free surface in this study), initial slab dip angle (70°in their vs. 35-40°in our models) and plate thickness (80 km in their compared to 96 km in our models).

Implications for Subduction Systems in Nature
Our experiments show a clear negative correlation between subducting plate length and v SP as well as subducting plate length and subduction partitioning. For a comparison between our experimental results and nature, we have chosen only those subducting plates that are relatively simple, like our model subducting plates, in that they contain one subduction zone, that this subduction zone forms some 30-50% of the plate's circumference, and that the plate does not border a collision zone. We have used the subduction zone data presented in (Schellart et al., 2007 to check for these criteria, which has led us to exclude a significant number of plates, including the Pacific, Eurasian, Australian, Indian, African, Arabian, North American, South American, Caribbean, and Antarctic plates. For example, we exclude the Pacific plate, because it is attached to three subduction zones (Aleutians-Alaska, Kamchatka to Mariana, and Tonga-Kermadec-Hikurangi). There are three plates that fit these criteria very well, namely, the Nazca, Cocos, and Juan de Fuca plates. We have also included the Philippine plate, which is attached to the Nankai-Ryukyu subduction zone. This plate also borders the Philippine Trench, but this trench is classified as an incipient subduction zone in Schellart et al. (2007) and therefore fits our criteria. Table 3 lists a number of characteristics of these plates and their subduction zones, including subduction zone width, average subducting plate length, and average subduction partitioning ratio. It is generally known that the absolute reference frame affects v SP , v T and subduction partitioning. We use the Indo-Atlantic hotspot reference frame (O'Neill et al., 2005), as this reference frame best agrees with a number of observations (Muller & Roest, 1992;Schellart, 2011) and with the concept of minimizing trench migration and viscous dissipation in the Earth's mantle . Apart from the Cascadia subduction zone, the other three show that decreasing plate length corresponds with increasing subduction partitioning, which is consistent with the experimental results. What is not consistent is that none of these subduction zones show a rollover slab geometry, while our models with low viscosity ratio and short plate length do show a rollover slab geometry. In addition, one may expect that, from the extremely short length of the Juan de Fuca plate, the Cascadia subduction zone should show a high partitioning ratio and a rollover slab geometry. The discrepancy between models and observations likely stems from several other physical parameters that are different between our experimental and natural subduction zones. For example, the natural subduction zones are much wider (1,400-7,400 km) than our experimental subduction zones (scaled to 800 km), and previous works have shown that slab width affects subduction partitioning and slab geometry . Additionally, the Mexico-Central America, Cascadia, and South America subduction zones are very old and have slab segments that penetrated into the lower mantle (Amaru, 2007;Gorbatov & Fukao, 2005;van der Meer et al., 2018), while our experimental rollover geometry forms in an early stage during upper mantle subduction, and our models do not allow for lower mantle Note. Subduction zone width and subduction partitioning are from Schellart et al. (2007Schellart et al. ( , 2010. Average subducting plate lengths are calculated from subducting plate surface area divided by subduction zone width. subduction due to the rigid 660-km discontinuity. In addition, the viscosity ratio also plays a role determining the subduction partitioning and slab geometry, as shown by our experiments and previous works Ribe, 2010;Schellart, 2008), and there is still a debate concerning the viscosity ratio in natural subduction systems. Conrad and Hager (1999) also found that there is no clear trend between plate length and v SP for natural examples. This can be explained from the fact that the kinematics and geometry of subduction zones depend on interactions between many physical parameters. Our study is a parametric study focusing mainly on the effect of plate length, which hinders the systematic comparison between models and nature because in nature, many physical parameters vary between different subduction zones, including slab width, slab age, slab thickness, subduction zone age, slab depth extent, overriding plate type (oceanic/continental), and lateral homogeneity/heterogeneity of the subducting plate. As such, we expect that it is difficult to extract a clear correlation between subducting plate length and v SP from observational data due to the many other parameters that also affect v SP .
Our models do show that a shorter plate promotes a higher v SP and trench advance, which facilitates occurrence of the rollover slab geometry. According to tomography data on slab morphology, most subduction zones in nature show a rollback slab geometry or steep slab geometry (e.g., Fukao et al., 2009;Li et al., 2008;Simmons et al., 2012;van der Meer et al., 2018), with only the slab at the India-Eurasia continental subduction zone showing a rollover slab geometry ( Van der Voo et al., 1999) although the India-Eurasia collision zone is relatively wide (~2,700 km wide) and the subducting Indian plate is relatively long. At the earlier stage of Neo-Tethys subduction northward under the continental Eurasian plate between~70 and 45 Ma, the subducting plate moved very fast, with a speed ranging between 10 and 20 cm/year (White & Lister, 2012). In addition, tectonic reconstructions show that the India-Eurasia collision boundary has moved northward (Replumaz et al., 2004;Replumaz & Tapponnier, 2003;Schellart et al., 2019), which indicates trench advance and implies slab rollover. Particularly, recent work (Schellart et al., 2019) estimated that the India-Eurasia collision zone has migrated approximately northward some 1,700 km over the last 50-55 Myr, with an average "trench" advance rate of about 3.2-3.5 cm/year. With an average convergence velocity of~5 cm/year over the last~50 Myr, this implies an average subduction velocity of 1.5-1.8 cm/year, and an average subduction partitioning of~3. Our experimental results indicated that such a high partitioning ratio coincides with trench advance and a rollover slab geometry, which is consistent with the observed rollover slab geometry (Figure 1b). The Indian plate has a plate length of the order 3,000-4,000 km, which is comparable to that of our Experiments 3 with rollover and 4 with rollback. As such, its moderate plate length could have partly played a role in promoting the fast plate motion, trench advance and rollover slab geometry. Nevertheless, as discussed in previous sections and as shown in previous works, the subducting plate does not have to be very short and a variety of plate lengths can coexist with a rollover geometry, since other physical parameters also provide a control on subduction kinematics and slab geometry, such as η SP /η M , slab width, mantle stratification, and upper mantle rheology Ribe, 2010;Schellart, 2008;. External forcing to drive rapid subducting plate motion such as plume push (Cande & Stegman, 2011) may also play a role. Therefore, we propose that as long as v SP /v S >~1.5 and sustains over several tens of million years, a slab rollover geometry will form.

Conclusions
Three-dimensional fully dynamic laboratory experiments of subduction were conducted at two different subducting plate to mantle viscosity ratios η SP /η M to study the effect of plate length on subduction behavior. Our geodynamic modeling results indicate that variation in plate length is one of the factors that can determine the style of subduction by controlling the partitioning of v S into v SP and v T during the free sinking phase, and thereby providing a control on the angle of the slab tip as it first touches the 660-km discontinuity (Figures 5a and 7b). The experimental results show two different subduction styles in models with a lower η SP /η M , whereas only one style is produced in models with a higher η SP /η M . A shorter plate with low η SP /η M favors a higher v SP and lower v T , promoting a larger (overturned) slab tip angle when it first touches the 660-km discontinuity, thereby causing the subduction system to adopt a rollover slab geometry with trench advance. However, if η SP /η M is higher, the slab bending is lower (larger bending radius) during the free sinking phase (Figure 7a) because the subducting plate's resistance to bending is higher. For these 10.1029/2020JB020514 Journal of Geophysical Research: Solid Earth models, although shorter plates produce a higher v SP , higher resistance by the high η SP /η M facilitates a larger bending radius, resulting in a lower slab tip dip angle when the slab tip approaches the 660-km discontinuity, thereby preventing the slab to form a rollover slab geometry and to advance, and instead forcing the slab to roll back and form a slab draping geometry.