Free Subduction Research Articles

The influence of viscous slab rheology on numerical models of subduction

Abstract. Numerical models of subduction commonly use diffusion and dislocation creep laws from laboratory deformation experiments to determine the rheology of the lithosphere. The specific implementation of these laws varies from study to study, and the impacts of this variation on model behavior have not been thoroughly explored. We run simplified 2D numerical models of free subduction in SULEC, with viscoplastic slabs following (1) a diffusion creep law, (2) a dislocation creep law, and (3) both simultaneously, as well as several variations of model 3 with reduced resistance to bending. We compare the results of these models to a model with a constant-viscosity slab to determine the impact of the implementation of different lithospheric flow laws on subduction dynamics. In creep-governed models, higher subduction velocity causes a longer effective slab length, increasing slab pull and asthenospheric drag, which, in turn, affect subduction velocity. Numerical and analogue models implementing constant-viscosity slabs lack this feedback but still capture morphological patterns observed in more complex models. Dislocation creep is the primary deformation mechanism throughout the subducting lithosphere in our models. However, both diffusion creep and dislocation creep predict very high viscosities in the cold core of the slab. At the trench, the effective viscosity is lowered by plastic failure, rendering effective slab thickness the primary control on bending resistance and subduction velocity. However, at depth, plastic failure is not active, and the viscosity cap is reached in significant portions of the slab. The resulting high slab stiffness causes the subducting plate to curl under itself at the mantle transition zone, affecting patterns in subduction velocity, slab dip, and trench migration over time. Peierls creep and localized grain size reduction likely limit the stress and viscosity in the cores of real slabs. Numerical models implementing only power-law creep and neglecting Peierls creep are likely to overestimate the stiffness of subducting lithosphere, which may impact model results in a variety of respects.

The Role of Slab Remnants in Modulating Free Subduction Dynamics: A 3‐D Spherical Numerical Study

AbstractSeismic tomography of Earth's mantle images abundant slab remnants, often located in close proximity to active subduction systems. The impact of such remnants on the dynamics of subduction remains underexplored. Here, we use simulations of multi‐material free subduction in a 3‐D spherical shell geometry to examine the interaction between visco‐plastic slabs and remnants that are positioned above, within and below the mantle transition zone. Depending on their size, negatively buoyant remnants can set up mantle flow of similar strength and length scales as that due to active subduction. As such, we find that remnants located within a few hundred km from a slab tip can locally enhance sinking by up to a factor 2. Remnant location influences trench motion: the trench advances toward a remnant positioned in the mantle wedge region, whereas remnants in the sub‐slab region enhance trench retreat. These motions aid in rotating the subducting slab and remnant toward each other, reducing the distance between them, and further enhancing the positive interaction of their mantle flow fields. In this process, the trench develops along‐strike variations in shape that are dependent on the remnant's location. Slab‐remnant interactions may explain the poor correlation between subducting plate velocities and subducting plate age found in recent plate tectonic reconstructions. Our results imply that slab‐remnant interactions affect the evolution of subducting slabs and trench geometry. Remnant‐induced downwelling may also anchor and sustain subduction systems, facilitate subduction initiation, and contribute to plate reorganization events.

Comparing the Dynamics of Free Subduction in Cartesian and Spherical Domains

AbstractThe effects of sphericity are regularly neglected in numerical and laboratory studies that examine the factors controlling subduction dynamics. Most existing studies have been executed in a Cartesian domain, with the small number of simulations undertaken in a spherical shell incorporating plates with an oversimplified rheology, limiting their applicability. Here, we simulate free‐subduction of composite visco‐plastic plates in 3‐D Cartesian and spherical shell domains, to examine the role of sphericity in dictating the dynamics of subduction, and highlight the limitations of Cartesian models. We identify two irreconcilable differences between Cartesian and spherical models, which limit the suitability of Cartesian‐based studies: (a) the presence of sidewall boundaries in Cartesian models, which modify the flow regime; and (b) the reduction of space with depth in spherical shells, alongside the radial gravity direction, which cannot be captured in Cartesian domains. Although Cartesian models generally predict comparable subduction regimes and slab morphologies to their spherical counterparts, there are significant quantitative discrepancies. We find that simulations in Cartesian domains that exceed Earth's dimensions overestimate trench retreat. Conversely, due to boundary effects, simulations in smaller Cartesian domains overestimate the variation of trench curvature driven by plate width. Importantly, spherical models consistently predict higher sinking velocities and a reduction in slab width with depth, particularly for wider subduction systems, enhancing along‐strike slab buckling and trench curvature. Results imply that sphericity must be considered for understanding the dynamics of Earth's wider subduction systems, and is already a significant factor for slabs of width 2,400 km.

Tectonic deformation at the outer rise of subduction zones

SUMMARY Fluids associated with subducting slabs play a crucial role in regulating the dynamics of water discharge, subsequent arc magmatism and intermediate-depth earthquakes in subduction zones. The incoming slab mantle hydration is primarily determined by deep normal faulting due to plate bending at the trench. However, the controlling factors on the outer rise faulting pattern, and the correlation between the inherited outer rise deformation and the intermediate-depth earthquakes, remain to be understood. Here we present high-resolution viscoelasto-plastic numerical models of free subduction for slab bending-related faulting prior to subduction. Our model results show that plastic weakening and friction coefficient of the slab mantle exhibit a significant impact on fault pattern, while plate age and elasticity have a minimal bearing for mature slabs. The brittle bending faults result in large positive pressure gradients in the vertical direction, facilitating seawater infiltrating into the subducting slabs, which corroborates previous numerical models. The faults reaching 15–30 km beneath the Moho coincide with the width of the double seismic zone in subduction zones. We anticipate that water pumped into the slab mantle along the faults, with decreasing water content along the depth, can explain the relatively sporadic lower plane earthquakes.

Geodynamic Modeling With Uncertain Initial Geometries

AbstractGeodynamic codes have become fast and efficient enough to facilitate sensitivity analysis of rheological parameters. With sufficient data, they can even be inverted for. Yet, the geodynamic inverse problem is often regularized by assuming a constant geometry of the geological setting (e.g., shape, location and size of salt diapirs or magma bodies) or approximating irregular bodies with simple shapes like boxes, spheres or ellipsoids to reduce the parameter space. Here, we present a simple and intuitive method to parameterize complex 3D bodies and incorporate them into geodynamic inverse problems. The approach can automatically create an entire ensemble of initial geometries, enabling us to account for uncertainties in imaging data. Furthermore, it allows us to investigate the sensitivity of the model results to geometrical properties and facilitates inverting for them. We demonstrate the method with two examples. A salt diapir in an extending regime and free subduction of an oceanic plate under a continent. In both cases, small differences in the model's initial geometry lead to vastly different results. Be it the formation of faults or the velocity of plates. Using the salt diapir example, we demonstrate that, given an additional geophysical observation, we are able to invert for uncertain geometric properties. This highlights that geodynamic studies should investigate the sensitivity of their models to the initial geometry and include it in their inversion framework. We make our method available as part of the open‐source software geomIO.

How Aseismic Ridges Modify the Dynamics of Free Subduction: A 3-D Numerical Investigation

The subduction of positively buoyant features has been implicated in the development of flat and shallow dipping slabs, the formation of cusps in trench geometry, and the cessation of associated arc magmatism. However, how such buoyant anomalies influence subduction dynamics to produce these different tectonic expressions remains debated. In this paper, using a series of multi-material 3-D simulations of free subduction, we investigate how linear buoyant ridges modify subduction dynamics, in particular downgoing plate velocities, trench motions and slab morphology. We examine the sensitivity of results to downgoing plate age (affecting buoyancy and strength), ridge buoyancy and ridge location along the trench, finding that buoyant ridges can locally change slab sinking and trench retreat rates, in turn modifying the evolution of slab morphology at depth and trench shape at the surface. In all cases examined, trench retreat is reduced, or switches to trench advance, where the ridge subducts. These effects depend strongly on downgoing plate age: on young, weak plates, the change in trench shape is more localised than on old, strong plates. Slab shallowing at the ridge only occurs for young plates, while the stronger and more negatively buoyant older plates pull down the ridge at a steeper angle than the rest of the slab. On old plates, ridges located near regions of trench stagnation or advance, which typically develop in wide slabs, have a stronger effect on trench and slab shape. The combined effects of buoyant feature location, subducting plate age and overriding plate properties can result in a range of responses: from mainly trench deformation, through local slab shallowing, to the formation of a flat slab, a variation in expressions also observed on Earth.

Fluid mechanics of free subduction on a sphere. Part 1. The axisymmetric case

To understand how a spherical geometry influences the dynamics of gravity-driven subduction of the oceanic lithosphere on Earth, we study a simple model of a thin and dense axisymmetric shell of thickness $h$ and viscosity $\eta _1$ sinking in a spherical body of fluid with radius $R_0$ and a lower viscosity $\eta _0$ . Using scaling analysis based on thin viscous shell theory, we identify a fundamental length scale, the ‘bending length’ $l_b$ , and two key dimensionless parameters that control the dynamics: the ‘flexural stiffness’ $St = (\eta _1/\eta _0)(h/l_b)^3$ and the ‘sphericity number’ $\varSigma = (l_b/R_0)\cot \theta _t$ , where $\theta _t$ is the angular radius of the subduction trench. To validate the scaling analysis, we obtain a suite of instantaneous numerical solutions using a boundary-element method based on new analytical point-force Green functions that satisfy free-slip boundary conditions on the sphere's surface. To isolate the effect of sphericity, we calculate the radial sinking speed $V$ and the hoop stress resultant $T_2$ at the leading end of the subducted part of the shell, both normalised by their ‘flat-Earth’ values (i.e. for $\varSigma = 0$ ). For reasonable terrestrial values of $\eta _1/\eta _0$ ( $\approx$ several hundred), sphericity has a modest effect on $V$ , which is reduced by $< 7\,\%$ for large plates such as the Pacific plate and by up to 34 % for smaller plates such as the Cocos and Philippine Sea plates. However, sphericity has a much greater effect on $T_2$ , increasing it by up to 64 % for large plates and 240 % for small plates. This result has important implications for the growth of longitudinal buckling instabilities in subducting spherical shells.

Effect of Buoyant Sediment Overlying Subducting Plates on Trench Geometry: 3D Viscoelastic Free Subduction Modeling

AbstractBuoyant trench sediment generates an isostatic restoring force that opposes the slab pull of subducting plates; however, the influence of this force on subduction dynamics is poorly understood. Here, we performed three‐dimensional free subduction simulations adopting a variety of sediment distributions along a trench to investigate the correlation between trench motions and heterogeneous buoyant forces. Two endmembers, sediment‐rich and sediment‐starved centers, induced convex and concave trenches, respectively. The trench curvatures obtained from the natural subduction zones are well constrained by our models over wide magnitudes (1–9 km) and wavelengths (100–400 km) of deficient and excess sediments. Conversely, a uniform sediment distribution leads to an extremely narrow range of trench curvatures. These results imply that trench sediment contributes significantly to deviated curvatures along trench strikes. Finally, we observed stress localization at shallow slab hinges that lie beneath abundant sediments, clarifying the relationship between sediment thickness and subduction earthquakes.

Stress and deformation mechanisms at a subduction zone: insights from 2-D thermomechanical numerical modelling

SUMMARY Numerous processes such as metamorphic reactions, fluid and melt transfer and earthquakes occur at a subducting zone, but are still incompletely understood. These processes are affected, or even controlled, by the magnitude and distribution of stress and deformation mechanism. To eventually understand subduction zone processes, we quantify here stresses and deformation mechanisms in and around a subducting lithosphere, surrounded by asthenosphere and overlain by an overriding plate. We use 2-D thermomechanical numerical simulations based on the finite difference and marker-in-cell method and consider a 3200 km wide and 660 km deep numerical domain with a resolution of 1 km by 1 km. We apply a combined visco-elasto-plastic deformation behaviour using a linear combination of diffusion creep, dislocation creep and Peierls creep for the viscous deformation. We consider two end-member subduction scenarios: forced and free subduction. In the forced scenario, horizontal velocities are applied to the lateral boundaries of the plates during the entire simulation. In the free scenario, we set the horizontal boundary velocities to zero once the subducted slab is long enough to generate a slab pull force large enough to maintain subduction without horizontal boundary velocities. A slab pull of at least 1.8 TN m–1 is required to continue subduction in the free scenario. We also quantify along-profile variations of gravitational potential energy (GPE). We evaluate the contributions of topography and density variations to GPE variations across a subduction system. The GPE variations indicate large-scale horizontal compressive forces around the trench region and extension forces on both sides of the trench region. Corresponding vertically averaged differential stresses are between 120 and 170 MPa. Furthermore, we calculate the distribution of the dominant deformation mechanisms. Elastoplastic deformation is the dominant mechanism in the upper region of the lithosphere and subducting slab (from ca. 5 to 60 km depth from the top of the slab). Viscous deformation dominates in the lower region of the lithosphere and in the asthenosphere. Considering elasticity in the calculations has an important impact on the magnitude and distribution of deviatoric stress; hence, simulations with increased shear modulus, in order to reduce elasticity, exhibit considerably different stress fields. Limiting absolute stress magnitudes by decreasing the internal friction angle causes slab detachment so that slab pull cannot be transmitted anymore to the horizontal lithosphere. Applying different boundary conditions shows that forced subduction simulations are stronger affected by the applied boundary conditions than free subduction simulations. We also compare our modelled topography and gravity anomaly with natural data of seafloor bathymetry and free-air gravity anomalies across the Mariana trench. Elasticity and deviatoric stress magnitudes of several hundreds of MPa are required to best fit the natural data. This agreement suggests that the modelled flexural behaviour and density field are compatible with natural data. Moreover, we discuss potential applications of our results to the depth of faulting in a subducting plate and to the generation of petit-spot volcanoes.

Subduction Initiation During Collision‐Induced Subduction Transference: Numerical Modeling and Implications for the Tethyan Evolution

AbstractThe collision‐induced subduction transference is a composite dynamic process including both the terrane collision/accretion and the subduction initiation (SI) at the neighboring passive margin. This process occurred repeatedly during the evolution of Tethyan systems, with multiple ribbon‐like continents or microcontinents drifting from Gondwana in the southern hemisphere and accreting to the Eurasian continent since Paleozoic. In the previous numerical studies, the dynamics of terrane collision and induced SI are investigated individually, which however need to be integrated to study the controlling factors and time scales of collision‐induced subduction transference. Systematic numerical models are conducted with variable properties of converging plates and different boundary conditions. The model results indicate that the forced convergence, rather than pure free subduction, is required to trigger and sustain the SI at the neighboring passive margin after terrane collision. In addition, a weak passive margin can significantly promote the occurrence of SI, by decreasing the required boundary force to reasonable value of plate tectonics. The lengths of subducted oceanic slab and accreting terrane play secondary roles in the occurrence of SI after collision. Under the favorable conditions of collision‐induced subduction transference, the time required for SI after collision is generally short within 10 Myr, which is consistent with the general geological records of Cimmerian collision and the following Neo‐Tethyan SI. In contrast, the stable Indian passive margin and absence of SI in the present Indian Ocean may be due to the low convergent force and/or the lack of proper weak zones, which remains an open question.

Subduction of non-Newtonian plates: Thin-sheet dynamics of slab necking and breakoff

Summary We use a hybrid boundary-integral/thin-sheet (‘BITS’) model to investigate the subduction of 2-D viscous sheets with composite diffusion creep/dislocation creep rheology, with a focus on the conditions required for slab necking and breakoff. To validate the model, we show that its predictions of the sinking speed of the slab follow a universal scaling law identical to one previously derived for purely Newtonian sheets. We obtain analytical expressions for the fiber stress resultant and bending moment of the sheet during deformation by pure stretching and pure bending, respectively, and show that the non-Newtonian dislocation creep component of the rheology significantly weakens the sheet for both types of deformation. We solve the BITS equations for two distinct situations: “free” subduction, in which the slab pulls an attached negatively buoyant plate without hindrance; and “arrested” subduction that slows or stops when a positively buoyant (continental) portion of the attached plate arrives at the trench. Strong lithospheric thinning is difficult to obtain in free subduction and requires dominantly non-Newtonian rheology, i.e. large (>5) values of the characteristic ratio λ of the Newtonian to the non-Newtonian viscosity. However, during arrested subduction strong thinning leading to breakoff occurs for much lower values of λ, and the point of maximum thinning is at shallower depths. These results are explained by a one-dimensional viscous dripping model (‘DRIP’) of a vertical slab with a composite rheology and an arbitrary kinematically prescribed time-dependent convergence rate U0(t). By injecting into DRIP convergence histories predicted by BITS, we find that DRIP reproduces closely the results of the more complicated BITS model. The DRIP model shows that the convergence rate controls slab breakoff in two distinct ways. On the one hand, the rate of lithospheric thinning is proportional to the accumulated convergence, i.e., the time integral of U0(t). The breakoff depth, on the other hand, is controlled by the convergence rate history, being shallower when U0 rapidly decreases (arrested subduction) after an initial period of oceanic subduction.

콤솔 멀티피직스를 이용한 점탄성 자유 섭입 수치 모형에 관한 벤치마크 연구: 2차원 라그랑지안 유한요소법 틀 내에서

콤솔 멀티피직스를 이용한 점탄성 자유 섭입 수치 모형에 관한 벤치마크 연구: 2차원 라그랑지안 유한요소법 틀 내에서

The effect of rheological approximations in 3-D numerical simulations of subduction and collision

Subduction and collision zones evolve differently from one another due to different rheological properties, different amounts of regional isostatic compensation, and the different mechanisms by which forces are applied to the convergent plates. The rheology of mantle and lithosphere is known to have the largest influence on the dynamics of subduction and continental collision. However, previous 3-D geodynamic models of subduction/collision processes have used various rheological approximations, making their results difficult to compare, since there is no clear understanding on the extent of these approximations on the dynamics. Here, we test the effect of rheological approximations in a geometrically simplified 3-D model in which the dynamics of subduction and collision is entirely driven by slab-pull (i.e., “free subduction”). The results exhibit a wide range of behaviors depending on the rheological law employed: from linear viscous (LV) to temperature-dependent visco-elasto-plastic (VEP) rheology that takes into account both diffusion and dislocation creep. Rheological approximations exert a strong effect on slab morphology, trench migration rates and topography in the upper plate. For example, we find that slab dynamics varies drastically between end member models. In LV approximations, slabs do not break-off, but instead buckle at the transition zone. The dynamics is dominated by low strain rates and diffuse stresses, and topography is controlled by lateral extrusion of material which homogenizes amplitudes in the upper plate. In simulations with non-linear VEP rheology, slab break-off occurs at the corners of the indenting continent, opening of a back-arc basin is fast and topography has a laterally heterogenous distribution of amplitudes. Opening of a back-arc basin is observed only in simulations with non-linear VEP rheology. When the rheology is approximated with multiple discrete layers, the dynamics also resembles more the non-linear VEP case.

Control of slab width on subduction‐induced upper mantle flow and associated upwellings: Insights from analog models

AbstractThe impact of slab width W (i.e., trench‐parallel extent) on subduction‐induced upper mantle flow remains uncertain. We present a series of free subduction analog models where W was systematically varied to upscaled values of 250–3600 km to investigate its effect on subducting plate kinematics and upper mantle return flow around the lateral slab edges. We particularly focused on the upwelling component of mantle flow, which might promote decompression melting and could thereby produce intraplate volcanism. The models show that W has a strong control on trench curvature and on the trench retreat, subducting plate, and subduction velocities, generally in good agreement with previous modeling studies. Upper mantle flow velocity maps produced by means of a stereoscopic particle image velocimetry technique indicate that the magnitude of the subduction‐induced mantle flow around the lateral slab edges correlates positively with the product of W and trench retreat velocity. For all models an important upwelling component is always produced close to the lateral slab edges, with higher magnitudes for wider slabs. The trench‐parallel lateral extent of this upwelling component is the same irrespective of W, but its maximum magnitude gets located closer to the subducting plate in the trench‐normal direction and it is more focused when W increases. For W ≤ 2000 km the upwelling occurs laterally (in the trench‐parallel direction) next to the subslab domain and the mantle wedge domain, while for W ≥ 2000 km it is located only next to the subslab domain and focuses closer to the trench tip, because of stronger poloidal flow in the mantle wedge extending laterally.

The role of viscoelasticity in subducting plates

AbstractSubduction of tectonic plates into Earth's mantle occurs when one plate bends beneath another at convergent plate boundaries. The characteristic time of deformation at these convergent boundaries approximates the Maxwell relaxation time for olivine at lithospheric temperatures and pressures, it is therefore by definition a viscoelastic process. While this is widely acknowledged, the large‐scale features of subduction can, and have been, successfully reproduced assuming the plate deforms by a viscous mechanism alone. However, the energy rates and stress profile within convergent margins are influenced by viscoelastic deformation. In this study, viscoelastic stresses have been systematically introduced into numerical models of free subduction, using both the viscosity and shear modulus to control the Maxwell relaxation time. The introduction of an elastic deformation mechanism into subduction models produces deviations in both the stress profile and energy rates within the subduction hinge when compared to viscous only models. These variations result in an apparent viscosity that is variable throughout the length of the plate, decreasing upon approach and increasing upon leaving the hinge. At realistic Earth parameters, we show that viscoelastic stresses have a minor effect on morphology yet are less dissipative at depth and result in an energy transfer between the energy stored during bending and the energy released during unbending. We conclude that elasticity is important during both bending and unbending within the slab hinge with the resulting stress loading and energy profile indicating that slabs maintain larger deformation rates at smaller stresses during bending and retain their strength during unbending at depth.

Evolution of 3-D subduction-induced mantle flow around lateral slab edges in analogue models of free subduction analysed by stereoscopic particle image velocimetry technique

We present analogue models of free subduction in which we investigate the three-dimensional (3-D) subduction-induced mantle flow focusing around the slab edges. We use a stereoscopic Particle Image Velocimetry (sPIV) technique to map the 3-D mantle flow on 4 vertical cross-sections for one experiment and on 3 horizontal depth-sections for another experiment. On each section the in-plane components are mapped as well as the out-of-plane component for several experimental times. The results indicate that four types of maximum upwelling are produced by the subduction-induced mantle flow. The first two are associated with the poloidal circulation occurring in the mantle wedge and in the sub-slab domain. A third type is produced by horizontal motion and deformation of the frontal part of the slab lying on the 660 km discontinuity. The fourth type results from quasi-toroidal return flow around the lateral slab edges, which produces a maximum upwelling located slightly laterally away from the sub-slab domain and can have another maximum upwelling located laterally away from the mantle wedge. These upwellings occur during the whole subduction process. In contrast, the poloidal circulation in the mantle wedge produces a zone of upwelling that is vigorous during the free falling phase of the slab sinking but that decreases in intensity when reaching the steady-state phase. The position of the maximum upward component and horizontal components of the mantle flow velocity field has been tracked through time. Their time-evolving magnitude is well correlated to the trench retreat rate. The maximum upwelling velocity located laterally away from the subducting plate is ∼18–24% of the trench retreat rate during the steady-state subduction phase. It is observed in the mid upper mantle but upwellings are produced throughout the whole upper mantle thickness, potentially promoting decompression melting. It could thereby provide a source for intraplate volcanism, such as Mount Etna in the Mediterranean, the Chiveluch group of volcanoes in Kamchatka and the Samoan hotspot near Tonga.

Dynamical similarity and density (non-) proportionality in experimental tectonics

Abstract The utility of analog laboratory models for tectonic processes relies on their dynamical similarity to their natural prototypes. Dynamical similarity is often thought to require that the density distribution in the model be a constant (position-independent) multiple of that in the prototype, a principle due to Hubbert (1937) . To clarify the status of this rule, we nondimensionalize the equations and boundary/initial conditions governing simple models of three paradigmatic processes: gravity tectonics, compressional tectonics, and free subduction. The results show that density proportionality, while compatible with dynamical similarity, is not always required by it in systems with negligible inertia, a category that includes most geological and tectonic processes. The density proportionality rule is therefore unnecessarily restrictive in many cases, implying that the range of analog materials that can be used to construct properly scaled models is wider than commonly recognized.

Effect of aseismic ridge subduction on slab geometry and overriding plate deformation: Insights from analogue modeling

We present analogue models simulating the subduction of a buoyant ridge oriented perpendicularly or obliquely with respect to the trench, beneath an advancing overriding plate. The convergence velocity is imposed by lateral boundary conditions in this experimental set. We analyze the three-dimensional geometry of the slab, the deformation and topography of the overriding plate. Experiments suggest that ridge subduction diminishes the dip of the slab, eventually leading to the appearance of a horizontal slab segment in case boundary conditions impose a rapid convergence. This result contrasts with that obtained in free subduction experiments, in which ridge subduction diminishes the convergence velocity which, in turn, increases the dip of the slab beneath the ridge. The slab dip decrease is accompanied by the indentation of the overriding plate by the ridge, resulting in arc curvature. Experiments suggest that indentation is larger for small convergence velocity and large slab dip. Ridge subduction also uplifts the overriding plate. Uplift first occurs close to the trench (~fore-arc area) and is accompanied by the flexural subsidence of the overriding plate behind the uplifted area (~back-arc subsidence). The uplifted area migrates within the overriding plate interiors following the appearance of a horizontal slab segment. These results are compared with natural examples of ridge subduction in the circum-Pacific area. They explain why ridge subduction may have contrasted effects on the overriding plate dynamics depending on the global conditions that constrain the converging system.

Numerical models of slab migration in continental collision zones

Abstract. Continental collision is an intrinsic feature of plate tectonics. The closure of an oceanic basin leads to the onset of subduction of buoyant continental material, which slows down and eventually stops the subduction process. In natural cases, evidence of advancing margins has been recognized in continental collision zones such as India-Eurasia and Arabia-Eurasia. We perform a parametric study of the geometrical and rheological influence on subduction dynamics during the subduction of continental lithosphere. In our 2-D numerical models of a free subduction system with temperature and stress-dependent rheology, the trench and the overriding plate move self-consistently as a function of the dynamics of the system (i.e. no external forces are imposed). This setup enables to study how continental subduction influences the trench migration. We found that in all models the slab starts to advance once the continent enters the subduction zone and continues to migrate until few million years after the ultimate slab detachment. Our results support the idea that the advancing mode is favoured and, in part, provided by the intrinsic force balance of continental collision. We suggest that the advance is first induced by the locking of the subduction zone and the subsequent steepening of the slab, and next by the sinking of the deepest oceanic part of the slab, during stretching and break-off of the slab. These processes are responsible for the migration of the subduction zone by triggering small-scale convection cells in the mantle that, in turn, drag the plates. The amount of advance ranges from 40 to 220 km and depends on the dip angle of the slab before the onset of collision.

Dynamics of free subduction from 3‐D boundary element modeling

In order better to understand the physical mechanisms underlying free subduction, we perform three‐dimensional boundary‐element numerical simulations of a dense fluid sheet with thickness h and viscosity η2 sinking in an ‘ambient mantle’ with viscosity η1. The mantle layer is bounded above by a traction‐free surface, and is either (1) infinitely deep or (2) underlain by a rigid boundary at a finite depth H + d, similar to the typical geometry used in laboratory experiments. Instantaneous solutions in configuration (1) show that the sheet's dimensionless ‘stiffness’ S determines whether the slab's sinking speed is controlled by the viscosity of the ambient mantle (S < 1) or the viscosity of the sheet itself (S > 10). Time‐dependent solutions with tracers in configuration (2) demonstrate a partial return flow around the leading edge of a retreating slab and return flow around its sides. The extra ‘edge drag’ exerted by the flow around the sides causes transverse deformation of the slab, and makes the sinking speed of a 3‐D slab up to 40% less than that of a 2‐D slab. A systematic investigation of the slab's interaction with the bottom boundary as a function of η2/η1 and H/h delineates a rich regime diagram of different subduction modes (trench retreating, slab folding, trench advancing) and reveals a new ‘advancing‐folding’ mode in which slab folding is preceded by advancing trench motion. The solutions demonstrate that mode selection is controlled by the dip of the leading edge of the slab at the time when it first encounters the bottom boundary.