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 (< or =1,500 km) retreat fast and develop a curved geometry, concave towards the mantle wedge side. Models with slabs intermediate in width ( approximately 2,000-3,000 km) are sublinear and retreat more slowly. Models with wide slabs (> or =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.

Plates are an integral part of the convection system in the fluid mantle, but plate boundaries are the product of brittle faulting and plate motions are strongly influenced by the existence of such faults. The conditions for plate tectonics are studied by considering brittle behaviour, using Byerlee’s law to limit the maximum stress in the lithosphere, in a mantle convection model with temperature‐dependent viscosity. When the yield stress is high, convection is confined below a thick, stagnant lithosphere. At low yield stress, brittle deformation mobilizes the lithosphere which becomes a part of the overall circulation; surface deformation occurs in localized regions close to upwellings and downwellings in the system. At intermediate levels of the yield stress, there is a cycling between these two states: thick lithosphere episodically mobilizes and collapses into the interior before reforming. The mobile‐lid regime resembles convection of a fluid with temperature‐dependent viscosity and the boundary‐layer scalings are found to be analogous. This regime has a well defined Nusselt number–Rayleigh number relationship which is in good agreement with scaling theory. The surface velocity is nearly independent of the yield stress, indicating that the ‘plate’ motion is resisted by viscous stresses in the mantle. Analysis suggests that mobilization of the Earth’s lithosphere can occur if the friction coefficient in the lithosphere is less than 0.03–0.13—lower than laboratory values but consistent with seismic field studies. On Venus, the friction coefficient may be high as a result of the dry conditions, and brittle mobilization of the lithosphere would then be episodic and catastrophic.

Abstract. Layered viscosity, temperature-dependent viscosity, and surface plates have an important effect on the scale and morphology of structure in spherical models of mantle convection. We find that long-wavelength structures can be produced either by a layered viscosity with a weak upper mantle or temperature-dependent viscosity even in the absence of surface plates, corroborating earlier studies. However, combining the layered viscosity structure with a temperature-dependent viscosity results in structure with significantly shorter wavelengths. Our models show that the scale of convection is mainly controlled by the surface plates, supporting the previous two-dimensional studies. Our models with surface plates, layered and temperaturedependent viscosity, and internal heating explain mantle structures inferred from seismic tomography. The models show that hot upwellings initiate at the core-mantle boundary (CMB) with linear structures, and as they depart from CMB, the linear upwellings quickly change into quasi-cylindrical plumes that dynamically interact with the ambient mantle and surface plates while ascending through the mantle. A linear upwelling structure is generated again at shallow depths (<200 km) in the vicinity of diverging plate margins because of the surface plates. At shallow depths, cold downwelling sheets form at converging plate margins. The evolution of downwelling sheets depends on the mantle rheology. The temperature-dependent viscosity strengthens the downwelling sheets so that the sheet structure can be maintained throughout the mantle. The tendency for linear upwelling and downwelling structures to break into plume-like structures is stronger at higher Rayleigh numbers. Our models also show that downwellings to first-order control surface plate motions and the locations and horizontal motion of upwellings. Upwellings tend to form at stagnation points predicted solely from the buoyancy forces of downwellings. Temperature-dependent viscosity greatly enhances the ascending velocity of developed upwelling plumes, and this may reduce the influence of global mantle flow on the motion of plumes. Our results can explain the anticorrelation between hotspot distribution and fast seismic wave speed anomalies in the lower mantle and may also have significant implications to the observed stationarity of hotspots.

Abstract.Small

[1] 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(Schellart, , 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).

Previous experimental studies of convection in fluids with temperature-dependent viscosity reached viscosity contrasts of the order of 105. Although this value seems large, it still might not be large enough for understanding convection in the interiors of Earth and other planets whose viscosity is a much stronger function of temperature. The reason is that, according to theory, above 104–105 viscosity contrasts, convection must undergo a major transition—to stagnant lid convection. This is an asymptotic regime in which a stagnant lid is formed on the top of the layer and convection is driven by the intrinsic, rheological, temperature scale, rather than by the entire temperature drop in the layer. A finite element multigrid scheme appropriate for large viscosity variations is employed and convection with up to 1014 viscosity contrasts has been systematically investigated in a 2D square cell with free-slip boundaries. We reached the asymptotic regime in the limit of large viscosity contrasts and obtained scaling relations which are found to be in good agreement with theoretical predictions.

We review the methods available for largedeformation simulations of geomaterials before presenting a La grangian integration point finite element method de signed specifically to tackle this problem. In our EL LIPSIS code, the problem domain is represented by an Eulerian mesh and an embedded set of Lagrangian in tegration points or particles. Unknown variables are computed at the mesh nodes and the Lagrangian parti cles carry history variables during the deformation pro cess. This method is ideally suited to model fluidlike behavior of continuum solids which are frequently en countered in geological contexts. We present bench mark examples taken from the geomechanics area.

[1] As high-performance computing facilities and sophisticated modeling software become available, modeling mantle convection in a three-dimensional (3-D) spherical shell geometry with realistic physical parameters and processes becomes increasingly feasible. However, there is still a lack of comprehensive benchmark studies for 3-D spherical mantle convection. Here we present benchmark and test calculations using a finite element code CitcomS for 3-D spherical convection. Two classes of model calculations are presented: the Stokes' flow and thermal and thermochemical convection. For Stokes' flow, response functions of characteristic flow velocity, topography, and geoid at the surface and core-mantle boundary (CMB) at different spherical harmonic degrees are computed using CitcomS and are compared with those from analytic solutions using a propagator matrix method. For thermal and thermochemical convection, 24 cases are computed with different model parameters including Rayleigh number (7 Â 10 3 or 10 5 ) and viscosity contrast due to temperature dependence (1 to 10 7 ). For each case, time-averaged quantities at the steady state are computed, including surface and CMB Nussult numbers, RMS velocity, averaged temperature, and maximum and minimum flow velocity, and temperature at the midmantle depth and their standard deviations. For thermochemical convection cases, in addition to outputs for thermal convection, we also quantified entrainment of an initially dense component of the convection and the relative errors in conserving its volume. For nine thermal convection cases that have small viscosity variations and where previously published results were available, we find that the CitcomS results are mostly consistent with these previously published with less than 1% relative differences in globally averaged quantities including Nussult numbers and RMS velocities. For other 15 cases with either strongly temperature-dependent viscosity or thermochemical convection, no previous calculations are available for comparison, but these 15 test calculations from CitcomS are useful for future code developments and comparisons. We also presented results for parallel efficiency for CitcomS, showing that the code achieves 57% efficiency with 3072 cores on Texas Advanced Computing Center's parallel supercomputer Ranger.