Dave A. May

Evolution and diversity of subduction zones controlled by slab width.

Subducting slabs provide the main driving force for plate motion and flow in the Earth’s mantle, and geodynamic, seismic and geochemical studies offer insight into slab dynamics and subduction-induced flow. Most previous geodynamic studies treat subduction zones as either infinite in trench-parallel extent (that is, two-dimensional) or finite in width but fixed in space. Subduction zones and their associated slabs are, however, limited in lateral extent (250–7,400 km) and their three-dimensional geometry evolves over time. Here we show that slab width controls two first-order features of plate tectonics—the curvature of subduction zones and their tendency to retreat backwards with time. Using three-dimensional numerical simulations of free subduction, we show that trench migration rate is inversely related to slab width and depends on proximity to a lateral slab edge. These results are consistent with retreat velocities observed globally, with maximum velocities (6–16 cm yr-1) only observed close to slab edges (<1,200 km), whereas far from edges (>2,000 km) retreat velocities are always slow (<2.0 cm yr-1). Models with narrow slabs (≤1,500 km) retreat fast and develop a curved geometry, concave towards the mantle wedge side. Models with slabs intermediate in width (∼2,000–3,000 km) are sublinear and retreat more slowly. Models with wide slabs (≥4,000 km) are nearly stationary in the centre and develop a convex geometry, whereas trench retreat increases towards concave-shaped edges. Additionally, we identify periods (5–10 Myr) of slow trench advance at the centre of wide slabs. Such wide-slab behaviour may explain mountain building in the central Andes, as being a consequence of its tectonic setting, far from slab edges.

Influence of trench width on subduction hinge retreat rates in 3-D models of slab rollback

Subduction of tectonic plates limited in lateral extent and with a free-trailing tail, i.e., “free subduction,” is modeled in a three-dimensional (3-D) geometry. The models use a nonlinear viscoplastic rheology for the subducting plate and exhibit a wide range of behaviors depending on such plate characteristics as strength, width, and thickness. We investigate the time evolution of this progressive rollback subduction, measure the accompanying return flow in the upper mantle, and quantify the plate kinematics. Due to the 3-D geometry, flow is allowed to accompany slab rollback around the lateral edges of the slab (the toroidal component), as opposed to 2-D geometry, where material is forced to flow underneath the slab tip (the poloidal component). A simple force balance is provided which relates the speed of backward trench migration to the resistive forces of generating flow and weakening the plate. Our results indicate most of the gravitational energy of the system (i.e., the negative buoyancy of the subducting slab) is converted into a toroidal flow (∼69%), a much smaller amount goes into weakening the plate (∼18%), and the remaining amount goes into driving flow parallel to displacement of the slab (∼13%). For the trench widths (W) we investigate (≤1500 km), a maximum trench retreat rate occurs for trenches 600 km wide, which is attributed to the interaction between a plate of finite width and the induced flow (which has a lengthscale in the horizontal direction). These numerical results quantitatively agree with comparable 3-D laboratory experiments using analogue models with a purely viscous plate material (Schellart, 2004a, 2004b), including correlations between increasing retreat rate with increasing plate thickness, trench width for maximum retreat rate (500 km), and estimated amount of slab buoyancy used to drive rollback-induced flow (∼70%). Several implications for plate tectonics on Earth result from these models such as rollback subduction providing a physical mechanism for ephemeral slab graveyards situated above the more viscous lower mantle (and endothermic phase transition) prior to a flushing event into the lower mantle (mantle avalanche).

Discretization errors and free surface stabilization in the finite difference and marker‐in‐cell method for applied geodynamics: A numerical study

The finite difference–marker-in-cell (FD-MIC) method is a popular method in thermomechanical modeling in geodynamics. Although no systematic study has investigated the numerical properties of the method, numerous applications have shown its robustness and flexibility for the study of large viscous deformations. The model setups used in geodynamics often involve large smooth variations of viscosity (e.g., temperature-dependent viscosity) as well large discontinuous variations in material properties (e.g., material interfaces). Establishing the numerical properties of the FD-MIC and showing that the scheme is convergent adds relevance to the applications studies that employ this method. In this study, we numerically investigate the discretization errors and order of accuracy of the velocity and pressure solution obtained from the FD-MIC scheme using two-dimensional analytic solutions. We show that, depending on which type of boundary condition is used, the FD-MIC scheme is a second-order accurate in space as long as the viscosity field is constant or smooth (i.e., continuous). With the introduction of a discontinuous viscosity field characterized by a viscosity jump (η*) within the control volume, the scheme becomes first-order accurate. We observed that the transition from second-order to first-order accuracy will occur with only a small increase in the viscosity contrast (η* ≈ 5). We have employed two methods for projecting the material properties from the Lagrangian markers onto the Eulerian nodes. The methods are based on the size of the interpolation volume (4-cell, 1-cell). The use of a more local interpolation scheme (1-cell) decreases the absolute velocity and pressure discretization errors. We also introduce a stabilization algorithm that damps the potential oscillations that may arise from quasi free surface calculations in numerical codes that employ the strong form of the Stokes equations. This correction term is of particular interest for topographic modeling, since the surface of the Earth is generally represented by a free surface. Including the stabilization enables physically meaningful solutions to be obtained from our simulations, even in cases where the time step value exceeds the isostatic relaxation time. We show that including the stabilization algorithm in our FD stencil does not affect the convergence properties of our scheme. In order to verify our approach, we performed time-dependent simulations of free surface Rayleigh-Taylor instability.

Inversion of fluvial channels for paleorock uplift rates in Taiwan

The transient response of erosion to changes in rock uplift rate leads to the preservation of rock uplift history in the long profiles of rivers. However, extracting this information is nontrivial as changes in channel steepness are the result of both spatial and temporal changes in rock uplift rate, as well as other factors such as climate and rock type. We exploit an analytical linear solution for river channel profile evolution in response to erosion and tectonic uplift to investigate the rock uplift history of Taiwan. The analytical approach allows us to solve the linear inverse problem, efficiently extracting rock uplift as a function of space and time, from digital elevation data. We assess the potential of fluvial topography to resolve rock uplift rates using three approaches: (1) a synthetic resolution test, (2) analysis of the forward model to demonstrate where in space and time the fluvial topography constrains rock uplift rate, and (3) interpretation of the model resolution matrix. Furthermore, the potential to analyze large data sets reduces the influence of stochastic processes such as landslides, small-scale river network reorganization, and also local lithological variability. In Taiwan, our analysis suggests that current rock uplift rates exceed erosion rates across much of the island and that there has been an increase in rock uplift rates since 0.5 Ma across the Central Range.

Development of a Stokes flow solver robust to large viscosity jumps using a Schur complement approach with mixed precision arithmetic

We develop an iterative solution technique for solving Stokes flow problems with smooth and discontinuous viscosity structures using a three dimensional, staggered grid finite difference discretization. Two preconditioned iterative methodologies are applied to the saddle point arising from the discrete Stokes problem. They consist of a velocity-pressure coupled approach (FC) and a decoupled, Schur complement approach (SC). Within both of these methods, we utilize either the scaled BFBt, or an identity matrix scaled by the local cell viscosity (LV) to define a preconditioner for the Schur complement. Additionally, we propose to use a mixed precision Krylov kernel to improve the convergence by reducing round-off error. In this approach, standard double precision is used during the application of the preconditioner, whilst higher precision arithmetic is used to define the matrix vector product, dot products and norms required by the Krylov method. In our Krylov kernel, we utilize quad precision arithmetic which is emulated via the double-double precision method. We consider several simplified geodynamic problems with a viscosity contrast to demonstrate the robustness and scalability of our solution methods. Through a careful choice of stopping conditions, we are able to quantitatively compare the residuals between the SC and FC approaches. We examine the trade-off relationship between the number of outer iterations required for convergence, and the computational cost per iteration, for the each solution methods. We find that it is advantageous to use the FC approach utilizing relaxed tolerances for solution of the sub-problems, combined with the LV preconditioner. We also observed that in general, the SC approach is more robust than FC and that BFBt is more robust than LV when used in our numerical experimental. In addition, our mixed precision method produces improved convergence rates of Arnoldi type Krylov subspace methods without a drastic increasing the computational time. The usage of a high precision Krylov kernel is found to be useful for the solver associated with the velocity sub-problem.

pTatin3D: high-performance methods for long-term lithospheric dynamics

Simulations of long-term lithospheric deformation involve post-failure analysis of high-contrast brittle materials driven by buoyancy and processes at the free surface. Geodynamic phenomena such as subduction and continental rifting take place over millions year time scales, thus require efficient solution methods. We present pTatin3D, a geodynamics modeling package utilising the material-point-method for tracking material composition, combined with a multigrid finite-element method to solve heterogeneous, incompressible visco-plastic Stokes problems. Here we analyze the performance and algorithmic tradeoffs of pTatin3Ds multigrid preconditioner. Our matrix-free geometric multigrid preconditioner trades flops for memory bandwidth to produce a time-to-solution > 2× faster than the best available methods utilising stored matrices (plagued by memory bandwidth limitations), exploits local element structure to achieve weak scaling at 30% of FPU peak on Cray XC-30, has improved dynamic range due to smaller memory footprint, and has more consistent timing and better intra-node scalability due to reduced memory-bus and cache pressure.

Discretization Errors in the Hybrid Finite Element Particle-in-cell Method

In computational geodynamics, the Finite Element (FE) method is frequently used. The method is attractive as it easily allows employment of body-fitted deformable meshes and a true free surface boundary condition. However, when a Lagrangian mesh is used, remeshing becomes necessary at large strains to avoid numerical inaccuracies (or even wrong results) due to severely distorted elements. For this reason, the FE method is oftentimes combined with the particle-in-cell (PIC) method, where particles are introduced which track history variables and store constitutive information. This implies that the respective material properties have to be interpolated from the particles to the integration points of the finite elements. In numerical geodynamics, material parameters (in particular the viscosity) usually vary over a large range. This may be due to strongly temperature-dependent rheologies (which result in large but smooth viscosity variations) or material interfaces (which result in viscosity jumps). Here, we analyze the accuracy and convergence properties of velocity and pressure of the hybrid FE-PIC method in the presence of large viscosity variations. Standard interpolation schemes (arithmetic and harmonic) are compared to a more sophisticated interpolation scheme which is based on linear least squares interpolation for two types of elements (

Influences of surface processes on fold growth during 3‐D detachment folding

Q_1P_0