Gabriele Morra

Cenozoic intraplate volcanism in central Asia is characterized by individual small volcanic provinces spread over three broad regions, thousand of kilometers long and hundred of kilometers wide, whose origin remains controversial. Seismological observation of high‐velocity anomalies at the base of the upper mantle, where the subducted Pacific plate has been sliding for tens of millions of years, prompted the “Big Mantle Wedge” hypoth- esis that some of the surface volcanism emerges from the rise of partially molten mantle rocks originating above the slab. Constraints from mineral physics, petrology, and geophysics have suggested that the wadsleyite layer at the top of the transition zone can store large amounts of volatiles. We show that while a 1‐km single diapir would rise too slowly and diffuse to the surrounding mantle, a sufficiently large group of diapirs would reach the surface because its rising speed increases with the square of the cluster size. We relate the regular distribution of the sur- face volcanism, characterized by cones that are about 200‐km distant, to the thickness of the volatile‐rich layer from which diapirs rise. Doing so we find that a 20‐km thick layer of diapirs on top of the transition zone would explain most of the observed volcanism. Our models suggest that decarbonation‐dehydration melting in the tran- sition zone is the main cause of the upwellings and of the surface volcanism in central Asia.

The present tessellation of the Earth's surface into tectonic plates displays a remarkably regular plate size distribution, described by either one (Sornette and Pisarenko, 2003) or two (Bird, 2003) statistically distinct groups, characterised by large and small plate size. A unique distribution implies a hierarchical structure from the largest to the smallest plate. Alternatively, two distributions indicate distinct evolutionary laws for large and small plates, the first tied to mantle flow, the second determined by a hierarchical fragmentation process. We analyse detailed reconstructions of plate boundaries during the last 200 Myr and find that (i) large and small plates display distinct statistical distributions, (ii) the small plates display little organisational change since 60 Ma and (iii) the large plates oscillate between heterogeneous (200–170 Myr and 65–50 Ma) and homogeneous (120–100 Ma) plate tessellations on a timescale of about 100 Myr. Heterogeneous states are reached more rapidly, while the plate configura- tion decays into homogeneous states following a slower asymptotic curve, suggesting that hetero- geneous configurations are excited states while homogeneous tessellations are equilibrium states. We explain this evolution by proposing a model that alternates between bottom- and top-driven Earth dynamics, physically described by fluid-dynamic analogies, the Rayleigh–Benard and Bénard–Marangoni convection, respectively. We discuss the implications for true polar wander (TPW), global kinematic reorganisations (50 and 100 Ma) and the Earth's magnetic field inversion frequency.

A proper determination of the lower-mantle viscosity profile is fundamental to understanding Earth geodynamics. Based on results coming from different sources, several models have been proposed to constrain the variations of viscosity as a function of pressure, stress and temperature. While some models have proposed a relatively modest viscosity variation across the lower mantle, others have proposed variations of several orders of magnitude. Here, we have determined the viscosity of ferropericlase, a major mantle mineral, and explored the role of the iron high-to-low spin transition. Viscosity was described within the elastic strain energy model, in which the activation parameters are obtained from the bulk and shear wave velocities. Those velocities were computed combining first principles total energy calculations and the quasi-harmonic approximation. As a result of a strong elasticity softening across the spin transition, there is a large reduction in the activation free energies of the materials creep properties, leading to viscosity undulations. These results suggest that the variations of the viscosity across the lower mantle, resulting from geoid inversion and postglacial rebound studies, may be caused by the iron spin transition in mantle minerals. Implications of the undulated lower mantle viscosity profile exist for both, down- and up-wellings in the mantle. We find that a viscosity profile characterized by an activation free energy of G∗(z0)∼300–400 kJ/mol based on diffusion creep and dilation factor δ = 0.5 better fits the observed high velocity layer at mid mantle depths, which can be explained by the stagnation and mixing of mantle material. Our model also accounts for the growth of mantle plume heads up to the size necessary to explain the Large Igneous Provinces that characterize the start of most plume tracks.

Large tectonic plates are known to be susceptible to internal deformation, leading to a range of phenomena including intraplate volcanism. However, the space and time dependence of intraplate deformation and its relationship with changing plate boundary configurations, subducting slab geometries, and absolute plate motion is poorly under- stood. We utilise a buoyancy driven Stokes flow solver, BEM-Earth, to investigate the contribution of subducting slabs through time on Pacific Plate motion and plate-scale deformation, and how this is linked to intraplate volcanism. We produce a series of geodynamic models from 62 to 42 Ma in which the plates are driven by the attached
subducting slabs and mantle drag/suction forces. We compare our modelled intraplate deformation history with those types of intraplate volcanism that lack a clear age pro- gression. Our models suggest that changes in Cenozoic subduction zone topology caused intraplate deformation to trigger volcanism along several linear seafloor struc- tures, mostly by reactivation of existing seamount chains, but occasionally creating new volcanic chains on crust weakened by fracture zones and extinct ridges. Around 55 Ma subduction of the Pacific-Izanagi ridge reconfigured the major tectonic forces acting on the plate by replacing ridge push with slab pull along its north-western perimeter, causing lithospheric extension along pre-existing weaknesses. Large scale deforma- tion observed in the models coincides with the seamount chains of Hawaii, Louisville, Tokelau, and Gilbert during our modelled time period of 62 to 42 Ma. We suggest that extensional stresses between 72 and 52 Ma are the likely cause of large parts of the formation of the Gilbert chain and that localised extension between 62 and 42 Ma could cause late-stage volcanism along the Musicians Volcanic Ridges. Our models demon- strate that early Cenozoic changes in Pacific plate driving forces only cause relatively minor changes in Pacific absolute plate motions, and cannot be responsible for the Hawaii-Emperor Bend (HEB), confirming previous interpretations that the 47 Ma HEB does not reflect an absolute plate motion event.

Tsunamigenic megathrust earthquakes, like the 2004 Sumatra-Andaman and 2011 Tohoku events, are the most dramatic consequences of subduction dynamics. The classical view is that megathrusts release elastic energy due to the rupture of a fault with a width of tens of kilometers in the down-dip direction and a length of hundreds to a thousand kilometers along the trench. However, recent research, particularly work on the Tohoku event, has suggested that the generation of huge tsunamis may require the release of gravitational energy as well as elastic energy [George et al., 2011]. Our growing understanding of the role of gravitational energy in generating tsunamis following megathrust earthquakes points to the need to reevaluate earthquake and tsunami hazard assessments.

We investigate the use of the Multipole-accelerated Boundary Element Method (BEM) and of the Singularity Method for studying the interaction of many bubbles rising in a volcanic conduit. Observation shows that the expression of volcanic eruption is extremely variable, from slow release of magma to catastrophic explosive manifestation. We investigate the application of the Fast Multipole Method to the solution of (i) the Boundary Element Formulation of the Stokes flow and of (ii) the particle formulation using the Stokeslets, the Green Function of the Stokes flow law, as a particle kernel. We show how these implementations allow for the first time to numerically model in a dynamic setting a very large number of bubbles, i.e few thousands with the BEM models, allowing investigating the feedback between the single bubble deformation and their collective evolution, and few hundred of thousands of bubbles with the particle approach. We illustrate how this method can be used to investigate the intense interaction of a large number of bubbles and suggest a framework for studying the feedback between many bubbles and a complex thermal nonlinear magmatic

In the last two decades it has been proposed several times that a non-monotonic profile might fit the average lower mantle radial viscosity. Most proposed profiles consist in a more or less broad viscosity hill in the middle of the mantle, at a depth roughly between 1200 km and 2000 km. Also many tomographic models display strong signals of the presence of “fast” material lying at mid mantle depths and a recent spectral analysis of seismic tomography shows a very clear transition for degree up to around 16 at a less than 1500 km depth. Finally latest works, both theoretical and experimental, on the high-to-low- spin transition for periclase, have suggested that the high-spin to low-spin transition of Fe++ might lie at the heart of all these observations. To verify the dynamical compatibility between possible mantle profile and observed tomographic images and compare them with possible mineral physics scenarios, such as the spin transition, we employ here a recently developed Fast Multipole-accelerated Boundary Element Method (FMM-BEM), a numerical approach for solving the viscous momentum equation in a global spherical setting, for simulating the interaction of an individual slab with a mid mantle smooth discontinuity in density and viscosity. We have focused on the complexities induced to the behaviour of average and very large plates O (2000–10,000 km), characteristic of the Farallon, Tethys and Pacific plate subducting during the Cenozoic, demonstrating that the a mid mantle density and/or viscosity discontinuity produces a strong alteration of the sinking velocity and an intricate set of slab morphologies. We also employ the Kula–Farallon plate system subducting at 60 Ma as a paradigmatic case, which reveals the best high resolution tomography models and clearly suggests an interaction with a strong and/or denser layer in the mantle. Our 38 models show that a plate might or might not penetrate into the lowest mantle and might stall in the mid lower mantle for long periods, depending on the radial profiles of density and viscosity, within a realistic range (viscosity 1, 10 or 100 times more viscous of the rest of the mantle, and a change of differential density in the range −2% to 2%), of a transitional layer of 200 km or 500 km. We conclude that a layer with high viscosity or negative density would naturally trigger the observed geodynamic snapshot. We finally propose a scenario in which the long time accumulation of depleted slabs in the mid mantle would give rise to a partially chemically stratified mantle, starting from the less prominent high-spin to low-spin contribution on the basis of mantle density and rheology.

We use the Multipole–Boundary Element Method (MP-BEM) to simulate regional and global geody- namics in a spherical 3-D setting. We first simulate an isolated subducting rectangular plate with length (Llitho) and width (Wlitho) varying between 0.5 and 2 times the radius of the Earth (REarth) and with viscosity hlitho varying between 100 and 500 times the upper mantle (hUM), sinking in a layered mantle characterized by lower-upper mantle viscosity ratio l = hLM/hUM varying between 1 and 80. In a mantle with small upper/ lower viscosity contrast (l ≅ 1), trench and plate motions are weakly dependent on Wlitho; plate motion is controlled by slab pull if Llitho ≤ REarth, while for longer plates plate speed strongly decreases because of the plate basal friction and flow reorganization. An increasing viscosity ratio l gradually breaks this pattern, and for l ≅ 10 combined with Wlitho ≈ REarth (and greater) trench advance and retreat are simultaneously observed. These results offer a first-order explanation of the origin of the size (Llitho ≈ Wlitho ≈ REarth) of the largest plates observed over the past 150 Myr. Finally, two global plate tectonic simulations are per- formed from reconstructed plates and slabs at 25 Ma before present and before 100 Ma, respectively. It is shown that MP-BEM predicts present plate kinematics if plate-mantle decoupling is adopted for the lon- gest plates (Llitho > REarth). Models for 100 Ma show that the slab-slab interaction between India and Izanagi plates at 100 Ma can explain the propagation of the plate reorganization from the Indian to the Pacific plate.

We present a novel dynamic approach for solid–fluid coupling by joining two different numerical methods: the boundary-element method (BEM) and the finite element method (FEM). The FEM results describe the thermomechanical evolution of the solid while the fluid is solved with the BEM. The bidirectional feedback between the two domains evolves along a Lagrangian interface where the FEM domain is embedded inside the BEM domain. The feedback between the two codes is based on the calculation of a specific drag tensor for each boundary on finite element. The approach is presented here to solve the complex problem of the descent of a cold subducting oceanic plate into a hot fluid-like mantle. The coupling technique is shown to maintain the proper energy dissipation caused by the important secondary induced mantle flow induced by the lateral migrating of the subducting plate. We show how the method can be successfully applied for modelling the feedback between deformation of the oceanic plate and the induced mantle flow. We find that the mantle flow drag is singular at the edge of the retreating plate causing a distinct hook shape. In nature, such hooks can be observed at the northern end of the Tonga trench and at the southern perimeter, of the South American trench.

Oceanic arcs and deep-sea trenches are the surface expressions of oceanic plates subducting into the Earth’s mantle. We use a new numerical technique for simulating the dynamical evolution of the lithosphere-mantle  interaction in order to assess the causes of arc curvature. We group the possible causes into two classes, external feedback between the migrating lithosphere and the secondary induced mantle flow, and internal heterogeneities within the lithosphere, e.g., owing to differences in cooling ages of the plate at the trench. We statistically assess that almost all arcs on the Earth can be described by these hypotheses. The method is also directly applied to the Tonga and Aleutian arcs, bringing new insights on the origin of their shapes.

It has come to our attention that the constitutive relationship used in the modeling of geodynamical flow problems with strongly variable physical properties, should have additional terms in the stress tensor, known in the literature as Korteweg stresses (K-stresses). These stresses arising at diffuse interfaces, which can best be explained in terms of density gradients, have already been mentioned in the literature for more than one hundred years, but have not received attention recently until a combination of experimental and numerical evidence have confirmed their existence. We will discuss the important potential role these new terms have for geophysics. It has many ramifications in geodynamics, ranging from mantle convection to earthquakes and magma fragmentation.

We present a novel approach for modeling subduction using a Multipole-accelerated Boundary Element Method (BEM). The present approach allows large-scale modeling with a reduced number of elements and scales linearly with the problem size. For the first time the BEM has been applied to a subduction model in a spherical planet with an upper-lower mantle discontinuity, in conjunction with a free-surface mesh algorithm.

We present the first application in geodynamics of a (Fast Multipole) Accelerated Boundary Element Method (Accelerated-BEM) for Stokes flow. The approach offers the advantages of a reduced number of computa- tional elements and linear scaling with the problem size. We show that this numerical method can be fruitfully applied for the simulation of several geodynamic systems at the planetary scale in spherical coordinates, and we suggest a general approach for modeling combined mantle convection and plate tectonics. The first part of the paper is devoted to the technical exposition of the new approach, while the second part focuses on the effect played by Earth curvature on the subduction of a very wide oceanic lithosphere (W = 6,000 km and W = 9,000 km), comparing the effects of two different planetary radii (ER = 6,371 km, 2ER = 2 9 6,371 km), corresponding to an ‘‘Earth-like’’ model (ER) and to a ‘‘flat Earth’’ one (2ER). The results show a distinct difference between the two models: while the slab on a ‘‘flat Earth’’ shows a slight undulation, the same subducting plate on the ‘‘Earth-like’’ setting presents a dual behavior characterized by concave curvature at the edges and by a folding with wavelength of the order of magnitude of 1,000 km at the center of the slab.

Abstract This PhD work is partof a largerproject called OrogenyI, thoughtto be the first step towards obtaining a self-consistentdynamic model of the formation of the Alps. The chainoriginated by Continental convergence, partof a dynamicprocessin which the subductionin the Mediter¬ ranean playeda major role. The aim of this specificPhD projectis to providea numerical tool forforward dynamicmodeling, to be used to simulate the Mediterraneantectonic scenarios.

We present a general framework to generate time-dependent global subduction history models from kinematic plate reconstructions and explore their associated coupled plate–mantle dynamic behaviour. Slabs are constructed by advecting material into the mantle by prescribing its radial velocity and following the absolute tangential motion of the subducting plate. A simple geodynamic scenario where plates and slabs define isopycnic and isoviscous regions in an homogeneous or layered mantle was explored using the boundary element method-based software BEMEarth. The resulting dynamic behaviour was used to predict the absolute plate motion directions for the present day and a particular mid-cretaceous (125 Ma) kinematic model. We show how the methodology can be used to compare and revise kinematic reconstructions based on their effect on the balance of plate driving forces and the resulting Euler poles of subducting plates. As an example we compare the Farallon plate dynamics at 125 Ma in a global model with two reconstructions in the context of the evolution of the Western North American Cordillera. Our results suggest a method to identify episodes of absolute plate motions that are inconsistent with the expected plate dynamics.► We develop a framework to generate time-dependent global subduction history models. ► Global convection simulations based on model input used to predict plate velocity. ► Effect of mantle layering in plate velocity fitting is explored. ► Simulated mid-cretaceous motions fit observations better than reconstructed ones. ► Methodology to constrain absolute plate motions is proposed.

Abstract We prepose here a formulation for studying a generic geothermal flow based on the use of the Fast Multipole Boundary Element Method for potential flow for extrapolating the non-gaussian stochastic representation of the background hydraulic transmissivity. We propose to use the Karhunen-Love (KL) expansion for representing the transmissivity, as it holds the property of having a minimum representation entropy.

[1] Subduction dynamics is strongly dependent on the geometry and rheology of the subducting slab and adjacent plates, as well as on the induced mantle flow driven by the evolution of tectonic configurations along subduction zones. However, these processes, and the associated plate tectonic driving forces, are difficult to study using time-dependent 3-dimensional computer simulations due to limitations in computing resources. We investigate these phenomena with a novel numerical approach, using BEM-Earth, a Stokes flow solver based on the Boundary Element Method (BEM) with a Fast-Multipole (FM) implementation. The initial BEM-Earth model configurations self-consistently determine the evolution of the entire lithosphere-mantle system without imposing additional constraints in a whole-Earth spherical setting. We find that models without an overriding plate overestimate trench retreat by 65% in a 20 m.y. model run. Also, higher viscosity overriding plates are associated with higher velocity subducting slabs, analogue to faster oceanic plates subducting beneath more rigid continental lithosphere. In our models poloidal flows dominate the coupling between the down-going and overriding plates, with trench-orthogonal length variations in overriding plates inducing flows at least ∼2× stronger than trench-parallel width variations. However, deformation in the overriding plate is related to its length and width, with narrower and longer plates extending more than wider and shorter plates.

Subduction dynamics is strongly dependent on the geometry and rheology of the subducting slab and adjacent plates, as well as on the induced mantle flow driven by the evolution of tectonic configurations along subduction zones. However, these processes, and the associated plate tectonic driving forces, are difficult to study using time-dependent 3-dimensional computer simulations due to limitations in computing resources.

Given the wealth and complexity of current plate reconstructions, a simple yet general methodology to compare them in relation to their global scale predicted geodynamics could be of great help revising them to improve their fit with the present structure of the mantle and with geological observations. The main drivers of large-scale mantle flow and plate tectonics are considered to be the pull and suction forces associated to cold material sinking into the mantle [1].

We present a general framework to generate time-dependent global subduction history models from kinematic plate reconstructions and explore their associated coupled plate–mantle dynamic behaviour. Slabs are constructed by advecting material into the mantle by prescribing its radial velocity and following the absolute tangential motion of the subducting plate.

Lithospheric plates bend at subduction zones where the vertical motions of the slabs are converted to surface plate motions. To understand the mechanics of plate bending we derive scaling laws for the deflection at the margin, i.e. radius and dip, from numerical models of a subducting viscoelastic plate. In such dynamic system we find that the buoyancy and the stiffness of the plates control the radius and the dip, as well as the plate motions toward the trench. This mechanical model successfully predicts the curvature of published three- dimensional laboratory and numerical models. For a thorough comparison with the observable, we have also implemented forces additional to the slab pull, such as the suction force and far-field stresses. By increasing or resisting the torque applied at the trench by the slab, these forces can largely rearrange the dip and the radius of slabs and the inherent plate motions, although they do not alter the observed anticorrelation between radius and dip. Similar inverse correlation relationship and dip-radius ranges are shown by most of the subduction zones analysed from a global compilation. Radii in the range of 100–350 km and dips of 30°–70° for slabs that extends to the bottom of the upper mantle are compatible with the models, and allow estimating an average lith- ospheric viscosity contrast of ~200 in the bending with respect to the ambient mantle. Radii and dips outside of this range are in good agreement with the trends and the magnitudes of models that include suction and far-field forces. In all these subduction zones, the correlation between dip, radius and plate velocity is found to be com- patible with that of the models, showing how relevant bending is for the dynamics of Earth.

Crustal faults and sharp material transitions in the crust are usually represented as triangulated surfaces in structural geological models. The complex range of volumes separating such surfaces is typically three-dimensionally meshed in order to solve equations that describe crustal deformation with the finite-difference (FD) or finite-element (FEM) methods. We show here how the Boundary Element Method, combined with the Multipole approach, can revolutionise the calculation of stress and strain, solving the problem of computational scalability from reservoir to basin scales. The Fast Multipole Boundary Element Method (Fast BEM) tackles the difficulty of handling the intricate volume meshes and high resolution of crustal data that has put classical Finite 3D approaches in a performance crisis. The two main performance enhancements of this method: the reduction of required mesh elements from cubic to quadratic with linear size and linear-logarithmic runtime; achieve a reduction of memory and runtime requirements allowing the treatment of a new scale of geodynamic models. This approach was recently tested and applied in a series of papers by [1, 2, 3] for regional and global geodynamics, using KD trees for fast identification of near and far-field interacting elements, and MPI parallelised code on distributed memory architectures, and is now in active development for crustal dynamics. As the method is based on a free-surface, it allows easy data transfer to geological visualisation tools where only changes in boundaries and material properties are required as input parameters. In addition, easy volume mesh sampling of physical quantities enables direct integration with existing FD/FEM code.

The most spectacular example of a plate convergence event on Earth is the motion of the Indian plate towards Eurasia at speeds in excess of 18 cm yr−1 (ref. 1), and the subsequent collision. Continental buoyancy usually stalls subduction shortly after collision, as is seen in most sections of the Alpine–Himalayan chain. However, in the Indian section of this chain, plate velocities were merely reduced by a factor of about three when the Indian continental margin impinged on the Eurasian trench about 50 million years ago. Plate convergence, accompanied by Eurasian indentation, persisted throughout the Cenozoic era1–3, suggesting that the driving forces of convergence did not vanish on continental collision. Here we estimate the density of the Greater Indian continent, after its upper crust is scraped off at the Himalayan front, and find that the continental plate is readily subductable. Using numerical models, we show that subduction of such a dense continent reduces convergence by a factor similar to that observed. In addition, an imbalance between ridge push and slab pull can develop and cause trench advance and indentation. We conclude that the subduction of the dense Indian continental slab provides a significant driving force for the current India– Asia convergence and explains the documented evolution of plate velocities following continental collision.

It is well accepted that subduction of the cold lithosphere is a crucial component of the Earth’s plate tectonic style of mantle convection. But whether and how subducting plates penetrate into the lower mantle is the subject of continuing debate, which has substantial implications for the chemical and thermal evolution of the mantle1,2. Here we identify lower-mantle slab penetration events by comparing Cenozoic plate motions at the Earth’s main subduction zones3 with motions predicted by fully dynamic models of the upper-mantle phase of subduction, driven solely by downgoing plate density4. Whereas subduction of older, intrinsically denser, lithosphere occurs at rates consistent with the model, younger lithosphere (of ages less than about 60 Myr) often subducts up to two times faster, while trench motions are very low. We conclude that the most likely explanation is that older lithosphere, subducting under significant trench retreat, tends to lie down flat above the transition to the high-viscosity lower mantle, whereas younger lithosphere, which is less able to drive trench retreat and deforms more readily, buckles and thick- ens. Slab thickening enhances buoyancy (volume times density) and thereby Stokes sinking velocity, thus facilitating fast lower- mantle penetration. Such an interpretation is consistent with seis- mic images of the distribution of subducted material in upper and lower mantle5,6. Thus we identify a direct expression of time-dependent flow between the upper and lower mantle.

The bending strength of subducting lithosphere plays a critical role in the Earth’s plate tectonics and mantle convection, modulating the amount of slab pull transmitted to the surface and setting the boundary conditions under which plates move and deform. However, it is the subject of a lively debate how much of the potential energy of the downgoing plate is consumed in bending the plate and how the lithospheric strength is defined during this process. We model the subduction of a viscoelastic lithosphere, driven solely by the downgoing plate’s buoyancy, freely sinking in a passive mantle, represented by drag forces. To investigate the dynamics of bending, (1) we vary the density and the viscosity profile within the plate from isoviscous, where strength is distributed, to strongly layered, where strength is concentrated in a thin core, and (2) we map the stress, strain, and dissipation along the downgoing plate. The effective plate strength during bending is not a simple function of average plate viscosity but is affected by rheological layering and plate thinning. Earth-like layered plates allow for the transmission of large fractions of slab pull (75-80%) through the bend and yield a net slab pull of Fnet = 1 to 6 10^12 N m^-1, which varies with the SP buoyancy of plates. In all models, only a minor portion of the energy is dissipated in the bending. Surprisingly, bending dissipation hardly varies with lithospheric viscosity because in our dynamic system, the plates minimize overall dissipation rate by adjusting their bending curvature. As a result, bending dissipation, FB, is mainly controlled by the bending moment work rate exerted by slab pull. We propose a new analytical formulation that includes a viscosity-dependent bending radius, which allows for assessment of the relative bending dissipation in the Earth’s subduction zones using parameters from a recent global compilation. This yields estimates of FB/FTOT < 20%. These results suggest that plates on Earth weakly resist bending, yet are able to propagate a large amount of slab pull.

Modelling subduction involves solving the dynamic interaction between a rigid (solid yet deformable) plate and the fluid (easily deformable) mantle. Previous approaches neglected the solid-like behavior of the lithosphere by only considering a purely fluid description. However, over the past 5 years, a more self-consistent description of a mechanically differentiated subducting plate has emerged. The key fea- ture in this mechanical description is incorporation of a strong core which provides small resistance to plate bending at subduction zones while simultaneously providing adequate stretching resistance such that slab pull drives forward plate motion. Additionally, the accompanying numerical approaches for sim- ulating large-scale lithospheric deformation processes coupled to the underlying viscous mantle flow, have been become available. Here we put forward three fundamentally different numerical strategies, each of which is capabable of treating the advection of mechanically distinct materials that describe the subducting plate. We demonstrate their robustness by calculating the numerically challenging problem of subduction of a 6000 km wide slab at high-resolution in three-dimensions, the successfuly achieve- ment of which only a few codes in the world can presently even attempt. In spite of the differences of the approaches, all three codes pass the simple qualitative test of developing an “S-bend” trench curva- ture previously observed in similar models. While reproducing this emergent feature validates that the lithosphere–mantle interaction has been correctly modelled, this is not a numerical benchmark in the traditional sense where the objective is for all codes to achieve exact agreement on a unique numerical solution. However, we do provide some quantitative comparisons such as trench and plate kinematics in addition to discussing the strength and weaknesses of the individual approaches. Consequently, we believe these developed algorithms can now be applied to study the parameters involved in the dynamics of subduction and offer a toolbox to be used by the entire geoscience community.

Understanding and explaining emergent constitutive laws in the multi-scale evolution from point defects, dislocations and two-dimensional defects to plate tectonic scales is an arduous challenge in condensed matter physics. The Earth appears to be the only planet known to have developed stable plate tectonics as a means to get rid of its heat. The emergence of plate tectonics out of mantle convection appears to rely intrinsically on the capacity to form extremely weak faults in the top 100km of the planet. These faults have a memory of at least several hundred millions of years, yet they appear to rely on the effects of water on line defects. This important phenomenon was first discovered in laboratory and dubbed ‘‘hydrolytic weakening’’. At the large scale it explains cycles of co-located resurgence of plate generation and consumption (the Wilson cycle), but the exact physics underlying the process itself and the enormous spanning of scales still remains unclear.
We present an attempt to use the multi-scale non-equilibrium thermodynamic energy evolution inside the deforming lithosphere to move phenomenological laws to laws derived from basic scaling quantities, develop self-consistent weakening laws at lithospheric scale and give a fully coupled deformation- weakening constitutive framework. At meso- to plate scale we encounter in a stepwise manner three basic domains governed by the diffusion/reaction time scales of grain growth, thermal diffusion and finally water mobility through point defects in the crystalline lattice. The latter process governs the planetary scale and controls the stability of its heat transfer mode.

It is much debated whether the forces associated with the downgoing plate, the overriding plate, passive or active mantle flow are dominant in controlling the paths of plates into the mantle. We investigate the dynamics and energetics for a free subduction system, driven solely by downgoing plate buoyancy, using a finite-element model of a viscoelastic plate with a free surface, sinking into a passive unbounded mantle represented by drag forces. Parameters are varied to study effects of an asthenosphere, ridge push, and a passive overriding plate, for a range of subducting plate viscosities and densities. Such a single, free plate achieves subduction mainly through trench retreat. Most of the energy dissipation occurs in driving the passive mantle response. As a result, the slab's sinking velocity is its Stokes velocity, determined by lithospheric buoyancy and mantle viscosity. The total subduction velocity and dip adjust to minimize bending dissipation in the lithosphere, and are affected by slab rheology as well as buoyancy. A low viscosity asthenosphere and ridge push facilitate plate advance, increasing plate dips and lowering subduction velocity, while suction and buoyancy of a work-free passive overriding plate decreases plate dips, thus increasing subduction and rollback velocities. However, the geometrical relation between the different parameters is the same in all model cases, because the slabs sink according to their Stokes velocity. The free subduction models thus provide a reference that can be used to distinguish the signature of downgoing plate buoyancy from that of other driving forces in global compilations of subduction parameters.

We use two-dimensional numerical experiments to investigate the long-term dynamics of an oceanic slab. Two problems are addressed: one concerning the influence of rheology on slab dynamics, notably the role of elasticity, and the second dealing with the feedback of slab-mantle interaction to be resolved in part 2. The strategy of our approach is to formulate the simplest setup that allows us to separate the effects of slab rheology (part 1) from the effects of mantle flux (part 2). Therefore, in this paper, we apply forces to the slab using simple analytical functions related to buoyancy and viscous forces in order to isolate the role of rheology on slab dynamics. We analyze parameters for simplified elastic, viscous, and nonlinear viscoelastoplastic single-layer models of slabs and compare them with a stratified thermomechanical viscoelastoplastic slab embedded in a thermal solution. The near-surface behavior of slabs is summarized by assessing the amplitude and wavelength of forebulge uplift for each rheology. In the complete thermomechanical solutions, vastly contrasting styles of slab dynamics and force balance are observed at top and bottom bends. However, we find that slab subduction can be modeled using simplified rheologies characterized by a narrow range of selected benchmark parameters. The best fit linear viscosity ranges between 5 􏰢 1022 Pa s and 5 􏰢 1023 Pa s. The closeness of the numerical solution to nature can be characterized by a Deborah number >0.5, indicating that elasticity is an important ingredient in subduction.

International Conference on Geophysics of Slab Dynamics;Jeju Island, South Korea, 19–22 August 2012 In a recent meeting on Jeju Island, South Korea, a community of researchers, half from East Asia and half from North America and Europe, came together to study subduction from complementary points of view—seismology, geodynamics, mineral physics, and tectonics—to enhance the synergy between subduction-related disciplines. The discussion focused on the predictability of the largest earthquakes, the origin of the interplate volcanism in Asia, the challenges for numerical modelers, and how the subduction process influences the evolution of the mantle.

We use a set of 2-D and 3-D numerical fluid dynamic experiments, modeled with different strain rate dependent rheologies (viscous, visco-plastic, power law) to ana- lyze the long-term dynamics of the subduction of an oceanic slab into an iso-viscous or stratified mantle. For the lithosphere a fluid dynamic approach has been bench- marked with our previous solid mechanical approach with the aim of overcoming the coherency problem of fluid dynamic calculations. The solid mechanical dichotomy Sstrong before failure and weak where it failsT has been cast into a specialized non- & cedil;linear fluid rheology. Analog 2-D and 3-D experiments are finally compared with the numerical experiments. 2-D numerical experiments are considered with and without free surface to investigate the limitations induced by a closed top boundary. The effect of asymmetric boundary conditions (with and without overriding plate) is analyzed with respect to the possibility of trench retreat. We clearly state the importance for the free surface analysis. 2-D experiments have inherent weaknesses: first they provide an unrealistic simulation of mantle flow (suppression of toroidal flow), second they give rise to the Sclosed boxT problem (interaction of the slab with a boundary, i.e. & cedil;660 km and the left and right box boundaries). 3-D numerical experiments permit to overcome these problems. A natural analysis of the behavior of the mantle flow during subduction and the three-dimensional behavior of the slab is thus possible. Physical observables like trench retreat and toroidal and poloidal flow are compared with the results of our companion analog 3-D experiments.

The aim of this work is to study the physics of subduction and to develop a dynamical model for back-arc extension in the Central Mediterranean. This is achieved combining 2-D numerical and 3-D laboratory models (Fig. 1) as well as geological and tomographic data of the Central Mediterranean. Concerning the physics of subduction two key problems have been addressed: the influence of rheology on slab dynamics and the feedback of slab-mantle deformation.

Deep geothermal from the hot crystalline basement has remained an unsolved challenge for the geothermal industry for the past 30 years. The original idea of a Hot Dry Rock (HDR) reservoir was to create fluid flow paths between an injection and extraction well deep into the crystalline hot rock by hydraulic fracturing. It became clear in early projects that rather than creating new hydraulic fractures, the existing natural fractures provided the flow paths, and their transmissibility was improved by stimulation.

This article has been retracted at the request of the Editor-in-Chief and Authors. Please see Elsevier Policy on Article Withdrawal (http://www.elsevier.com/locate/withdrawalpolicy).Reason: after publication, errors were discovered in the plate-motion database that it was based on. This dataset was an updated version of the dataset presented in Sdrolias and Muller (2006), provided to us by the first author. The errors in this version were in the away-from-trench/towards-trench assignment for subduction zones with back-arcs and also due to the fact that the next generation plate model had only partially been completed. These errors affect the conclusions about seismic coupling. They also change some of the points in most of the other plots, and although this does not invalidate the other conclusions, the discussion to reach them would be altered.