Nonlinear viscoplasticity in ASPECT: benchmarking and applications to subduction

. ASPECT (Advanced Solver for Problems in Earth’s ConvecTion) is a massively parallel ﬁnite element code originally designed for modeling thermal convection in the mantle with a Newtonian rheology. The code is characterized by modern numerical methods, high-performance parallelism and extensibility. This last characteristic is illustrated in this work: we have extended the use of ASPECT from global thermal convection modeling to upper-mantle-scale applications of subduction. Subduction modeling generally requires the tracking of multiple materials with different properties and with non-linear viscous and viscoplastic rheologies. To this end, we implemented a frictional plasticity criterion that is combined with a viscous diffusion and dislocation creep rheology. Because ASPECT uses compositional ﬁelds to represent different materials, all material parameters are made dependent on a user-speciﬁed number of ﬁelds. The goal of this paper is primarily to describe and verify our implementations of complex, multi-material rheology by reproducing the results of four well-known two-dimensional benchmarks: the indentor benchmark, the brick experiment, the sandbox experiment and the slab detachment benchmark


Introduction
Earth is a complex dynamic system that deforms on a wide range of spatial and temporal scales.To obtain realistic predictions of this system from numerical simulations, it is key to capture the relevant aspects of this deformation behavior.
Here we are concerned with the longer geological timescales of the subduction of lithospheric plates into the mantle.On such timescales, rock deformation is mostly nonelastic and characterized by unrecoverable solid-state creep and brittle-plastic failure (Ranalli, 1995;Karato, 2008;Burov, 2011).Strain-rate-dependent viscous deformation through the mechanism of solid-state creep is dominated by linear (Newtonian) diffusion creep and various forms of nonlinear high-and low-temperature dislocation creep (e.g., Ranalli, 1995;Burov, 2011).Plastic yielding occurs when large differential stresses cause rocks to fail beyond the creep regime by local brittle fracture or, at higher temperatures, through ductile homogeneous material flow (Ranalli, 1995;Karato, 2008).
The implementation of plastic yielding into numerical modeling software entails the definition of a yield criterion that the maximum stress must satisfy (Davis and Selvadurai, 2002).Several different plastic yield criteria, such as the Mohr-Coulomb, Drucker-Prager or the Griffith-Murrell criteria (see Braun, 1994;Braun and Beaumont, 1995;Davis and Selvadurai, 2002;Kachanov, 2004, and references therein), are commonly used.These formulations introduce a pressure dependence (frictional plasticity) in the yield criterion.Whereas failure behavior is similar between different rock types and depends primarily on pressure (Burov, 2011), deformation in the viscous creep regime (when stresses are below the plastic yield strength) requires the implementation of rheological descriptions varying with rock type, pressure, temperature, strain rate and other factors such as grain size and water content (Burov, 2011).The implementation of plastic failure and viscous creep complicates solving the governing equations of flow problems due to the nonlinear dependence of the so-called effective viscosity on the solution variables strain rate, pressure and temperature (Gerya, 2010).However, the necessity of using viscoplastic rheologies for simulating natural deformation processes, particularly of the lithosphere, is generally accepted.
To this list we can now add the recent open-source code ASPECT (Advanced Solver for Problems in Earth's Con-vecTion; Kronbichler et al., 2012), which was originally designed for modeling thermal convection in the mantle.AS-PECT is a massively parallel finite element code that is based on state-of-the-art numerical methods, such as highperformance iterative and direct solvers and adaptive mesh refinement, to solve problems of both compressible and incompressible flow.It builds on tried and well-tested libraries such as deal.II (Bangerth et al., 2007;Arndt et al., 2017), Trilinos (Heroux et al., 2005;Heroux and Willenbring, 2012) and p4est (Burstedde et al., 2011) and is under constant development (Bangerth et al., 2017a;Dannberg and Heister, 2016;Rose et al., 2017;Heister et al., 2017).
However, ASPECT originally did not include modeling with multiple nonlinear viscoplastic materials as needed for long-term tectonics modeling, for example.Therefore, we implemented and benchmarked a frictional plasticity (Drucker-Prager) criterion that can be combined with a viscous creep rheology (diffusion, dislocation or composite creep) for any number of materials, allowing for fully thermomechanically coupled viscoplastic flow, on which we here report.There are two papers that use ASPECT that employ a simpler, one-material viscoplastic rheology for planetary convection (Tosi et al., 2015;Zhang and O'Neill, 2016).
Here we focus on benchmarking our implementations in light of lithospheric deformation and these implementations as well as example model setups have or will become part of ASPECT, together with extensive documentation, providing hands-on applications of the code.We show that our viscoplastic rheology description enables the extension of applications beyond thermal mantle convection to detailed lithospheric subduction modeling.
We first present the algorithms underpinning the ASPECT code and our additions pertaining to rheology and compositional fields (Sect.2).We then verify our implementations in Sect. 3 using four benchmarks of increasing complexity: the indentor benchmark (Thieulot et al., 2008;Thieulot, 2014), the brick experiment (Lemiale et al., 2008;Kaus, 2010), the numerical sandbox (Buiter et al., 2006) and the slab detachment benchmark (Schmalholz, 2011;Hillebrand et al., 2014).Finally, in Sect. 4 we present two 3-D subduction applications to showcase the new suite of possibilities made available through our additions and adaptations, and we discuss our overall results in Sect. 5.

Methods
A short summary of the governing equations solved by AS-PECT is given in Sect.2.1 (for more information the reader is referred to Kronbichler et al., 2012;Heister et al., 2017).Section 2.2 lists our specific additions to the code.

Governing equations
ASPECT can solve for both compressible and incompressible flow, but here we focus on the latter, adopting the Boussinesq approximation and assuming an infinite Prandtl number (i.e., inertial term is omitted).Heat production is not incorporated.This results in the following equations of conservation of momentum (Eq.1), mass (Eq.2) and energy (Eq.3): where density ρ = ρ 0 (1 − α(T − T 0 )).Other symbols are explained in Table 1.Artificial diffusivity ν h is used to prevent oscillations due to the advection of the temperature field.It is calculated according to the entropy viscosity method of Guermond et al. (2011), as described in Kronbichler et al. (2012).Similar to the description of temperature, distinct sets of material parameters are represented by compositional fields that are advected with the flow.For each field c i , this formulation introduces an additional advection equation (Eq.4) to the system of Eqs.(1)-(3) described above.As these equations contain no natural diffusion, artificial diffusivity ν h is again introduced to stabilize advection:

Solving the governing equations
ASPECT solves the equations above using the finite element method: the domain is discretized into quadrilateral (in 2-D) or hexahedral (in 3-D) finite elements and the solution (velocity, pressure, temperature and compositional fields) is expanded using Lagrange polynomials as interpolating basis functions.Default settings employ second-order polynomials for velocity, and first-order polynomials for pressure (Q 2 Q 1 elements, e.g., Donea and Huerta, 2003), and secondorder polynomials for temperature and composition.Unless stated otherwise, these default polynomial degrees are used in the following.The linearized Stokes system is solved in a procedure involving the iterative Flexible GMRES (FGM-RES) solver with an inexact right preconditioner.For details on the construction of the preconditioner, see Kronbichler et al. (2012).A cheap Stokes solve option in which the preconditioner employs only one V cycle is available.The number of such FGMRES iterations before switching to the more expensive preconditioner is set to 0 in this paper, unless stated otherwise.The GMRES method with an incomplete LU decomposition preconditioner is used for the temperature and composition systems.Nonlinearities in the rheology are resolved with Picard-type (fixed-point) iterations, iteratively updating the velocity and pressure, strain rate and viscosity (Ismail-Zadeh and Tackley, 2010) until the relative nonlinear residual for iteration i has fallen below a user-set tolerance (default value of 10 −6 ), or the user-specified maximum number of iterations (NIs) is reached.The initial residual ||A(x 0 )x 0 − b|| 2 is computed with zero velocities and a lithostatic pressure profile calculated at the center horizontal coordinate.x contains the velocity and pressure solutions of the previous iteration, b represents the right-hand side of the Stokes equations and A is the Stokes part of the system matrix.

Nonlinear rheologies
The ASPECT code is divided into different modules for boundary conditions, initial conditions, mesh refinement etc.Each module comprises of several plug-ins providing different implementations (e.g., constant vs. space-and timedependent boundary conditions), to which the user can add its own if more functionality is needed.Rheologies are implemented within the so-called Material model module.Plug-ins in this module must provide functions that compute the viscosity, density, thermal conductivity, thermal diffusivity, specific heat and the thermal expansion coefficient at the quadrature points.The solution variables T , P and c i as well as the derived strain rate ˙ (u) and the position are available to compute these material properties.This then provides a straightforward way of implementing nonlinear rheologies, which we have taken advantage of.Deformation of materials on longer timescales is predominantly defined by brittle fracture or viscous creep in terms of diffusion and dislocation creep at relatively low stresses (Karato, 2008).We thus implement three basis rheologies that can be combined into more complex ones: 1. grain boundary or bulk diffusion creep 2. power-law dislocation creep 3. plastic yielding.
Rheologies 1 and 2 can be conveniently formulated with one equation (Karato and Wu, 1993;Karato, 2008): where in case of diffusion creep, n = 1 and m > 0, while for dislocation creep n > 1 and m = 0. See is defined as ˙ e = 1 2 ˙ ij ˙ ij .We simplify Eq. ( 5) by defining prefactor B as 1 Ab m 1 n and add a scaling factor β to easily tune the effective viscosity: The superscript "df" here indicates diffusion creep and "dl" indicates dislocation creep.Plastic yielding (rheology 3) is implemented by locally rescaling the effective viscosity in such a way that the stress does not exceed the yield stress, also known as the viscosity rescaling method (Willett, 1992;Kachanov, 2004).The effective plastic viscosity is thus given by where σ y is the yield value.In our implementation it is defined by the Drucker-Prager criterion (Davis and Selvadurai, 2002): where dilatancy is neglected for simplicity.In case the internal friction angle φ is zero, this criterion reverts back to the von Mises criterion in 2-D.In 2-D it is set to equal the Mohr-Coulomb criterion, while in 3-D it circumscribes the Mohr-Coulomb yield surface (de Souza Neto et al., 2008).Both types of viscous creep act simultaneously (Karato, 2008) under the same deviatoric stress, so the contributions of diffusion µ df eff and dislocation µ dl eff creep to the effective viscosity are harmonically averaged into a composite viscosity (van den Berg et al., 1993): To combine plastic yielding and viscous creep, we assume they are independent (parallel) processes (Karato, 2008), i.e., the mechanism resulting in the lowest effective viscoplastic viscosity is favored: However, for a smoother transition between the different deformation regimes (which should be easier for the numerical scheme to solve), we also experimented with a harmonic average (following Ismail-Zadeh and Tackley, 2010): Because of the strain rate dependence of viscosity and the lack of an initial guess for the strain rate for the first time step, a user-defined initial viscosity µ init is adopted for each compositional field, or an initial uniform strain rate ˙ init is set.We find that the values of µ init and ˙ init can significantly affect the compute time of the first time step.During subsequent time steps, the strain rate of the previous time step is used as an initial guess for the iterative process.
The final effective viscosity µ vp eff is capped by the userdefined minimum viscosity µ min and maximum viscosity µ max to avoid extremely low or high viscosity values due to possible velocity anomalies feeding back into the rheology as well as large viscosity jumps and thus ensure stability of the numerical scheme: We have successfully run the models presented here with overall viscosity contrasts of up to 7 orders of magnitude.Such a range covers the mantle viscosity profiles suggested in most literature, for example as summarized in Cizkova et al. (2012), and we assume that viscosities higher than µ max do not change the behavior significantly.

Multiple compositional fields
Lithospheric geodynamic models often require the specification of materials with different properties, for example a light and weak upper crust versus a denser and stronger lithospheric mantle.To provide the functionality needed for geodynamic modeling, all major material properties of our Material model plug-in depend on any number of fields, as defined by the user (composition-dependent parameters are denoted with an asterisk in Table 1).The use of multiple compositional fields raises the question of how to average their properties (in our case viscosity, specific heat, thermal conductivity, thermal expansivity and density).We have implemented the four averaging schemes commonly referred to in the literature (e.g., Deubelbeiss and Kaus, 2008;Schmeling et al., 2008) for computing the viscosity used in Eq. ( 13) or ( 14): (17) where "nc" is the total number of compositional fields c i in the domain.Note that each field c i is initialized with values on the interval [0,1] and capped values 0 ≤ c i ≤ 1 are used for averaging, as compositional field values may come to slightly exceed this interval over time despite artificial diffusion (Eq.4).The µ vp eff i value is obtained by evaluating Eqs. ( 11) or (12) using the material constants of composition i.The other material properties are arithmetically averaged or, in case the viscosity averaging method is set to Eq. ( 18), averaged using this infinity norm.
The methods above have been shown to affect model results in the context of subduction: Schmeling et al. (2008) showed that the subduction process can be up to 3 times faster between one averaging method and the other, and the effect of mesh resolution on subduction evolution varies per method as well.Unless stated otherwise, we use the infinity norm rule in this paper; for a discussion of this choice, see Appendix A.

Nonlinear rheology benchmarks
To test and verify our implementation of multi-material viscoplastic rheologies, we performed four 2-D experiments: the indentor benchmark, the brick experiment, the numerical sandbox and the slab detachment experiment.The experiments increase in the number of materials and in the complexity of the rheology used, as outlined in Table 2. Consequently, each experiment highlights different parts of the implementation and the functionalities of ASPECT.
All experiments were conducted on an in-house computer consisting of 1 Dell PE-R515 master node and 15 Dell PE-C6145 compute servers made up of 2×4 AMD Opteron 6136 CPUs with QLogic InfiniBand QDR interconnect.ASPECT was compiled using GCC 4.9.2.(Davis and Selvadurai, 2002;Kachanov, 2004;Thieulot et al., 2008).Dark red arrows indicate the prescribed punch velocity v p .The shaded area inside CDE has a resulting velocity of v CDE = v p , while velocities in the lightest shaded areas are . Pressure at point I is P I = σ y (1 + π ) and P ABC = P EFG = σ y .

The indentor benchmark
In the indentor benchmark, a rigid indentor "punches" a rigid-plastic half space.The exact solution to this boundary value problem is given by slip-line field theory (Davis and Selvadurai, 2002;Kachanov, 2004;Thieulot et al., 2008, Appendix B).The analytical solution (Fig. 1) is characterized by three observations: 1.The angles of the shear bands stemming from the edges of the indented area are 45 • .
2. The pressure at the surface in the center of the punch (I) and the pressure in triangles ABC and EFG are P I = σ y (1 + π ) and P ABC = P EFG = σ y , respectively.

The velocity magnitude in areas CDE and ABDC
, respectively.

Model setup
The numerical setup of the instantaneous indentor benchmark comprises a 2-D unit square of purely plastic von Mises material, i.e., its yield value σ y is independent of pressure and remains constant.The material's upper boundary is punched along a distance p by prescribing an inward vertical velocity v p on the otherwise open (stress-free) boundary (see Fig. 2a).The horizontal component of velocity along p is either set to zero or left free to implement the so-called "rough" and "smooth" punch (Lliboutry, 1987;Lee et al., 2005;Thieulot et al., 2008), respectively, where the smooth punch assumes a frictionless contact between the punched medium and the indentor.Model and numerical parameters of the performed indentor experiments are presented in Table 3.

Model results compared to the analytical solution
, where v is the maximum velocity difference at the edges of the punch (see Fig. 1) and x the element size.The surface pressure normalization parameter indicates whether the pressure at the surface is normalized to be zero on average (ASPECT's default) or not.agree with the analytical solution according to criteria 1-3 listed above.Outside the slip lines, viscosity is uniformly high.The low-viscosity shear bands fit the analytical slip lines well (Fig. 3a and f) and stem from the edges of the indentor at a 45 • angle with respect to the top of the medium (Fig. 3b and g).For a rough punch, measurements of pressure in point I and velocity in points K and L deviate from the analytical solution by about 15 and 1 %, respectively, but the block-like behavior of triangle CDE is evident (Fig. 3c).The analytical solution is reproduced with errors < 0.14 % for a smooth punch, but the velocity vectors in Fig. 3h show some horizontal motion of triangle CDE and the velocity field is more diffuse.
When using 200 cheap Stokes iterations for the smooth punch, results are not changed, but wall time is about 1.6 times longer.Using harmonic averaging of the material properties as discussed in Heister et al. (2017) increases the velocity error for the smooth punch to ∼ 1 %, but reduces wall time about 4.7 times.Loosening the linear Stokes solver tolerance by 1 order of magnitude to 10 −8 reduces the wall time of the rough punch by a factor of 1.6, while keeping the velocity error < 1 %.

Discussion
ASPECT successfully reproduces the analytical solution of Prandtl for the rigid-plastic indentor benchmark, a problem with mixed boundary conditions and a nonlinear rigid-plastic rheology with overall viscosity contrasts of 6 orders of magnitude.
The indentor experiment performed also shows a tradeoff between accuracy in pressure and velocity measurements and the rigid-plastic-like behavior of the medium (compare left and right column of Fig. 3).This same dichotomy is seen however in other studies that performed the experiment.Note, for example, the continuous velocity vectors in Huh et al. (1999).Also, Thieulot (2014) shows that the pressure under the punch improves from a 15 % error to a mere 0.5 % error when switching from a rough to a smooth footing.

The brick experiment
As brittle failure in rocks is more appropriately described by pressure-dependent plasticity than by the perfectly plastic deformation (Gerbault et al., 1998) used in the punch problem, our material model plug-in includes frictional plasticity.The brick benchmark has been used to investigate the numerical stability of shear band angles θ and their dependence on the internal angle of friction φ by Kaus (2010) and references therein.Three theoretical relationships have been proposed (Vermeer, 1990): where ψ is the dilation angle (assumed to be zero in our case of incompressibility).

Model setup
In our instantaneous version of the brick benchmark, a viscous-frictional plastic medium with a small viscous inclusion at the bottom boundary (Fig. 4) is either compressed or extended (Lemiale et al., 2008;Kaus, 2010)

Viscous inclusion
Constant density ρ 2700 kg m −3 Linear viscous viscosity µ 10 20 Pa s The background strain rate resulting from the boundary conditions of 2vx Lx = 4×10 −11 40 000 = 10 −15 s −1 can be used as the initial strain rate.The maximum strain rate over all mesh resolutions can be estimated from v x = 4×10 −11 40 000 1024 = 1.024 × 10 −12 s −1 .Yield stress σ y during the first iteration will be minimal at the top of the domain (zero pressure) for the highest friction angle, so that the minimum viscosity over all runs will be 1.7 × 10 19 Pa s.The linear viscous viscosity of the medium is set to 10 25 Pa s, which the maximum viscosity will not exceed.From the variation in friction angle (0 to 30 • ) and lithostatic pressure (0-270 MPa), together with the background strain rate, an estimate for the initial viscosity can be made (1.7-8.5 × 10 22 Pa s).40 × 10 km domain is set to free slip and the top boundary is stress free (material is free to flow in or out).Other domain characteristics and material parameters are given in Table 4.
The angle of internal friction φ (Eq.8) is varied from 0 to 30 • to test the pressure dependency of the implemented plasticity criterion.It is expected that the resultant shear band angle θ varies with the internal friction angle.Shear band angles are automatically computed from the location of the maximum Frobenius norm of the strain rate √ ˙ : ˙ at x = 17.4 km, x = 19.4km, x = 20.6 km and x = 22.6 km (Kaus, 2010).

Model results compared to the theoretical solution
Figure 5 depicts the measured shear band angle versus the supplied internal friction angle for 21 runs in both the compressional and tensional regime.The constant and uniform elemental resolution of the runs varies from 256 × 64 to 1024×256 elements.Lower-resolution runs were performed, but do not resolve the viscous inclusion well (≤ two elements) and are not shown (see instead Fig. 12 of Kaus, 2010).We monitor the residual as a measure of convergence, but the maximum number of nonlinear iterations (NIs; see Sect.2.1.2) is fixed at 1000.The red symbols in Fig. 5 indicate runs for which the residual did not drop below the convergence criterion u = 10 −4 after 1000 iterations, as is evident from the corresponding red lines in Fig. 6.In fact, the higher the internal friction angle, the more iterations are needed to reach a particular residual tolerance, and Fig. 6 also shows that higher internal friction angle runs stall at higher residuals.This coincides with a greater deviation from the theoretical Coulomb solution.Note that in these higher angle runs multiple shear bands are generated and an asymmetry between the right and left shears develops, as shown in Fig. 7.The spurious shear bands can occur mainly at the top of the domain, as several curved pieces forming one shear band or as complete additional shears.Despite these difficulties, there is a clear trend of measured angles verging from Arthur to Coulomb angles for an increasing resolution.For example, the average deviation from the theoretical Coulomb angle decreases from 2.8 to 0.8 • in extension when going from a resolution of 256 × 64 elements to 1024 × 256 elements.
To estimate the effect of adaptive mesh refinement on the shear band angles, we ran additional tests with 3 × 333 nonlinear iterations at increasing refinement levels, with refinement based on gradients in the velocity, viscosity or strain rate and different fractions of cells that are refined.These simulations indicate a maximal variation of 5 • in shear band angle compared to results for a uniform mesh of 512 × 128 elements.Varying the initial viscosity of the viscoplastic medium from µ min to µ max for a uniform mesh of 512×128 elements leads to the same shear band angles for well-behaved residual runs (see black lines in Fig. 6), while for higher internal angle of friction runs, a variation of maximally 3 • is found.

Discussion
Testing the pressure dependency of our plasticity formulation with the brick benchmark, shear band angles were found to increase with internal friction angle φ as expected, almost all falling within the theoretical values of Arthur and Coulomb.Moreover, with increasing mesh resolution, the angles approach the Coulomb theoretical angle θ = 45 ± φ 2 and the variation in error with respect to Coulomb angles decreases.This tendency towards Coulomb angles for higher mesh resolution was also reported by Lemiale et al. (2008), Kaus (2010) and Buiter (2012) because at higher resolution the viscous inclusion is better resolved.Kaus finds that at least 10 to 20 elements are required horizontally within the inclusion to obtain Coulomb angles.This corresponds to our two highest resolutions.Choi and Petersen (2015) recently showed that  Interestingly, internal angles of friction larger than 15-20 • lead to irregular convergence behavior that stalls at higher relative residuals; these runs also show shear band angles further away from the theoretical Coulomb angle and multiple additional shears.These additional shears do often coincide with the theoretical angle; see for instance Fig. 7b.Why these latter shears are not dominant (highest strain rate) would require further investigation.Note that the models of Kaus (2010) also show multiple shear bands and that these spurious bands and stalling of convergence are shown to be the result of the dynamic pressure dependence of the Drucker-Prager yield criterion by Spiegelman et al. (2016).
However, we have shown that we consistently obtain shear band angles between Arthur and Coulomb theoretical angles at sufficient resolution and that these angles verge to Coulomb angles with increasing resolution.This happens despite the fact that our implementation is relatively basic: it does not include softening of cohesion or of the internal angle of friction as in Kaus (2010) and Buiter (2012), nor does it have a sophisticated guess of the initial stress state (Lemiale et al., 2008) or an incremental build-up of the prescribed boundary velocity (Kaus, 2010).

The sandbox extension experiment
Buiter et al. ( 2006) compared numerical and analog models of shortening and extension.We reproduce the numerical sandbox extension experiment, which was originally run with six different numerical codes and compared to the analog results described in Schreurs et al. (2006).This time-dependent experiment has previously been repeated by Thieulot (2011) and -in a symmetrical version -by Gerya et al. (2013).

Model setup
The analog sandbox (Buiter et al., 2006) consists of a basal layer of weak, viscous silicone overlain by brittle sand.The sand is extended by the movement of the right vertical wall and the connected basal plate extending from the wall to, initially, the center of the domain.We model this setup with three compositional fields (Fig. 8): (1) a viscous basal layer, (2) an overlying Drucker-Prager dynamic pressuredependent plastic sand layer and (3) a low-viscous sticky-air layer (Crameri et al., 2012) on top.Extension is driven by a prescribed horizontal velocity on the right vertical boundary and the right half of the lower boundary, mimicking the effect  should not exceed the maximum strain rate estimate vx x = 6.94×10 −6 0.00039 = 1.78 × 10 −2 s −1 .The maximum strain rate can be used together with the minimum yield strength (zero pressure) to compute the minimum viscosity of the sand of about 227 Pa s.
of the moving basal sheet.Basal friction is not taken into account.This approach is appropriate since Buiter et al. (2006) have shown that the nature of the basal contact is less important than the interaction of the velocity discontinuity (VD in Fig. 8) and the silicone.The rest of the bottom boundary has zero velocity (without smoothing of the velocity discontinuity) and this no-slip area increases as the velocity discontinuity moves to the right of the domain.The left boundary is free slip and the top boundary is open.Adaptive mesh refinement (AMR) is applied based on the effective strain rate field or on viscosity and density to obtain a maximum local refinement of 0.39 × 0.39 mm along the shear bands (Fig. 10).The model is run until 2 cm of extension has occurred, equalling 2880 s of model time.Material properties and other model parameters are listed in Table 5.

Model results
The results for 1 and 2 cm of extension of the sandbox are presented in Fig. 9.The evolution of the model closely resembles what was found by Buiter et al. (2006).The initially symmetric system forms two conjugate shear zones stemming from the velocity discontinuity imposed by the velocity boundary conditions.With ongoing extension, the silicone layer distributes deformation and the shear bands spread to the edges of the layer.After 1 and 2 cm of extension, the angle of the shear bands left and right of the velocity discontinuity (see Fig. 9d) measure left ∼ 51, right ∼ 62 and left ∼ 54, right ∼ 60 • .With time the system becomes more and more asymmetric (compare left and right column of Fig. 9).The left side of the domain is at rest, while the outer right footwall moves at the prescribed velocity and we observe that the sticky air properly accommodates the movement of the sand.
The viscosity field (Fig. 9h) is very irregular and displays sharp gradients up to 7 orders of magnitude.Figure 10 demonstrates viscosity-and density-based AMR: refinement is localized in the low-viscosity shear bands, following the evolution of deformation.Through AMR, the total (velocity, pressure, temperature, composition) number of degrees of freedom (DOFs) is limited to on average ∼ 475 000 DOFs, instead of the ∼ 1 383 000 DOFs (∼ 527 000 velocity DOFs) for a uniform resolution, thus decreasing the required computational resources by half (for the same number of cores).

Discussion
The evolution of the numerical sandbox model -a model with AMR, high viscosity contrasts, large deformation and complex boundary conditions -compares well with those shown in Buiter et al. (2006) and Thieulot (2011).Although the shear band angles to the right of the velocity discontinuity of 62 and 60 • after 1 and 2 cm of deformation fall just outside the ranges found by Buiter et al. (2006) and Thieulot (2011) of 45-55 and 45-53 • , respectively, they lie within the theoretical Arthur-Coulomb angles of 54-63 • for a friction angle of 36 • (Vermeer, 1990; see also Sect.3.2).Even when considering that the codes in Buiter et al. (2006) add strain softening by decreasing the friction angle from 36 to 31 • , for which the range of Arthur-Coulomb angles would be 52.75-60.50• , their measured angles lie mostly below the Arthur angle.In the previous section, we demonstrated that our shear band angles fall within the Arthur-Coulomb range and verge towards Coulomb angles with increasing mesh resolution.Differences in the measured angle of shear bands are therefore not surprising.A lower-resolution run (maximum of 0.78 mm resolution, not shown) results in angles to the right of the discontinuity of 60 and 55 • .
Similar to Buiter et al. (2006), Thieulot (2011) and Buiter (2012), we observe an increase in the number of shear bands with resolution as well as a decrease in their width.As ex-plained by Spiegelman et al. (2016), this lack of internal length scale is caused by the singularities in strain rate and pressure deriving from the model setup (e.g., sharp corners of the silicon layer and the discontinuous velocity boundary condition) that are resolved better at higher resolutions, thereby decreasing shear band width.

The slab detachment benchmark
Slab detachment, or break-off, in the final stage of subduction is often invoked to explain geophysical and geological observations such as tomographic images of slab remnants and exhumed ultra-high-pressure rocks (see for example Wortel and Spakman (2000) and references therein).Due to the increased interest in the process of slab tearing, it has recently been the subject of several numerical modeling studies (e.g., Gerya et al., 2004;Andrews and Billen, 2009;Burkett andBillen, 2009, 2010;van Hunen and Allen, 2011;Duretz et al., 2012Duretz et al., , 2014)).Numerical modeling of slab detachment is computationally challenging due to the mesh-resolution dependency of the strain (van Hunen and Allen, 2011) and the gradual decrease in monitoring particle density, the gradual overlapping of level sets (Hillebrand et al., 2014) or the gradual thinning of compositional fields (as is the case here) in the detachment area.Here we test ASPECT with the slab detachment model of Schmalholz (2011), which considers a simplified geometry of detachment by viscous necking of a vertical lithospheric slab of nonlinear rheology in a linearly or nonlinearly viscous mantle.It has been extended to 3-D by von Tscharner et al. (2014).

Model setup
The 2-D detachment model geometry is outlined in Fig. 11; both the lithosphere and the mantle are represented by a compositional field.Schmalholz (2011) prescribes a nonlinear viscosity in the subducting lithosphere given by where η 0 = 4.75 × 10 11 Pa s 1 n and n = 4. Conversion to Eq. ( 5) results in the parameters listed in Table 6.Here we only reproduce the constant mantle viscosity case of Schmalholz (2011), with µ mantle = 10 21 Pa s.We use a set of 20 tracers placed on the outline of the slab to track the width and depth of necking.The necking width and depth, as well as time, are normalized by the initial slab width of 80 km and the characteristic time t c = 7.1158 × 10 14 s, respectively (Schmalholz, 2011).

Results and comparison
The evolution of the detachment model is shown in Fig. 12: after about 20 Myr, the slab is fully necked and detached.Although a thin line of lithosphere composition is still visible, its value is less than 0.5 and thus the infinite norm igwww.solid-earth.net/9/267/2018/Solid Earth, 9, 267-294, 2018  nores this contribution to the viscosity, allowing for full detachment.Upon detachment, slab pull is removed and thus viscosity reaches µ max throughout the remaining lithosphere (Fig. 12c and d).
From comparison of the red and black lines in Fig. 13, it can be seen that the width of the necked zone through time agrees very well with the results of Schmalholz (2011).Only after t = 0.8, do our results start to deviate due to our use of tracers for measuring the necking width.The two other lines illustrate the effect of using harmonic averaging for the compositional fields' contribution to viscosity (blue line) or for averaging the viscosity over the elements (yellow line): although these averaging methods reduce the wall time by a factor of 3 and 2, respectively, they also result in much faster necking.This agrees with the findings in Appendix A.

Discussion
The first three benchmarks focused on plastic rheologies.The detachment benchmark involves modeling of a highly nonlinear, power-law viscosity, which is often used in subduction modeling.Our observed model evolution compares well with that of Schmalholz (2011) and other codes (Hillebrand et al., 2014).It also demonstrates the effective splitting of a compositional field into two bodies and how the rheology interacts with the compositional fields through localization of deformation at the slab hinge.It should be noted that the particular geometry of the slab with its sharp corners results in a mesh dependence of the solution.Differences in model evolution can also arise from the particular viscosity and material averaging method applied.

Viscoplastic subduction models
Lastly, we consider a geodynamical application of the viscoplastic rheology implemented: the spatiotemporal evolution of 3-D subduction.So far, no ASPECT applications to 2-D or 3-D regional subduction have been published.To demonstrate ASPECT's promise in this field, here we present 3-D models of free-plate intraoceanic subduction.Recent 3-D subduction models have been applied in the study of along-strike effects such as oblique convergence (Malatesta et al., 2013), toroidal flow (Schellart and Moresi, 2013), varying lithospheric structure (Mason et al., 2010;van Hunen and Allen, 2011;Capitanio and Faccenda, 2012;Duretz et al., 2014), slab width (Stegman et al., 2006(Stegman et al., , 2010;; PECT and Schmalholz (2011).The ASPECT necking width is calculated from the 20 tracer positions.Because the tracers above the necking zone no longer move after detachment, the thus-calculated width stagnates after t 0.9.

Model setup
We discuss two models; the first is an adaptation of the freeplate model of Schellart and Moresi (2013), which considers no temperature effects and features constant viscosities except for a viscoplastic crustal layer.The second model is an extension of the first, where we add a temperature field and a viscosity dependent on temperature, pressure and strain rate, resulting in a nonlinear viscoplastic thermomechanically coupled system.

Model 1
The first model comprises two free plates: an overriding plate (OP) and a subducting plate (SP) (see Fig. 14).The SP is made up of a crustal layer of non-frictional (von Mises) viscoplastic rheology and a mantle lithospheric layer of constant viscosity (see Table 8 for actual values).The of the SP extends into the mantle for 200 km.The OP is of one composition and has a constant viscosity, as does the surrounding mantle.These four compositions are each represented by a compositional field.

Model 2
The second model augments the first with an adjacent plate (AP) separated from the other plates by a 20 km thick weak zone (WZ).The SP and OP are extended and their thickness is no longer uniform, but based on the temperature field, as if they originate from ridges situated at the left and right vertical domain boundaries.The initial temperature distribution (Fig. 15) in the plates is computed according to the age-based plate cooling model, for a mantle temperature T a of 1593 K and surface temperature T 0 = 293 K (Eq.4.2.24 of Schubert et al., 2001; κ = 10 −6 m 2 s −1 ).Thickness of the plate is defined at the temperature T for which T a −T T a −T 0 = 0.1.The maximum thickness of the plate (for time to infinity) is set to 125 km, based on Schubert et al. (2001).At the trench, both plates have the same thickness as in model 1.The adjacent plate (AP) has a fixed thickness of 100 km for an age of ∼ 60 Myr.From a depth of 125 km, a linear temperature increase is prescribed everywhere to a bottom temperature of 1771 K.The mantle and overriding plate compositions are of nonlinear rheology, with which model 2 differs strongly from model 1.The specific rheological parameters for diffusion and dislocation creep and Drucker-Prager plasticity can be found in Table 8 and the representative viscosity profiles in Fig. 16.

Model 1
Figure 17 depicts the evolution of the subduction system of model 1 over time.The pull of the slab extending into the mantle results in a plastically weakened subduction fault zone through high strain rates.Although this allows for decoupling from the surface, mechanical coupling is strong enough for the OP to move towards the trench.There the OP thickens and a small portion of OP material is entrained with the slab.Within the first ∼ 7 Myr, the velocity of the plates steadily increases about 1 order of magnitude and reaches several centimeters per year.Similar to Schellart and Moresi (2013), poloidal and toroidal flow can be observed Table 8.Three-dimensional subduction model parameters (based on Hirth andKohlstedt, 2003, andRanalli, 1995).The AMR based on composition and viscosity follows the outline of the plates as they move through the mantle (Fig. 17), resulting in a local resolution of roughly 3 km while the mantle is resolved with ∼ 50 km elements.

Model 2
The SP of model 2 steepens for the first 12 Myr (Fig. 18) until full-fledged subduction starts.Subduction is much slower than in model 1: while in model 1 subduction is completed by 35 Myr, rollback of the slab in model 2 only sets in around 50 Myr.Simultaneously with slab rollback, the subduction channel is weakened by asthenospheric inflow into the gap between the OP and SP.Formation of this gap is probably initiated by the fact that the OP does not completely release itself from the lateral boundary (see snapshot at 49 Myr).Even though the slab reaches the bottom boundary around 50 Myr and starts to shallow, it does not lie flat on the 660 km boundary, but continues to hover above it.Also note the halo of increased strain rate around the tip of the slab and the smallscale strain rate features in the mantle compared to model 1.

Discussion
The evolution of model 1 strongly resembles that of Schellart and Moresi (2013), on whose setup the model is based: the slab first sinks freely until it reaches the 660 km bottom boundary and then starts draping while rolling back.The absence of folding of the slab at the 660 km boundary and the differences in timing of the aforementioned events are probably due to the weak, linearly viscous layer of the SP that is left out here (and the lesser extent of the domain in the y direction).This leaves the plate stronger and more resistant to folding than for Schellart and Moresi (2013).Note that our viscoplastic SP of model 2 also does not show draping at the 660 km boundary.In switching from model 1 to this thermomechanically coupled model 2, we found changes to the setup were necessary to avoid subduction of the plate at locations other than the slab tip (i.e., the sides and back of the plate would subduct as well) due to high mantle temperatures.Therefore, we added the adjacent plate and transform fault.By locating the plate ridges at the left and right vertical boundaries, free motion of the plates perpendicular to the trench is still enabled.Mesh resolution here was reduced over time because the refinement strategy chosen focused on the compositional fields, which moved away from the boundaries.This unfortunately increased the coupling of the OP plate to the left boundary, limiting the plate's ability to move.
The subduction evolution of model 1 and 2 in Figs. 17 and 18 clearly differs.Models in the Appendix of Schellart and Moresi (2013) have shown that the addition of adjacent plates in itself does not affect the geometry of the SP over time (although velocities are affected).A test with a uniform viscosity adjacent plate for our model 1 corroborates this.Therefore, the differences derive from the temperature, pressure, strain rate and composition-dependent rheology.In- deed, a snapshot (Fig. 19) of the viscosity field of the SPs and OPs shows that the viscosity of the plates of model 2 is about an order of magnitude higher (cutoff at 10 24 Pa s).As both the SP and the OP of model 2 keep growing at the trench, they also experience more mantle drag than in model 1.Moreover, the model 2 slab tip is surrounded by a higher-viscosity area due to local cooling of the mantle.These rheological differences slowing down the subduction process are also evident from the strain rate in Figs. 17 and 18, showing a much weaker slab and mantle in model 1.A full investigation into the differences between mechanical and thermomechanical viscoplastic models is beyond the scope of this paper.
More elaborate models of subduction should incorporate phase changes and latent heat effects as well as adiabatic and shear heating.This is also possible with ASPECT and we include an example of such a 2-D model in Appendix B.

Discussion
The four benchmarks shown using our viscoplasticity implementations in ASPECT either reproduce the available analytical solution (Sect.3.1) or compare well with theory (Sect.3.2) or the results of other codes (Sect.3.2-3.4).Thus verifying our implementations, the four benchmarks allowed us to set up a 4-D model of oceanic subduction, exemplifying the functionality that our implementations have added.
It should be noted that although the rheology described in this paper is often applied in numerical modeling, more elaborate laws have been proposed.For example, Gerya and Yuen (2007) included dilatant materials and Choi and Petersen (2015) argue that numerical models should incorporate an initially associated plastic flow rule that evolves into a non-associated flow rule with increased slip to assure persistent Coulomb shear band angles while avoiding unlimited dilatation.The inclusion of an intrinsic length scale in the plasticity formulation would work to remove the mesh dependence of the rheology (e.g., strain gradient plasticity; Fleck and Hutchinson, 2001).Another addition would be to include plastic softening or hardening -changes in the yield surface due to the accumulated strain.Considering creep flow laws, improvements could be made by adding Peierls creep, a dislocation mechanism acting at low temperatures and/or high stresses (e.g., in parts of the slab; Karato, 2008;Duretz et al., 2011;Garel et al., 2014).Other authors such as Farrington et al. (2014) and Fourel et al. (2014) have investigated the effect of incorporating elasticity into models of lithosphere subduction, demonstrating that although observables such as dip angle, slab morphology and plate motion are not affected, an elastoviscous or elastoviscoplastic rheology leads to different viscosities in the hinge of the slab.
Incorporating more realistic nonlinear rheologies such as described in this paper creates the necessity for additional nonlinear iterations within a single time step.Also, we have seen that at higher mesh resolutions, more of such iterations are required to converge the solution.This greatly increases model run time and therefore it is important to implement a more efficient nonlinear solving strategy than the Picard iterations currently used by ASPECT.The more sophisticated Newton solver (see for example Popov and Sobolev, 2008;May et al., 2015;Rudi et al., 2015;Kaus et al., 2016;Wilson et al., 2017) will help achieve faster convergence.Convergence behavior has also been suggested to improve from including elasticity (Kaus, 2010), but especially dynamic pressure-dependent plasticity remains difficult to converge for both Picard iterations and Newton solvers (Spiegelman et al., 2016).
Nonlinear rheologies also affect the linear solver by introducing large viscosity gradients.Different strategies to reduce the increased computational time and under-and overshooting of the numerical approximation of the resulting pressure gradient are available in ASPECT.For one, one can reduce the linear tolerance (while making sure the results do not change significantly), as was shown in Sect.3.1.Secondly, a cheap Stokes solver can be employed, although this does not help for each model setup (compare Sect. 3.1 and 3.4).Thirdly, averaging the contributions of the compositional field to the viscosity and other material properties in a specific point reduces the sharpness of viscosity boundaries, making the problem easier to solve, but with the choice of averaging method affecting the model evolution (Sect.3.4 and Appendix A).Lastly, averaging of material properties such as viscosity and density over each element reduces pressure oscillations (Heister et al., 2017), but can also influence the model evolution as was shown in Sect.3.4.

Conclusion and outlook
Numerical modeling of intricate geodynamic processes such as crust and lithosphere deformation and plate subduction encompasses challenges at different levels.For one, the 4-D nature of the subduction process requires state-of-the-art numerical methods to efficiently handle the parallel computations necessary for such large problem sets.Secondly, models should incorporate realistic (non)linear rheologies to mimic nature as close as possible.Thus, the need arises for algorithms that can solve highly nonlinear equations and deal with large viscosity contrasts effectively.Thirdly, far-field effects of mantle flow and plate motion cannot be ignored, and neither can topography building, resulting in a demand for complex boundary conditions such as open boundaries (Chertova et al., 2012) and free surfaces (Kaus et al., 2010;Crameri et al., 2012;Rose et al., 2017).
In this paper, we have shown that the open-source code ASPECT is up to these challenges.Building on its modern, massively parallel numerical methods, here we have outlined our basic additions that enable viscoplastic crust and lithosphere modeling.We then tested and verified the algorithms with four different benchmarks that are well-known in the geodynamic modeling community.Last, we highlighted the possibilities arising from the adaptations with 4-D thermomechanically coupled viscoplastic models of interoceanic subduction, showing that ASPECT is a serious contender in the field of lithospheric subduction modeling.
The continued development of ASPECT based on the needs of its expanding user and developer community ensures ever-growing capacities and possibilities.Important recent additions are a full free surface (Rose et al., 2017), the formation and migration of partial melt (Dannberg and Heister, 2016), active particles (Gassmöller et al., 2016), a discontinuous Galerkin method for advection (He et al., 2017), and a Newton solver (Fraters et al., 2017).Also, the extensive user manual (Bangerth et al., 2017a) accompanying all developments is a great asset for new and current users.In consequence, opportunities for future research are reinforced and a firm foundation is provided for ASPECT in the geodynamics community. ) where β is the compressibility 1 Phase changes (changes of one compositional field into another with different material properties) are implemented by extending our, now compressible, multicomponent viscoplastic material model with a depth-dependent transition function (e.g., Christensen and Yuen, 1985): X represents the fraction of the new phase, P t , T t and z t are the reference pressure, temperature and depth of the transition, respectively, and d the transition half-width in terms of pressure (assuming P = P lith = ρgz).The phase function derivatives in Eq. ( B3) are computed as in the latent heat plug-in of the Material model module.
Open boundary conditions are newly implemented as a plug-in to the Traction boundary conditions module by prescribing the traction as (Chertova et al., 2012) with n the outward normal to the domain boundary.P lith is the lithostatic pressure calculated by numerical integration along the domain boundary of the density (i.e., the temperature and composition) of the previous time step.

B1 Model setup
The compressible subduction model considers oceancontinent subduction in a domain of 1600 km by 1000 km (Fig. B1).The model includes phase changes around 410 and 660 km depth; see Table B1.An inward plate velocity is prescribed on the upper part of the left vertical boundary for the first 13 Myr at 4 cm yr −1 , while prescribed outward mantle velocities compensate for this volume increase.After 13 Myr, material is free to move in or out of the domain through open boundary conditions on the left boundary.On the right boundary, no-slip conditions are switched to a prescribed velocity profile after 8 Myr, with a plate inflow of 2 cm yr −1 .Apart from the subducting plate crust, which is of linear viscosity, all materials are nonlinear viscoplastic.At 200 km, crustal material is transformed to mantle material (as is done, for example, by Androvicova et al., 2013).Initial temperature is based on an adiabatic profile in the mantle and linear profiles in the plates.

B2 Model results
Figure B2 shows the evolution of the compressible subduction model in terms of viscosity.When the left boundary is opened at 13 Myr, just after the slab has reached the 410 km phase boundary, sinking velocities increase to 10.5 cm yr −1 .Upon reaching the 660 km phase boundary (with a 2-fold viscosity increase), the slab slows down to about 5 cm yr −1 .The tip of the slab is impeded by the 660 km boundary and moves at around 1 cm yr −1 along the boundary.The change in flow through the left boundary is depicted in Fig. B3: when the left domain boundary is opened up by prescribing stresses instead of a fixed velocity profile, first of all a downward shift of the transition from in-to outflow is seen.Moreover, the subducting plate velocity initially increases up to 11 cm yr −1 , together with an increase in transition and lower-mantle velocity.When the slab reaches the 660 km phase boundary around 16 Myr, lower-mantle outflow goes up rather uniformly, but in-and outflow above decreases.When the slab tip moves over the 660 km boundary, all flow decreases in magnitude -gradually in the mantle and uniformly in the lower mantle.The phase changes are clearly expressed in the density fields: the density isocontours in Fig. B2 show positive topography in the slab at the 410 km discontinuity, while the 660 km transition in the slab occurs deeper.The deformation of the surface at 18 Myr shows a lowered subducting plate (about 1 km below the initial top boundary), a topographic rise of continental crust above the subducting slab (of maximally 8 km) and an elevated overriding continental plate (about 4 km).

Figure 1 .
Figure 1.Prandtl's analytical solution of a rigid die indenting a rigid-plastic half space(Davis and Selvadurai, 2002;Kachanov, 2004;Thieulot et al., 2008).Dark red arrows indicate the prescribed punch velocity v p .The shaded area inside CDE has a resulting velocity of v CDE = v p , while velocities in the lightest shaded areas are v ABDC = v EDFG =

Figure 3
Figure3shows the model results for a rough (left column) and smooth (right column) punch.The solutions obtained

Figure 2 .
Figure 2. The indentor benchmark model setup: a unit square with free-slip vertical and no-slip lower boundaries.The punch area has a prescribed vertical velocity v p ; the rest of the upper boundary is open.
lists Hill's solution for the smooth punch, other numerical studies do not recover this slip-line geometry in 2-D either.In fact, our Prandtl shear band geometry compares well with results ofGerbault et al. (1998, Fig. 6b, smooth),Huh et al. (1999, Fig. 1,

Figure 3 .
Figure 3.The punch benchmark results after 500 NIs for a rough punch (left column) and a smooth punch (right column).(a, f) Viscosity field with analytical slip lines.(b, g) Strain rate norm ( √ ˙ : ˙ ) with measured shear band angles.(c, h) Velocity magnitude with velocity vectors along the surface of the domain and velocity measurements in points K and L. (d, i) Pressure field.(e, j) Pressure along the surface of the domain (colored line) and analytical solution values π + 1 and 1 (grey lines).Rough punch: P I = 4.7382 and P H = P J = 0.6224.Smooth punch: P I = 4.1415 and P H = P J = 0.9999.

Figure 4 .
Figure 4.The brick benchmark model setup after Kaus (2010): a rectangular domain with a prescribed inward or outward horizontal component of velocity on the vertical boundaries (the vertical component is left free).The upper boundary is open, while the bottom boundary is free slip.A small viscous inclusion of 800 × 400 m is placed at the bottom of the domain.

Figure 5 .
Figure5.Measured shear band angle versus angle of internal friction for models in tension and compression.Resolution runs from 256×64 (light grey line), to 512×128 (grey) to 1024×256 elements (dark grey).All models were run for 1000 NIs; red symbol runs correspond to the red line runs in Fig.6.The black lines represent the theoretical angles of Coulomb (solid), Arthur (dashed) and Roscoe (dotted).We have corrected one of the automated shear band angle measurements manually -that of the 1024 × 256 element extension case with φ = 25 • -because it was computed using two different shear bands.

Figure 6 .
Figure 6.Measured velocity residual ( u = ||U i −U i−i || sup ||U i−1 || sup with U the velocity solution; van den Berg et al., 1993; Kaus, 2010) versus the number of nonlinear iterations for models of extension.Elemental resolution is 512×128 elements.Black lines represent runs with well-behaved convergence.

Figure 7 .
Figure 7. Strain rate norm fields for (a) φ = 0 • and (b) φ = 30 • for a 512×128 elemental resolution.The models were run in extension for 1000 NIs.Black lines indicate the theoretical Coulomb angle θ C .Measured shear band angles θ M are also given.

Figure 8 .
Figure8.The sandbox experiment setup afterBuiter et al. (2006).VD represents the velocity discontinuity moving at the same speed as is prescribed on the right and bottom boundary.Left of the VD, the bottom boundary is no slip.The left vertical boundary is set to free slip, while the top is open.The silicone layer measures 10 × 0.5 cm.
a sandbox scale, as they are used in the model.An estimate for the initial strain rate can be made from the boundary conditions ˙ init ≈ vx Lx = 6.94×10 −6 0.2 = 3.47 × 10 −5 s −1 and

Figure 9 .
Figure 9. Results after ∼ 1 and ∼ 2 cm of extension of the numerical sandbox.(a, e) The three compositional fields.Note the asymmetric depression of the sand surface.(b, f) Frobenius norm of the strain rate √ ˙ : ˙ .(c, g) Total pressure field.(d, h) Viscosity field.

Figure 10 .
Figure10.Numerical grid after ∼ 2 cm of extension of the numerical sandbox.Adaptive mesh refinement and coarsening based on the viscosity and density leads to a minimum resolution of 6.25 × 6.25 mm and a maximum resolution of 0.39 × 0.39 mm.

Figure 11 .Figure 12 .
Figure 11.The detachment benchmark model setup of Schmalholz (2011): a symmetric system of nonlinear viscous lithosphere (Eq.19) with a vertical slab extending into a linear viscous mantle.The top and bottom boundaries are free slip, while the vertical boundaries are no slip.Along the outline of the slab are placed 20 passive tracers to track the necking of the slab.

Figure 13 .
Figure 13.Nondimensional necking width versus time for AS-PECT and Schmalholz (2011).The ASPECT necking width is calculated from the 20 tracer positions.Because the tracers above the necking zone no longer move after detachment, the thus-calculated width stagnates after t 0.9.

Figure 14 .
Figure14.Three-dimensional subduction model setups (also, see Table7).(a) Model 1: a free overriding plate (OP) and a subducting plate (SP) with a trench at x = 2400 km.At start-up, the slab extends 200 km into the mantle (measured vertically) at an angle of 29 • .The 100 km thick SP consists of two compositional layers, the OP of only one composition.All boundaries are free slip, except the no-slip bottom boundary.(b) Model 2: a 58.8 Myr old adjacent plate (AP) is added, separated from the OP and SP by a weak zone (WZ) of 20 km width.The initial temperature distribution of the SP and OP is based on the plate cooling model dependent on age, which increases from the ridge situated at the left, respectively right, vertical domain boundary, up to an age of 16 Myr for the OP and 60 Myr for the SP, resulting in thicknesses of 50 and 100 km at the trench, respectively.An adiabat of 0.25 Kkm −1 is prescribed in the mantle.The bottom temperature is fixed at 1728 K, the top at 293 K.

Figure 15 .
Figure 15.Temperature distribution with depth and distance from the ridge for (a) the OP and (b) the SP.

Figure 16 .
Figure 16.Viscosity profiles at 0 Myr for the OP (at x = 2700 km) and SP (at x = 2200 km) of models 1 and 2. For comparison, the viscosity profiles derived by Cizkova et al. (2012) are included.

Figure 17 .
Figure 17.Model 1 -strain rate Frobenius norm √ ˙ : ˙ over time together with SP and OP isocontours.Also shown is the adaptive mesh following the SP into the mantle.

Figure 18 .
Figure 18.Model 2 -strain rate Frobenius norm √ ˙ : ˙ over time together with SP, SP crust and OP isocontours.Also shown is the adaptive mesh following the SP into the mantle.

Figure 19 .
Figure 19.Viscosity snapshots of (a) model 1 and (b) model 2 at comparable moments in the respective subduction evolutions.

Figure A2 .
Figure A2.Slab tip depth versus model time for four different averaging methods of the contribution of the compositional fields to viscosity.Colors indicate the averaging method, while one color goes from light to dark with local resolution, which varies from 256 × 64 elements to 2048 × 512 elements.Minimum resolution is always 128 × 32 elements.(a) Case 1.The dashed red line model has a resolution varying from 128 × 128 to 2048 × 2048 elements.Shaded areas represent results of Schmeling et al. (2008, Fig. 6).(b) Case 2.

Figure A3 .
Figure A3.Viscosity field for each averaging method at 2048 × 512 element local resolution at similar moments in the subduction evolution.

Figure B1 .
Figure B1.Compressible subduction model setup: subduction of an oceanic plate of 80 km thickness underneath a 100 km thick continental plate is initiated with an 80 km slab.Different compositional fields are used to describe the oceanic and continental crusts and the upper mantle (UM), transition zone (TZ), and lower mantle (LM) material.A 2-times viscosity increase is also included at the 660 km phase boundary.The top boundary is a true free surface; the right vertical boundary has a prescribed in-and outflow as indicated from 8 Myr onward.A similar flow (4 cm yr −1 inflow) is prescribed on the left boundary until 13 Myr, when the boundary is opened.

Figure B2 .Figure B3 .
Figure B2.Compressible subduction evolution in terms of viscosity.Density is contoured at 3700 and 4200 kg m −3 , demonstrating phase boundary topography in the slab.

Table 1 .
Definition of symbols.

Table 2 .
Characteristics of performed experiments.
* nc: number of compositions.None of the benchmarks include temperature effects in the rheology.

Table 3 .
The indentor benchmark model parameters.
An estimate of the maximum strain rate and thus minimum viscosity can be made from µ min ≈ σy

Table 4 .
The brick benchmark model parameters.

Table 5 .
The sandbox experiment model parameters.

Table 6 .
The detachment benchmark model parameters.

Table A2 .
Wall time for time step 2000 of the self-consistent subduction benchmark for different viscosity averaging methods using 28 cores.