Using open sidewalls for modelling self-consistent lithosphere subduction dynamics

Introduction Conclusions References


Introduction
In the past decades, numerical modelling of lithosphere subduction has advanced considerably by incorporating coupling between plates, between plates and mantle, and by incorporating the complexity of detailed subduction zone processes (see Gerya, 2011, for a review and references therein).Up to now modelling of regional subduction evolution is still being performed within spatially bound modelling domains in 2-D or 3-D (e.g.Quinquis et al., 2011;Jadamec and Billen, 2012).The limited spatial domain particularly requires prescribing boundary conditions on the vertical sidewalls of the domain.These conditions are an important influence on the development of the model interior (Quinquis et al., 2011;Duretz et al., 2011Duretz et al., , 2012;;Ueda et al., 2012).The usual attempt to reduce possible sidewall influence is by moving these far away from where subduction occurs by using a sufficiently large aspect ratio of model length to depth (e.g.Cizkova et al., 2012).Boundary conditions on the vertical sidewalls can be no-slip (no flow at the boundary), free slip (impermeable; no flow through the boundary and zero traction in tangential direction), or open to some particular form of through-flow.
Free slip is the most commonly used boundary condition while open boundaries have been mostly limited to completely prescribed in-and outflow (e.g.van Hunen et al., 2000;Baes et al., 2010;Quinquis et al., 2011), or periodic Published by Copernicus Publications on behalf of the European Geosciences Union.
conditions requiring that the through-flow at one side is the mirror image of through-flow on the other (e.g.Enns et al., 2004;Capitanio et al., 2010).
Open boundaries for which the horizontal in-and outflow are defined by a fully internally developed flow, have hardly been used and are the main topic of the present paper.Such open boundaries basically prescribe a hydrostatic pressure condition on the boundary preventing the model to collapse while horizontal in and outflow is free, in the sense that it is driven by the internal dynamics and the usual condition of incompressible flow.Among the range of boundary conditions used, open boundaries may fit best to real-mantle flow conditions surrounding subduction zones.We know of only one example (Quinteros et al., 2010) of Eulerian modelling with non-periodic open boundaries.
Our aim in this paper is to investigate the benefits of using open boundaries as compared to using closed (free slip) conditions at the sidewalls of a two-dimensional (2-D) model domain.We focus on modelling self-consistent subduction driven by internal buoyancy and boundary stress conditions only, i.e. no kinematics are prescribed, in the presence of an overriding (oceanic) plate.Our focus will be on the effects of boundary conditions and model aspect ratios on subduction and mantle evolution.As our results show strong differences between using free slip and open boundaries, we are considering first order aspects only.Our results also show that with open sidewalls increasing the model aspect ratio does not change the overall evolution of subduction and mantle flow.In contrast, closed boundaries keep influencing the evolution of the model even for large model size of 6000 km by 1000 km.The primary reason is that closed sidewalls basically cause return flows from both sides towards the centre of the model, which feeds back artificially into the evolving subduction process.We expect this also to hold for 3-D models despite the larger degree of freedom to develop lateral flow.

Model setup
We model self-consistent, internally driven, lithosphere subduction in the presence of an overriding plate in a 2dimensional Cartesian geometry.Our main focus is on how distinctly different boundary conditions on the two sidewalls, open versus impermeable boundaries, affect subduction evolution.We will also investigate whether increasing the aspect ratio of the model domain from 3:1 (3000 km × 1000 km) to 6:1 (6000 km × 1000 km) is of influence, particularly, in reducing any observed adverse effect of boundary conditions.
The boundary condition for the top and the bottom of the box is free slip (impermeable).The surface condition will not allow for modelling topography but, as discussed later in Sect.2.3, we impose a low-viscosity top layer of crust, which will partly compensate for the impermeability condition by allowing for lithosphere bending prior to subduction.For the right and left sides, different types of boundary conditions were implemented: open boundaries or free-slip boundaries.Open boundaries are implemented by constraining zero tangential velocity on the boundary and by imposing a lithostatic pressure condition for the normal stress on the boundary: σ n = P lith .
This allows for horizontal in-and outflow purely driven by the internal dynamics of the model.The pressure condition prevents that the model collapses.As discussed in the introduction, open boundaries are hardly used in subduction modelling, but provide a more natural simulation of the mantle outside the model domain than the more common free slip, impermeable boundary conditions.The free-slip condition prevents material transport through the boundary and forces the flow parallel to the boundary.

Governing equations
We adopt the Boussinesq approximation comprising three coupled equations, namely mass conservation of an incompressible viscous fluid, the Stokes equation describing force balance, and the heat equation, which here only takes into acount heat diffusion and heat advection: (for explanation of symbols see Table 1).This system of equations is solved numerically using the finite element modelling package SEPRAN (Segal and Praagman, 2005).The mesh element size varies from 1.5 km in the trench region to 20 km at the bottom of the model.Advection of the low viscosity material defining the crust and wedge is performed with a Lagrangian tracer technique where material properties are defined on tracers that are advected with the flow.Tracers are distributed initially only over the top 200 km of our domain where we use them to define rheological properties for a low viscosity top layer and wedge.

Rheological model
A composite rheology is used, which comprises dislocation and diffusion creep and a viscosity maximum η max (Fig. 1).The effective viscosity η eff is determined as , 3, 313-326, 2012 www.solid-earth.net/3/313/2012/with η max = 10 24 Pas limiting the effective viscosity in the coldest parts of the lithosphere, with the viscosity due to diffusion creep, and dislocation(power law) creep where ˙ is the second invariant of the strain-rate tensor, A diff , A dis are diffusion and dislocation creep viscosity prefactors, µ is the shear modulus, b is Burgers vector, d is the grain size, m is the grain size exponent, V diff,dis and E diff,dis are activation volume and activation energy for diffusion and dislocation creep, respectively, P is the lithostatic pressure and T -temperature (Table 1).Parameters are taken for wet olivine (Karato et al., 2001), values are given in Table 1.Activation volumes, energy and grain size were chosen to fit seismic studies and postglacial rebound estimations of upper mantle and astenosphere viscosities (Kaufmann, 2000;Burgmann and Dresen, 2008;Simmons et al., 2006).

First order phase changes
Our models include the two major phase transitions at approximately 410 km and 660 km depth.The values of the Clapeyron slope and density contrast are given in Table 1.These parameters are chosen following Billen (2010).The 410 km phase change contributes to the buoyancy force and increases slab pull.The phase change at 660 km has a positive buoyancy effect on cool material resisting slab penetration to the deeper mantle.We ignore thermal effects associated with the phase changes.Phase transitions are parameterized in the model with the phase-transition function: where , w is the half-width of kth transition zone, z tr and T tr are the reference depth and temperature of the phase transition, respectively, γ k is the Clapeyrone slope, T is the temperature (Cristensen and Yuen, 1985;van Hunen, 2001).

The starting configuration leading to the initial buoyancy field
To enable internally-driven subduction, first an initial buoyancy distribution is created by kinematically forced subduction (10 cm yr −1 ) along an arc-shaped fault extending from 0 to 300 km depth.The lithosphere temperature distribution prior to subduction is determined from the equation of cooling of a semi-infinite half-space (Turcotte and Schubert, 2002) for a lithosphere age of 100 My.Boundary temperature conditions are T = 273 K at the surface, T = 2000 K at the bottom.On the side boundaries we prescribe a stationary temperature profile during the subduction process.This also defines the thermal structure of the overriding plate.
We implemented a 10 km thick weak layer (10 19 Pas) on the top of the subducted plate (Han and Gurnis, 1999;Manea and Gurnis, 2007;Cizkova et al., 2007;Behounkova and Cizkova, 2008;Babeyko and Sobolev, 2008;Quinteros et al., 2010).This mimics a subduction channel and proves sufficient for initiating and maintaining subduction in our modelling.An accretionary wedge of weak crust is formed above the subduction zone and has the advantage to prevent artificial rheological coupling between the subducting and overriding plates.The presence/absence of the accretionary wedge only leads to small differences in dip angle and stress field of the trench zone and does not affect the overall subducting slab evolution for the boundary conditions we consider in our models, we conclude from various tests.The focus in our paper is on the large-scale evolution of the subduction system linked to various boundary conditions and aspect ratios of the model domain rather than on the detailed evolution of the plate boundary region.
Figure 1a shows the initial rheology field for the model with open boundaries.In the oldest part of the slab, both dislocation and diffusion creep give high viscosity values, which are limited here by η max = 10 24 Pas.In the asthenosphere viscosity decreases to values of 10 19 Pas, below which it increases to 0.5 × 10 21 Pas in the transition zone and to 10 22 Pas at the top of the lower mantle.The tip of the slab in this initial configuration shows thickening due to mantle resistance, which is also visible in the τ 22 component of the stress distribution (Fig. 2).The starting configuration, rheology and flow field, using free-slip sidewalls is illustrated in Fig. 1b.The two types of boundary conditions, closed or open, lead to a different internal flow field and velocity gradient field on which viscosity depends.Therefore, the starting configuration, particularly the viscosity field, depends on the boundary conditions used and is different when using open sidewalls or free-slip sidewalls.
The dominant deformation mechanisms acting in the initial model are shown in Fig. 3. Diffusion creep(red) is dominant below the asthenosphere and away from the slab, where dislocation creep is active.Small red regions beneath the overriding plate correspond to low strain rate regions.In the core of the slab and overriding plate, the viscosity is limited

Results of numerical modelling
We first focus on models with open versus closed boundary conditions and on model domains with different aspect ratios.This concerns end-member models driven by slab buoyancy only.Subsequently, we incorporate far-field effects additionally constraining the motion of the upper plate.These are imposed by means of a normal-stress boundary condition (τ 11 ) acting on the strong lithosphere at the sidewall.This condition allows investigating stationary or advancing subduction using open boundaries of which we will show several results.
After the initial 3.5 My of kinematically forced subduction has provided an initial buoyancy configuration as discussed in Sect.2.5, the forcing is removed and the internal dynamics take over in driving subduction.We label the models using "O" to denote open, "C" for closed, "R" for "spreading ridge", and "3" and "6" to denote aspect ratios 3:1, and 6:1, respectively.

Open versus closed vertical boundaries and aspect ratio of 3:1
For the models with aspect ratio 3:1, Fig. 4 shows the velocity field and rheology structure of four models (columns) at 2 My intervals (rows) and with different boundary conditions on the left and right sidewalls.Figure 4a depicts the evolution of the model OO3 with two open side boundaries.During the initial stage, about 3.5-4 My from the beginning of the subduction, we observe a strong horizontal flow associated with the subducting plate.This flow pattern bends into the subduction zone following the slab.During this early stage, the upper plate does not move appreciably.However, in the next 2 My slab rollback starts and forces an overall left directed horizontal flow of the upper plate and its underlying asthenosphere.The subducting plate is still being pulled into the subduction zone, but with decreasing speed.Around 5.5 My the slab reaches the base of the upper mantle and further subduction meets resistance due to slab interaction with the 660 km phase transition and with the increased mantle viscosities of the lower mantle.This is followed by the onset of increased slab retreat.The overriding plate attains velocities of 10-15 cm yr −1 , while the advance velocity of the subducting plate drops down to small values (2-3 cm yr −1 ).During the subsequent subduction evolution, we observe strong horizontal left-directed laminar flow concentrated in the asthenosphere and characterized by low viscosity (red colours) in Fig. 4a.This development proves to be characteristic for models with open boundaries and driven by slab buoyancy only.
Figure 4b illustrates the evolution of a model with a closed left and an open right boundary (CO 3 ).The closed left boundary effectively fixes the upper plate to the side.During the whole evolution of the subduction process, we observe a low velocity magnitude (more than 10 times lower than for model OO3).As a result also the viscosity structure is different from model OO3.The minimum viscosities for the asthenosphere, corresponding to lower strain rate, are higher (around 1 order of magnitude), and the rheological width of the asthenosphere is reduced considerably.After the slab reaches the top of the lower mantle, we observe no more trench retreat and the velocity magnitude drops to almost zero.This occurs because the left part of the domain effectively forms a closed volume at this stage in which the slab is blocking flow to the right part of the domain.This is a pure 2-D effect.In 3-D modelling flow is developing at the (open) lateral edge of the slab, which allows material to move away sideways from below the slab (e.g.OzBench et al., 2008;Stegman et al., 2010;Funiciello et al., 2003;Piromallo et al., 2006).
Figure 4c shows the evolution of a model with closed boundaries on both the left and right sides.The closed left boundary again fixes the upper plate, but at the closed right boundary a spreading ridge is allowed to develop in the upper right corner of the model (CCR3) (e.g.Enns et al., 2004).The spreading ridge enables the lithosphere plate to separate from the vertical boundary by allowing hot asthenosphere to flow upward.This facilitates more free lateral movement of the upper plate and is initiated by defining a warm weak zone at the boundary in the initial temperature field.For this model we also observe much lower flow magnitudes (1 cm yr −1 after arriving of the slab at the 660 km boundary) compared to OO3 and, correspondingly, a different evolution of the viscosity structure.The slab behaviour and trench rollback are similar to the model CO 3 except one difference: in model CCR3 a gap between the subducted slab and the overriding plate forms.This gap is filled with asthenosphere material as can be seen from Fig. 4c  Lastly, Fig. 4d illustrates the evolution of a model with closed boundaries but with spreading ridges in both upper corners (CRCR3) allowing both lithosphere plates to separate from the sidewall.Although the 3:1 aspect ratio was wide enough to enable the development of large horizontal flow in the model with open boundaries OO3, this proves more complicated for model CRCR3 with closed boundaries in combination with free plates.We observe strong boundary influence on the motion of the plates and the flow field around the slab.This flow is directed upward close to the vertical boundaries and pushes both plates toward the subduction zone.For this model with closed boundaries and free plates, we observe a lower velocity for the overriding plate (5-6 cm yr −1 on average) compared to OO3 model and higher viscosities for the asthenosphere.The rollback process develops slower than in model OO3, but the evolution of the slab and the overriding plate, as well as the evolution of rollback, are more similar to the model OO3 than to models CO 3 and CCR3.No back-arc basin evolves in model CRCR3.The reason is that, while rolling back, the subducting plate is still advancing as can be inferred from the slab reaching larger depths, which reduces the sinking velocity and no gap appears between upper and lower plate.

Comparison of the end-member (OO3 and CRCR3) models
From the velocity vector plots, Fig. 4a and d, a clear difference in the magnitude of flow between the models OO3 and CRCR3 is observed, illustrating the strong impact of the choice of boundary conditions.Despite the smaller rollback rate, we observe higher horizontal velocities of the subducting plate for the closed boundary case CRCR3, as a result of the ridge push from the spreading centres at both sides of the subduction zone.The flow patterns in both models are distinctly different: in model OO3 a strong lateral flow is created by slab rollback pushing material out of the model domain.This sets up a high-magnitude channel flow in the low-viscous asthenosphere (peak velocity 16.5 cm yr −1 ).For the model with closed boundaries, the flow aligns with the slab and splits up into two cells with upward limbs near the side boundaries.Besides the differences in magnitude of the velocity field, we notice also a difference in rheological structure.For model OO3, at the start of the subduction process, the minimum viscosities in the asthenosphere are around 10 19 Pas.For model CRCR3 these are by one order of magnitude larger.These viscosity values decrease with time, but the asthenosphere in model CRCR3 stays narrower than in model OO3.This feature is related to the dominant deformation mechanism in the asthenosphere, which is strain-rate dependent dislocation creep (Fig. 2).Overall, viscosity values in model CRCR3 are larger during the entire subduction process.
The evolution of the subduction angle is also different between the two end member models OO3 and CRCR3: for model OO3 it gradually decreases while for model CRCR3 it increases with time.This is related to the large difference in the speed of trench retreat while the deep part of the slab is not moving backward.In OO3 this changes the average slab dip and subduction angle.In model CRCR3, the slab penetrates deeper into the lower mantle during the subduction process as a result of a smaller rate of trench retreat creating a steeper average subduction angle.

Comparison of models with different aspect ratios
One possible way to reduce the influence of the sidewalls conditions on the evolution of lithosphere subduction is to increase the width of the domain, but at increased cost of computations.To investigate this we increased the domain width to 6000 km doubling the aspect ratio to 6:1. Figure 5 shows results for models with open and closed sidewalls with spreading ridges, labelled OO6 and CRCR6, respectively.The flow fields are illustrated here by plotting the instantaneous streamlines, which show again significant differences between both models.
For model CRCR6, the velocity magnitude is smaller (after 14 My of subduction maximum velocities are 6.5 and 3.5 cm yr −1 for OO6 and CRCR6 respectively) and the lowest viscosities beneath the slab are at least 10 times higher than in model OO6.Similar observations were made for model OO3 and CRCR3.Figure 6 shows vertical profiles of horizontal velocity computed at different distances from the left side of the models OO6 and CRCR6 after 14 My from the beginning of the subduction process.These profiles clearly illustrate the difference between the flow regimes for these two models.Particularly for model OO6, it illustrates the channel flow regime of the asthenosphere as well as plate-like behaviour (van den Berg et al., 1991;Turcotte and Schubert, 2002).Flow velocities in the asthenosphere channel are much higher for the open boundary case, while plate velocities are higher for the closed boundary case with free plates.Apparently the slab rollback process produces a significant pressure gradient that drives laminar flow in the astenosphere.
In contrast, model CRCR3 shows weaker slab rollback associated with only little asthenosphere return flow from the far-field.Flow patterns around the slab in model CRCR6 are different from those in OO6 due to the different retreat velocity and the return flow resulting from the closed sidewalls.In model OO6 with open boundaries streamlines cross the slab and leave the domain, illustrating that the boundary does not obstruct the slab migration.For the closed boundary model CRCR6 the flow tends to follow the slab then deflects downward and forms closed streamlines.In this model the return flow in the upwelling limbs of the convective cells at both sides contributes to the convergence of the two plates, which puts the subduction channel under lateral compression.For models with a thin crustal weak layer this lateral compression may lead to locking between the subducting and overriding plates.We conclude that using a larger aspect ratio does not reduce the differences in subduction evolution and overall flow field between models with open and closed boundaries.

Rollback velocity and the overall magnitude of flow speed
The almost uniform velocities within the subducting and overriding plates allow us to focus on the speed of slab rollback measured from the trench position through time.Figure 7 shows the speed of rollback for the models with aspect ratio 3:1 and 6:1 and different side boundary conditions.For the models OO3 and CRCR3, yellow and green curves respectively, Fig. 7 shows the difference between rollback speed amounting to a factor of nearly 3 after 8 My of subduction.We observe common trends in the development of the rollback speed: (1) until the slab reaches the 410-km phase change, the speed of rollback increases with a factor of 2; (2) next subduction rollback slows down until the slab reaches the 660-km phase change and the top of the lower mantle; (3) Lastly, slab rollback continues while the tip of the slab is hanging in the highly viscous lower mantle, i.e. without an increase in the penetration depth of the slab.
For model OO3, we observe an increase in rollback speed during this last stage of subduction, which is linked to the trench coming closer to the open boundary.This latter effect is not observed in the 1:6 aspect ratio model OO6 (Fig. 7, blue curve).For model OO6, the rollback speed stabilises after the slab tip gets stuck in the high viscosity lower mantle.Model CRCR6 (Fig. 7, red curve) with free plates has a particular evolution from initially no rollback to strongly increasing speeds peaking around 8 My.This evolution is dominated by the detaching of the overriding plate from the right boundary while low-viscosity mantle material starts filling the gap. Figure 7 demonstrates in a different way the large effects of boundary conditions, open versus closed, which cannot be reduced using larger aspect ratio.
Another large difference between the results obtained for different aspect ratios of the modelling domain concerns the overall magnitude of flow speed of the lithosphere plates, subduction speed and mantle flow.Figure 7 shows that for the 3:1 model with open boundaries the average rollback speed is roughly 15 cm yr −1 , whereas for the model with aspect ratio 6:1 it is around 5 cm yr −1 .A similar reduction in flow speed characterizes the mantle flow.In all aspect-ratio cases the initial driving slab buoyancy is the same.These flow speed reductions can be tied to the longer length of the lithosphere plates in the 6:1 models causing a longer sublithospheric frictional shear zone where a significant part of the mechanical energy of the system is dissipated.
For models with open boundaries, we can, in an approximate way, compensate for the reduced (increased) effect of viscous dissipation in the 3:1 (6:1) models by a scaling relation applied to the velocity field and which accounts for the effect of additional viscous dissipation in an extended computational domain.To determine this velocity scaling, we investigated aspect ratios of 3:1, 4:1, 5:1 and 6:1, where an aspect ratio of 6:1 would correspond to the ocean halfwidth of 3000 km of the subducting plate.The scaling procedure is explained in more details in Appendix A. The scaling procedure compensates the overall velocity field for the effect of bottom-side traction of that part of the plate located outside the model domain, which allows to account for bottom-side traction of the whole plate from the spreading ridge to the trench.In combination with intraplate stress imposed at the sidewalls, it facilitates modelling of a natural subduction process with correct plate length and ridge push within a smaller model domain.
This scaling procedure has been applied during the computation of the OO3 model to approximate flow speed results as observed for the 6:1 aspect ratio.Figure 8 shows scaled rollback velocities of the OO3 model together with original (unscaled) velocities for both the OO3 and OO6 models.These results illustrate the feasibility of approximate upscaling the numerical results for a larger domain.For the models with closed boundaries, we cannot apply velocity scaling due to lateral variations of the flow close to the side boundaries.For models with closed boundaries not only the velocity magnitude but also the flow pattern is changing with decreasing/increasing domain size.In this case comparison of the subduction dynamics in models with different aspect ratio is not meaningful.

Constraining the motion of the lithosphere plates
In the open boundary models OO3 and OO6, the motion of the subducting and overriding plate are entirely controlled by the buoyancy of the subducting slab.We invariably observe (relatively fast) slab rollback in these models, whereas on Earth not all subduction zones show strong rollback and advancing trenches are also proposed (e.g.Funiciello et al., 2008;Schellart et al., 2008).OO3 and OO6 are in fact endmember models as on Earth the global coupling between plates may impose far-field control on the velocity of both overriding and subducting plate.We devised a number of experiments to investigate the combination of far-field control and local slab buoyancy on the evolution of subduction us- To avoid prescribing plate velocities, we can devise a more general implementation of the far-field control by imposing an intraplate stress as a normal traction on the open boundary from the surface down to the base of the lithosphere.We applied this to the upper plate only.Some examples are presented here demonstrating that open boundaries in combination with intraplate stress constraints can lead to strongly reduced slab rollback, (temporary) stationary subduction, or even advancing trenches as compared to the end-member "free" in/outflow models OO3 or OO6.
We varied the applied intraplate stress within reasonable limits for subduction zones (up to 52 MPa, Lithgow-Bertelloni and Guynn, 2004).The results of the subduction models with different intraplate stress values applied to the upper plate are presented in Fig. 10 where a pull is exerted of 18 MPa (10A), 36 MPa (10B), and 52 MPa (10C) on the upper plate.With increasing value of the intraplate stress we generally observe decreasing trench retreat.In more detail Fig. 10b shows an initial stage of stationary subduction, while Fig. 10c exhibits an initial phase of trench advance (to the right).These initial stages differ when as a result of a relatively short slab, the slab pull is still small and the pull on the upper plate is able to force stationary subduction or even slab advance.When the slab touches the 660-km boundary, there is a short episode of trench-stationary subduction, and after a few My trench retreat develops in the model (Funiciello et al., 2003).These results demonstrate that open boundaries are not restrictive on modelling rollback, stationary subduction or trench advance.In addition they again demonstrate the strong dependence of the evolution of trench motion and slab morphology on boundary conditions.

Discussion and conclusions
In this paper we set out to investigate the merits of using open sidewalls in 2-D modelling of subduction evolution as opposed to the more common impermeable free-slip condition.The particular implemented condition is to maintain lithostatic pressure at the boundaries while flow perpendicular to the boundary is free.The internal buoyancy in combination with normal stress conditions (pull/push) on the cross sectional area of the two lithosphere plates, is driving the flow.The absence of kinematic boundary conditions leads to a fully dynamic, self-consistent evolution of the internal dynamics of the model.Simulating a weak upper crust (10 km thick, 10 19 Pas) allowed for modelling continuous subduction without prescribing a particularly subduction channel geometry.We observed that changing the aspect ratio for models with open boundaries did not change the general flow patterns and subduction evolution, except for a general drop in flow speed amplitude for which we derived an approximate scaling procedure based on energy dissipation.
The modelling results obtained with the usual free-slip condition (no horizontal flow) at the sidewalls are in strong contrast with the results obtained for open boundaries.We observed in all experiments a strongly deviating subduction evolution due to the unavoidable influence of return flows induced by the free-slip boundaries and an order of magnitude difference in flow-dependent effective viscosity.If the lithosphere plates are allowed to move away from the boundary, which we implemented to relax the free-slip condition on the motion of the plates, upward return flows at either boundary have an important adverse effect on subduction evolution, by forcing additional convergence, and model evolution in general.These effects could not be sufficiently reduced by taking a larger 6:1 aspect ratio (simulating the half-size of modern oceans), whereas for open boundaries it was found that the flow pattern and subduction evolution was basically independent of the aspect ratio.Also for closed sidewalls a general drop in flow speed was observed for larger aspect ratio models.However, the free-slip boundaries prevented application of a useful scaling procedure.
While the width of the model domain was varied we kept the depth of the model constant.In our rollback models, the slab is draping on the 660-boundary as a result of increasing viscosity in combination with a decrease in negative buoyancy resulting from the 660 km phase transition.In this model scenario, the effect of the bottom boundary will not be very strong.For the investigation of slab behaviour in the lower mantle a model with much larger depth is required.But, for such models we still expect that the open boundaries would allow using a smaller lateral extent of the domain (2000-3000 km) in comparison to free-slip models to reduce the effect of the sidewall boundary conditions.
Although our experiments are in 2-D, we do expect that using open boundaries in 3-D modelling of subduction evolution (e.g.Jadamec and Billen, 2010;Piromallo et al., 2006;van Hunen et al., 2011) may proof beneficial.In 3-D, flow patterns have larger degrees of freedom as toroidal flow can be excited, e.g.around slab edges, and perhaps remote freeslip boundaries may suffice, but as no material is allowed to leave or enter the model unfavourable effects of free-slip sidewalls cannot be excluded.
Our primary conclusion is that open boundaries lead to the most natural boundary conditions for modelling realistic subduction evolution in 2-D and, we expect, 3-D as they avoid adverse effects of impenetrable walls, which are not present in Earth's mantle either.Open boundaries can be combined with plate push or pull conditions, or kinematic conditions on the plates at the sidewalls, to simulate the far-field control of global plate tectonics.As also noted by Capitanio et al. (2010), andvan Dinther et al. (2010), the tectonic style of subduction can be strongly controlled by far-field effects on both upper and lower plate.We notice a strong interaction between boundary conditions, internal flow field, and the viscosity field when using non-linear strain-rate dependent rheology.For closed boundaries this feeds back into an effective viscosity of one order of magnitude larger in the asthenosphere, as compared to the asthenosphere viscosity in open-boundary models.
Open boundary flow, if any, is restricted to be perpendicular to the sidewall.This condition, as any type of boundary condition, poses constraints on the internal flow.But the advantage we experience from using open boundaries is that we can deal with lateral in-and outflow of the lithosphere in a more general way (in 2-D and current 3-D models), e.g. by prescribing intraplate stress, or even just by imposing kinematic constraints on the side walls, while lateral flow below the lithosphere can fully develop, in the sense that it is not influenced by (possibly remote) impermeable vertical boundaries.
Other advantages are the independence of the aspect ratio of the model domain, which allows for smaller models with increased resolution for modelling detail.An approximate scaling procedure can be used to tune the overall flow speed amplitude to levels consistent with the mantle outside the model domain as far as buoyancy inside the model would drive motions outside the model.
2. Estimate the corresponding power, dissipated in the virtual extensions of the model, the regions labeled 1 and 2 in Fig. A2, This is done by uniform lateral extrapolation of the dissipation function profiles at the left and right hand side boundaries of the interior domain.
3. A velocity scaling factor is defined as out + P (j ) in ) P (1) in (A4) 4. Apply the scaling factor to the velocity field, U (j +1)  = f j U (j ) . (A5) This procedure is repeated until convergence | f j +1 − f j |≤ 10 −6 , which is typically obtained within a few iterations.Despite of the fact that the values for the dissipation in the external regions represent less than 15-20 % of the total dissipation in the model, this iterative procedure leads to a significant reduction in the velocity field magnitude, due to the non-linear rheology.
The above-mentioned procedure was tested on several models with different aspect ratios.The results of the velocity scaling for these models are shown in Fig. A3.
The results of the velocity scaling for these models are shown in Fig. A3a where maximum velocity values in the domain for the model with aspect ratio 6 and scaled models with smaller aspect ratios are presented.The maximum velocities are typically observed in the asthenosphere below the overriding plate.The maximum evolution time is 18 My since, for the model with the smallest aspect ratio, the trench has reached the left boundary at this stage.The timeaveraged difference in maximum velocities is ∼6 % for the 4:1 and 5:1 model and 14 % for 3:1.The maximum difference in velocity magnitude is, 20 % for the 5:1 model 30 % for the 4:1 model and 45 % for the 3:1 model.Without velocity scaling the difference in maximum velocities between 3:1 and 6:1 aspect ratio models after 12 My of evolution reaches approximately 70 % and increases over time as illustrated in Fig. A3b.

Fig. 1 .
Fig. 1.The rheology and flow field after the initial 3.5 My of kinematically forced subduction consistent with (A) -open sidewalls, (B) -closed sidewalls.The colour scale shows 10-logarithm of effective viscosity.White lines corresponds to the approximate position of phase transition zones at 410 km and 660 km.

Fig. 2 .
Fig. 2. Stress component τ 22 for the initial model with open boundaries after 3.5 My of kinematically forced evolution (Fig.1a).Note the down-dip compression in the slab.

Fig. 3 .
Fig. 3. Dominant deformation mechanism in the initial model with open side boundaries of Fig. 1a.Different colours correspond to regions where the individual deformation mechanism are dominant.

Fig. 4 .
Fig. 4. Evolution of the subduction process for model OO3 with open boundaries, model CO 3 closed left and open right boundary, model CCR3 with closed right and left boundaries with spreading centre on the right boundary and model CRCR3 with closed boundaries.Arrows show the direction and magnitude of flow field.Identical scaling of the velocity vectors applies to all cases.
mimicking the formation of a backarc basin.While in model CO 3 , the open right boundary allows for modest inflow of asthenosphere material with flow speeds comparable to the overlying upper plate, this does not happen in model CCR3.When the slab reaches the bottom of the upper mantle material exchange with the left part of the model is mostly blocked as in model CO 3 .However, the closed right boundary now forces an asthenospheric return flow in response to slab rollback with larger magnitude than the free upper plate can attain.As a consequence back-arc basin opening develops.

Fig. 5 .
Fig. 5. Snapshots of the effective viscosity and flow pattern for 2 subduction models, on the right-model domain with ratio 6:1, closed boundaries (CRCR6), on the left-model domain with ratio 6:1, open boundaries (OO6).

Fig. 6 .
Fig. 6.Velocity profiles for the model with open boundaries, OO6 (solid lines) and model with closed boundaries, CRCR6 (dashed lines).Profiles were taken after 14 My from the beginning of the subduction process, distances from the left side of the model domain indicated in the legend.

Fig. 7 .
Fig. 7. Speed of rollback for the model with open boundaries and model with closed boundaries: model OO3 (yellow line), model CRCR3 (green line), model OO6 (blue line) and model CRCR6 (red line).

Fig. 8 .
Fig. 8. Evolution of the speed of rollback for two models: model OO3 with aspect ratio 3:1 (blue line), and model OO6 with aspect ratio 6:1 (red line) and scaled speed of rollback for the model OO3 with aspect ratio 3:1 (green line)

Fig. 9 .
Fig. 9. Evolution of the subduction process for model OOF3 with open boundaries and fixed overriding plate (on the left) and model OC3 with closed right boundary (on the right).Arrows show the direction and magnitude of flow field.

Fig. 10 .
Fig. 10.The evolution of the models with open boundary conditions with different intraplate stress values on the overriding plate: (A) model with intraplate stress 18 MPa; (B) model with intraplate stress 36 MPa; (C) model with intraplate stress 52 MPa.Vertical black lines represent the initial trench position.

Fig. A1 .
Fig. A1.Energy dissipation for the model with open boundaries and aspect ratio 3:1 at time = 3.5 My.

Table 1 .
Parameters of the model.