Implementing nonlinear viscoplasticity in ASPECT : benchmarking and applications to 3 D subduction modeling

ASPECT (Advanced Solver for Problems in Earth’s ConvecTion) is a massively parallel finite 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 nonlinear 5 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 fields to represent different materials (and, since recently, tracers), all material parameters are made dependent on a user-specified number of fields. The goal of this paper is primarily to describe and validate 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 10 sandbox experiment and the slab detachment benchmark. Furthermore, we aim to provide hands-on examples for prospective users by demonstrating the use of multi-material viscoplasticity with three-dimensional, thermomechanical models of oceanic subduction, putting ASPECT on the map as a community code for high-resolution, nonlinear rheology subduction modeling.


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 non-elastic 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).numerical sandbox (Buiter et al., 2006) and the slab detachment benchmark (Schmalholz, 2011;Thieulot et al., in prep.).
Finally, in Section 4 we present two 3D subduction applications to showcase the new suite of possibilities made available through our additions and adaptations, and we discuss our overall results in Section 5.

Methods
A short summary of the governing equations solved by ASPECT is given in Section 2.1 (for more information the reader is referred to Kronbichler et al., 2012).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 (1), mass (2) and energy (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 (4) to the system of equations ( 1)-(3) described above.As these equations contain no natural diffusion, artificial diffusivity ν h is again introduced to stabilise advection: Note that as of 2016 it is also possible to use active as well as passive tracers in ASPECT (version 1.4.0).

Solving the governing equations
ASPECT solves the above equations using the finite element method: the domain is discretized into quadrilateral/hexahedral 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, temperature and composition and first order polynomials for pressure (Q 2 Q 1 elements, e.g.Donea and Huerta, 2003).Unless stated otherwise, these default polynomial degrees are used in the following.The linearized Stokes system is solved in a procedure involving the iterative FGMRES solver with an inexact right preconditioner.For details on the construction of the preconditioner, see Kronbichler et al. (2012).The CG method with an incomplete LU decomposition preconditioner is used for the temperature ).The initial residual ||A(x 0 )x 0 − b|| 2 is computed with zero velocities.

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 time dependent 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 ˙ (u), T , P and c i as well as 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 at 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

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 Table 1 for the definition of used symbols.The effective deviatoric strain rate is defined as ˙ e = 1 2 ˙ ij ˙ ij .In this paper we report model set-up values for simplified prefactor B, defined as 1 Ab m 1 n , and add a scaling factor β to Eq. ( 5) to easily tune the effective viscosity: 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 (VRM) (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 internal friction angle φ is zero, this criterion reverts back to the von Mises criterion.
Both types of viscous creep act simultaneously (Karato, 2008) under the same deviatoric stress, so their contributions 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 timestep, 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 timesteps, the strain rate of the previous timestep is used as an initial guess for the iterative process.
The final effective viscosity µ vp ef f is capped by the user-defined minimum viscosity µ min and maximum viscosity µ max to avoid extreme excursions and to ensure stability of the numerical scheme: We have successfully run the models presented here with µ max /µ min 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 (viscosity, density and other).
We have implemented the four averaging schemes commonly referred to in the literature (e.g.Deubelbeiss and Kaus, 2008;Schmeling et al., 2008): where nc is the total number of compositional fields c i in the domain.
The above methods 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 three 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 infinite norm rule in this paper; for a discussion of this choice, see Appendix A.

Nonlinear rheology benchmarks
To test and validate our implementation of multi-material viscoplastic rheologies, we performed four 2D 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 with ∼ 1, 000 cores (2.34GHz, 2.88Tb RAM memory, Qlogic InfiniBand).Wall times quoted can have changed with versions of ASPECT newer than those used for the described experiments.

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).The analytical solution (Fig. 1) is characterized by 3 observations: 1.The angles of the shear bands stemming from the edges of the indented area are 45 degrees.
2. The pressure at the surface in the centre of the punch (I) and the pressure in triangles ABC & EFG is P I = σ y (1+π) and P ABC = P EF G = σ y , respectively.

The velocity magnitude in areas CDE and ABDC
2 , respectively.

Model set-up
The numerical set-up of the instantaneous indentor benchmark comprises a 2D 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
Figure 3 shows the model results for a rough (left column) and smooth (right column) punch.The obtained solutions 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&e) and stem from the edges of the indentor at a 45 degrees angle with respect to the top of the medium (Fig. 3b&f).For a rough punch, pressure and velocity measurements deviate from the analytical solution by about 15% and 2%, respectively, but the block-like behavior of triangle CDE is evident (Fig. 3c).The analytical solution is exactly reproduced for a smooth punch (errors < 0.2%), but the velocity vectors in Fig. 3g show some horizontal motion of triangle CDE and the velocity field is more diffuse.This trade-off is as expected, because the horizontal component of surface velocity is left free for the smooth punch, while it is fixed to zero for the rough punch.

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 performed indentor experiment also shows a trade-off between accuracy in pressure and velocity measurements and the rigid-plastic like behaviour 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
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 set-up
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).Compression and extension are prescribed through kinematic boundary conditions on the vertical domain walls.The bottom boundary of the 40 × 10km 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.4km, x = 19.4km,x = 20.6km and x = 22.6km (Kaus, 2010).Lower resolution runs were performed, but do not resolve the viscous inclusion well ( 2 elements) and are not shown (see instead Fig. 12 of Kaus, 2010).We monitor the relative residual as a measure of convergence, but the number of Nonlinear Iterations (NI, see Section 2.1.2) is fixed at 1,000.The red symbols in Fig. 5 indicate runs for which the residual is not monotonously decreasing (after the first peak in residual), as is evident from Fig. 6 (black versus red lines).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 1, 024 × 256 elements.

Discussion
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 plugin includes frictional plasticity.Testing this pressure-dependent plasticity 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 consistent Coulomb angles can be achieved by an (initially) associated flow law (where φ = ψ).
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 Solid Earth Discuss., doi:10.5194/se-2017-9,2017 Manuscript under review for journal Solid Earth Discussion started: 9 February 2017 c Author(s) 2017.CC-BY 3.0 License.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 6 different numerical codes and compared to the analog results described in Schreurs et al. (2006).This experiment has previously been repeated by Thieulot (2011) and -in a symmetrical versionby Gerya et al. (2013).

Model set-up
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 centre of the domain.We model this set-up with 3 compositional fields (Fig. 8): 1) a viscous basal layer, 2) an overlying Drucker-Prager dynamic pressure-dependent 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 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 open.Adaptive Mesh Refinement (AMR) is applied based on the effective strain rate field to obtain a maximum local refinement of 0.39 × 0.39mm along the shear bands (Fig. 10).The model is run until 2cm of extension has occurred, equalling 2, 880s of model time.Material properties and other model parameters are listed in Table 5.

Model results
The results for 1 and 2cm 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 2cm of extension, the angle of the shear bands 5 left and right of the velocity discontinuity (see Fig. 9d) measure ∼ 51 • ,∼ 62 • and ∼ 54 • ,∼ 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 10 the strain rate 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 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).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.These are numerical effects tied to finite element models that should be taken into consideration when interpreting and comparing model results.

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

Model set-up
The 2D 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 Pas 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 Pas.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 80km 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 20My, 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 ignores 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. 12e&f).
It can be seen from the graph in Fig. 13 that the width of the necked zone through time agrees very well with the results of Schmalholz (2011).Only after t = 0.6, our results deviate slightly, showing faster necking.This could be the result of the different averaging technique used here, as the infinite norm ignores the higher-viscosity contribution of the slab with respect to the other averaging methods (see Eq. ( 15)-( 18) and Appendix A).ASPECT 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.

Discussion
The three previous benchmarks focussed 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 5 with that of Schmalholz (2011) and other codes (Hillebrand et al., 2014;Thieulot et al., in prep.).It also demonstrates the 17 Solid Earth Discuss., doi:10.5194/se-2017-9,2017 Manuscript under review for journal Solid Earth Discussion started: 9 February 2017 c Author(s) 2017.CC-BY 3.0 License.effective splitting of a compositional field into two bodies and how the rheology interacts with the compositional fields through localisation of deformation at the slab hinge.It should be noted that for this type of rheology the solution is mesh-resolution dependent as well, and iterative convergence should be monitored as for plastic rheologies.Differences in model evolution can also arise from the particular viscosity averaging method applied, as will be discussed in Thieulot et al. (in prep.), who present a 3D detachment benchmark based on Schmalholz (2011).

Model set-up
We discuss two models; the first is an adaptation of the free-plate model of Schellart and Moresi (2013) and 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 temperature, pressure and strain rate dependent viscosity, 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 tip of the SP extends into the mantle for 200km.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 fixed, 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 1, 593K and surface temperature   7).(a) Model 1: a free Overriding Plate (OP) and a Subducting Plate (SP) with a trench at x = 2, 400km.At start-up, the slab extends 200km into the mantle (measured vertically) at an angle of 29 • .The 100km thick SP consists of 2 compositional layers, the OP of only 1 composition.All boundaries are free slip, except the no slip bottom boundary.(b) Model 2: A 58.8My old Adjacent Plate (AP) is added, separated from the OP and SP by a Weak Zone (WZ) of 20km 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 16My for the OP and 60My for the SP, resulting in thicknesses of 50km and 100km at the trench, respectively.An adiabat of 0.25Kkm −1 is prescribed in the mantle.The bottom temperature is fixed at 1, 728K, the top at 293K.

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 10 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 The AMR based on composition and viscosity follows the outline of the plates as they move through the mantle (Fig. 17

Model 2
The SP of model 2 steepens for the first 12My (Fig. 18) until full-fledged subduction starts.Subduction is much slower than in model 1: while in model 1 subduction is completed by 35My, rollback of the slab in model 2 only sets in around 50My.
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 49My).Even though the slab reaches the bottom boundary around 50My and starts to shallow, it does not lie flat on the 660km boundary, but continues to hover above it.Also note the halo of increased strain rate around the tip of the slab, and the small-scale strain rate features in the mantle compared to model 1.

Discussion
The evolution of model 1 strongly resembles that of Schellart and Moresi (2013)  of the slab at the 660km 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 660km boundary.In switching from model 1 to this thermomechanically coupled model 2, we found changes to the set-up 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 increased the coupling of the OP plate to the left boundary unfortunately, 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).Therefore, the differences derive from the temperature, pressure, strain rate and composition dependent rheology.Indeed, 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 (cut-off at 10 24 Pas).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 higherviscosity 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 2D model in Appendix B.

Discussion
The four benchmarks shown using our viscoplasticity implementations in ASPECT either reproduce the available analytical solution (Section 3.1), or compare well with theory (Section 3.2) or the results of other codes (Sections 3.2-3.4,see also Tosi 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, 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.Another addition would be to include plastic softening or hardening -changes in the yield surface due to the accumulated strain.As of December 2016, the ASPECT developer version provides a way to use either tracer particles or compositional fields to track the strain and weaken both the cohesion and friction angle (ASPECT, 2016).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 efficient Newton iterations (see for example Popov and Sobolev, 2008;May et al., 2015;Rudi et al., 2015) could help prevent nonlinear rheologies forming a bottleneck.Such iterations have been implemented in ASPECT (Fraters et al., in prep.) and are currently being incorporated in the release version.

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 4D 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 arises the need 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, we have outlined here our basic additions that enable viscoplastic crust and lithosphere modeling.We then tested and validated the algorithms with four different benchmarks well-known in the geodynamic modeling community.Last, we highlighted the possibilities arising from the adaptations with 4D thermomechanically coupled viscoplas-Solid Earth Discuss., doi:10.5194/se-2017-9,2017 Manuscript under review for journal Solid Earth Discussion started: 9 February 2017 c Author(s) 2017.CC-BY 3.0 License.tic 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 evergrowing 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) and active particles (Gassmöller et al., 2016, submitted).Also, the extensive user manual (Bangerth et al., 2016) 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.

Code availability
The plasticity formulation has become part of the ASPECT distribution.A pull request for a plugin that calculates a constant lithostatic pressure profile needed for open boundary traction conditions is available from GitHub and will be incorporated in ASPECT.Input parameter files to reproduce the benchmarks will be incorporated as well.

Data availability
There is no data to be made available.
Appendix A: Self-consistent subduction and compositional averaging As discussed in Section 2.2.2, the choice of averaging method for models of multiple compositions can significantly influence the results (Deubelbeiss and Kaus, 2008;Schmeling et al., 2008).Based on a number of experiments, we have chosen the infinite norm as method of choice for this paper.One of the models on which this choice is based, is that of Schmeling et al. (2008), on which we will elaborate below.

A1 Model set-up
The 2D linear viscous model is composed of three compositions: the mantle, subducting lithosphere and sticky-air to allow for surface topography build-up and detachment of the lithosphere from the top boundary (Fig. 20, Table 9).The subducted part of the lithosphere supplies the force to start subduction.The four different averaging methods in Eq. ( 15)-( 18)) are tested at different resolutions.

A2 Model results
The evolution of subduction is summarized in a plot of the slab tip depth over time in Fig. 21.It is clear that the effect of averaging method dominates over that of resolution, but that both are significant.As do Schmeling et al. (2008)   compressible formulation of the governing equations is used, including shear heating, adiabatic heating and latent heat: where β is the compressibility 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 10 (e.g.Christensen and Yuen, 1985):

B1 Model set-up
The compressible subduction model considers ocean-continent subduction in a domain of 1, 600km by 1, 000km (Fig. 23).The model includes phase changes around 410km and 660km depth, see  At 200km, 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

Wall time 16h
Lithosphere

Figure 3 .
Figure 3.The punch benchmark results after 500 NI for a rough punch (left column) and a smooth punch (right column).(a) & (e): Viscosity field with analytical slip lines.(b) & (f): Strain rate field with measured shear band angles.(c) & (g): Velocity field with velocity vectors along the surface of the domain and velocity measurements in points K and L. (d) & (h): Pressure along the surface of the domain (colored line) and analytical solution values π + 1 and 1 (grey lines).Rough punch: PI = 4.7348 and PH = PJ = 0.6983.Smooth punch: PI = 4.1413 and PH = PJ = 0.9997.

Figure 4 .
Figure 4.The brick benchmark model set-up 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 × 400m is placed at the bottom of the domain.

Figure 5
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 1, 024×256 elements.

Figure 5 .
Figure 5. 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 1, 024 × 256 elements (dark grey).All models were run for 1, 000 NI; red symbols indicate runs with convergence behavior that is not monotonous.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 el.extension case with φ = 25 • -because it was computed using two different shear bands.

Figure 6 .
Figure 6.Measured relative residual 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 fields for (a) φ = 0 • and (b) φ = 30 • for a 512 × 128 elemental resolution.The models were run in extension for 1, 000 NI. Black lines indicate the theoretical Coulomb angle θC .Measured shear band angles θM are also given.

Figure 8 .
Figure 8.The sandbox experiment set-up after Buiter et al. (2006).VD: 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.5cm.

Figure 10 .
Figure 10.Numerical grid after ∼ 2cm of extension of the numerical sandbox.Adaptive mesh refinement and coarsening based on the strain rate and density leads to an elemental resolution varying from 512 × 128 to 32 × 8 elements.
3.3.3DiscussionThe evolution of the numerical sandbox model -a model with AMR, high viscosity contrasts, large deformation and complex boundary conditions -compares well with those shown inBuiter et al. (2006) andThieulot (2011).Although the right shear band angles of 62 • and 60 • after 1 and 2cm 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 Section 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 converge to Coulomb angles with increasing mesh resolution.Differences in the measured angle of shear bands are therefore not surprising.A lower resolution run (maximum of 256 × 64 elemental resolution, not shown) results in angles to the right of the discontinuity of 60 • and 55 • .
Solid Earth Discuss., doi:10.5194/se-2017-9,2017 Manuscript under review for journal Solid Earth Discussion started: 9 February 2017 c Author(s) 2017.CC-BY 3.0 License.(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 3D by von Tscharner et al. (2014).

Figure 11 .
Figure 11.The detachment benchmark model set-up 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 boundary are free slip, while the vertical boundaries are no slip.20 Passive tracers are placed along the outline of the slab.

Figure 12 .
Figure 12.Detachment benchmark model evolution showing the viscosity field over time.After about 20My, necking is complete and the remaining lithosphere reaches a high, uniform viscosity.

Figure 13 .
Figure 13.Nondimensional necking width versus time for ASPECT and Schmalholz (2011).The ASPECT necking width is calculated Solid Earth Discuss., doi:10.5194/se-2017-9,2017 Manuscript under review for journal Solid Earth Discussion started: 9 February 2017 c Author(s) 2017.CC-BY 3.0 License.T 0 = 293K (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 Ta−T Ta−T0 = 0.1.The maximum thickness of the plate (for time to infinity) is set to 125km, 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 100km for an age of ∼ 60My.From a depth of 125km, a linear temperature increase is prescribed everywhere to a bottom temperature of 1, 771K.The WZ and AP are of constant viscosity, but all other compositions are of nonlinear rheology, with which model 5 2 differs strongly from model 1.The specific rheological parameters for diffusion and dislocation creep and Drucker-Prager plasticity can be found in Table8.

Figure 14 .
Figure 14.3D Subduction model set-ups (also, see Table7).(a) Model 1: a free Overriding Plate (OP) and a Subducting Plate (SP) with a

Figure 15 .
Figure 15.Temperature distribution with depth and distance from the ridge for (a) the OP and (b) the SP.
), resulting in a local resolution of roughly 3km while the mantle is resolved with ∼ 50km elements.20 Solid Earth Discuss., doi:10.5194/se-2017-9,2017 Manuscript under review for journal Solid Earth Discussion started: 9 February 2017 c Author(s) 2017.CC-BY 3.0 License.
, on whose set-up the model is based: the slab first sinks freely until it reaches the 660km bottom boundary and then starts draping while rolling back.The absence of folding 21 Solid Earth Discuss., doi:10.5194/se-2017-9,2017 Manuscript under review for journal Solid Earth Discussion started: 9 February 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 17 .
Figure 17.Model 1 -strain rate field 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 field 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.
et al. (2015) for a benchmark with ASPECT using a different viscoplastic formulation).Thus validating our implementations, Solid Earth Discuss., doi:10.5194/se-2017-9,2017 Manuscript under review for journal Solid Earth Discussion started: 9 February 2017 c Author(s) 2017.CC-BY 3.0 License.they allowed us to set up a 4D model of oceanic subduction, exemplifying the functionality that our implementations have added.

Figure 21 .
Figure 21.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 2, 048×512 elements.Minimum resolution is always 128×32 elements.The dashed red line model has a resolution varying from 128×128 to 2, 048 × 2, 048 elements.Shaded areas represent results of Schmeling et al. (2008, Fig. 6).

Figure 22 .
Figure 22.Viscosity field for each averaging method at 2, 048×512 elements local resolution at similar moments in the subduction evolution.

Figure 23 .
Figure 23.Compressible subduction model set-up: Subduction of an oceanic plate of 80km thickness underneath a 100km thick continental plate is initiated with an 80km slab.Different compositional fields are used to describe the oceanic and continental crust, and the upper mantle, transition zone and lower mantle material.A 2-times viscosity increase is also included at the 660km phase boundary.The top boundary is a true free surface, the right vertical boundary has a prescribed in/outflow as indicated.A similar flow (4cmyr −1 inflow) is prescribed on the left boundary until 13My, when the boundary is opened.
Figure24shows the evolution of the compressible subduction model in terms of viscosity.When the left boundary is opened at 13My, just after the slab has reached the 410km phase boundary, sinking velocities increase to 10.5cmyr −1 .Upon reaching the 660km phase boundary (with a two-fold viscosity increase), the slab slows down to about 5cmyr −1 .The tip of the slab is impeded by the 660km boundary and moves at around 1cmyr −1 along the boundary.The change in flow through the left boundary is depicted in Fig.25: 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 11cmyr −1 , together with an increase in transition and lower mantle velocity.When the slab reaches the 660km phase boundary around 16My, lower mantle outflow goes up rather uniformly throughout the lower mantle, but in/outflow above decreases.When the slab tip moves over the 660km 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.24show positive topography in the slab at the 410km discontinuity, while the 660km transition in the slab occurs deeper.The deformation of the surface at 18My shows a lowered subducting plate (about 1km below the initial top boundary), a topographic rise of continental crust above the subducting slab (of maximally 8km) and an elevated overriding continental plate (about 4km).

Table 11 .
An inward plate velocity is prescribed on the upper part of the left vertical boundary for the first 13My at 4cmyr −1 , while prescribed outward mantle velocities compensate for this volume increase.After 13My, 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 8My with a plate inflow of 2cmyr −1 .Apart from the subducting plate crust, which is of linear viscosity, all materials are nonlinear viscoplastic.

Table 2 .
Characteristics of performed experiments Hillebrand et al. (2014) †nc: number of compositions.None of the benchmarks include temperature effects in the rheology

Table 5 .
The sandbox experiment model parameters Parameters are on sandbox scale, as they are used in the model.36 Solid Earth Discuss., doi:10.5194/se-2017-9,2017 Manuscript under review for journal Solid Earth Discussion started: 9 February 2017 c Author(s) 2017.CC-BY 3.0 License.

Table 6 .
The detachment benchmark model parameters