We model deformation in the subduction channel (SC) and overriding plate as two-dimensional incompressible, Newtonian, isotropic, Stokes viscoelastic creeping flow under the influence of gravity in a Lagrangian finite element grid. SC flow in the trench-parallel direction is assumed to be small compared to flow along a plane normal to the trench axis. In this case, a 2D treatment can capture the predominant forces acting on the channel (i.e., the shear forces applied by the subducting plate and gravity-derived forces). Even if SC and upper plate materials may present short-scale (i.e., grain size) anisotropies in viscosity, we further assume, like many others, that there is no such anisotropy when averaged over longer scales. Since mantle and crustal motions can be considered as highly viscous creeping over geologic timescales, all inertial forces are negligible during creep. We thus use the Stokes approximation to fluid flow. Finally, rocks can typically recover part of their deformation once stress disappears, so we use a viscoelastic rather than a purely viscous rheology. For further simplicity in these initial models, we also assume a linear shear stress-shear rate relationship. A Newtonian rheology for the SC is consistent with the (micro)structural record of exhumed metamorphic rocks which indicates that deformation is localized in low-strength shear zones dominated by dissolution precipitation creep or fluid-assisted granular flow (Gerya and Stöckhert, 2002). The model geometry is a several hundred-meter-thick SC at the interface between the upper plate and subducting slab. This geometry is inferred from imaging of SCs as ~100-1000-m-thick low-velocity layers (e.g., Sage et al., 2006). The lower model domain boundary is 100 km, while the eastern model domain boundary is 400 km. The Moho and upper crust-lower crust geometries are prescribed according to the Andean density model of Tassara et al. (2006). Triangle C (Shewchuk, 1996) is used to generate ~105-element unstructured two-dimensional Delaunay triangular meshes after each time step for two subdomains: the SC and the upper plate. A ~250 m horizontal resolution is achieved along the free surfaces. We follow deformation of the atmosphere-land and seafloor free surfaces using 10 kyr time steps, wherein deformation results from shear stresses imposed on the SC top by the subducting slab, by flow in the SC, by buoyancy forces originating from differences in density between SC and forearc materials, and by the weight of the water column in the case of the seafloor surface. We use MILAMIN, a MATLAB-based finite element method mechanical solver (Dabrowski et al., 2008), modified to include the viscoelasticity terms. The incompressibility constraint is achieved through an iterative penalty method. A relatively small penalty factor k is used, guaranteeing a good condition number, and incompressibility of flow is achieved by Powell/Hestenes/Uzawa-type iterations (Dabrowski et al., 2008), that is also known as an Augmented Lagrangian formulation.

In order to run the code, it is necessary to upload to the local environment the SuitSparse package functions, as well as the MUTILS 4.2 functions (both not included here).