Modelling shear bands in a volcanic conduit: Implications for over-pressures and extrusion-rates

Shear bands in a volcanic conduit are modelled for crystal-rich magma flow using simplified conditions to capture the fundamental behaviour of a natural system. Our simulations begin with magma crystallinity in equilibrium with an applied pressure field and isothermal conditions. The viscosity of the magma is derived using existing empirical equations and is dependent upon temperature, water content and crystallinity. From these initial conduit conditions we utilize the Finite Element Method, using axi-symmetric coordinates, to simulate shear bands via shear localisation. We use the von Mises visco-plasticity model with constant magma shear strength for a first took into the effects of plasticity. The extent of shear bands in the conduit is explored with a numerical model parameterized with values appropriate for Soufriere Hills Volcano, Montserrat, although the model is generic in nature. Our model simulates shallow (up to approximately 700 in) shear bands that occur within the upper conduit and probably govern the lava extrusion style due to shear boundaries. We also model the change in the over-pressure field within the conduit for flow with and without shear bands. The pressure change can be as large as several MPa at shallow depths in the conduit, which generates a maximum change in the pressure gradient of 10's of kPa/m. The formation of shear bands could therefore provide an alternative or additional mechanism for the inflation/deflation of the volcano flanks as measured by tilt-metres. Shear bands are found to have a significant effect upon the magma ascent rate due to shear-induced flow reducing conduit friction and altering the over-pressure in the upper conduit. Since we do not model frictional controlled slip, only plastic flow, our model calculates the minimum change in extrusion rate due to shear bands. However, extrusion rates can almost double due to the ton nation of shear bands, which may help suppress volatile loss. Due to the paucity of data and large parameter space available for the magma shear strength our model results can only allow for a qualitative comparison to a natural system at this stage. (c) 2007 Elsevier B.V. All rights reserved

Studying the influence of a solid shell on lava dome growth and evolution using the level set method

A finite element formulation of the level set method, a technique to trace flow fronts and interfaces without element distortion, is presented to model the evolution of the free surface of a spreading flow for a highly viscous medium on a horizontal surface. As an example for this class of problem we consider the evolution of an axisymmetric lava dome. Equilibrium configurations of lava domes have been modelled analytically as brittle shells enclosing pressurized magma. The existence of the brittle shell may be viewed as a direct consequence of the strong temperature dependence of the viscosity. The temperature dependence leads to the formation of a thin predominantly elastic-plastic boundary layer along the free surface and acts as a constraint for the shape and flow of the lava dome. In our model, we adopt Iverson's assumption that the thin boundary layer behaves like an ideal plastic membrane shell enclosing the ductile interior of the lava dome. The effect of the membrane shell is then formally identical to a surface tension-like boundary condition for the normal stress at the free surface. The interior of the dome is modelled as a Newtonian fluid and the axisymmetry equations of motion are formulated in a Eulerian framework. We show that the level set is an effective tool to trace and model deforming interfaces for the example of the free surface of a lava dome. We demonstrate that Iverson's equilibrium dome shapes are indeed steady states of a transient model. We also show how interface conditions in the form of surface tension involving higher order spatial derivative (curvature) can be considered within a standard finite element framework

Magma Flow Instabilities in a Volcanic Conduit: Implications for Long-Period Seismicity

Silicic volcanic eruptions are typically accompanied by repetitive Long-Period (LP) seismicity that originates from a small region of the upper conduit. These signals have the capability to advance eruption prediction, since they commonly precede a change in the eruption vigour. Shear bands forming along the conduit wall, where the shear stresses are highest, have been linked to providing the seismic trigger. However, existing computational models are unable to generate shear bands at the depths where the LP signals originate using simple magma strength models. Presented here is a model in which the magma strength is determined from a constitutive relationship dependent upon crystallinity and pressure. This results in a depth-dependent magma strength, analogous to planetary lithospheres. Hence, in shallow highly-crystalline regions a macroscopically discontinuous brittle type of deformation will prevail, whilst in deeper crystal-poor regions there will be a macroscopically continuous plastic deformation mechanism. This will result in a depth where the brittle-ductile transition occurs, and here shear bands disconnected from the free-surface may develop. We utilize the Finite Element Method and use axi-symmetric coordinates to model magma flow as a viscoplastic material, simulating quasi-static shear bands along the walls of a volcanic conduit. Model results constrained to the Soufrière Hills Volcano, Montserrat, show the generation of two types of shear bands: upper-conduit shear bands that form between the free-surface to a few 100 metres below it and discrete shear bands that form at the depths where LP seismicity is measured to occur corresponding to the brittle-ductile transition and the plastic shear region. It is beyond the limitation of the model to simulate a seismic event, although the modelled viscosity within the discrete shear bands suggests a failure and healing cycle time that supports the observed LP seismicity repeat times. However, due to the paucity of data and large parameter space available these results can only be considered to be qualitative rather than quantitative at this stage

Free surfaces modelling based on level sets

We use a finite element formulation of the level set method to model the evolution of the free surface of axi-symmetric spreading flows of highly viscous media on a horizontal plane. We consider specifically the growth of a lava dome as an example however similar problems also occur in flows involving the spreading of molten metals or ceramics. Here we restrict ourselves on constant viscosity fluids for simplicity. In real lavas or melts the viscosity is highly temperature dependent. This manifests itself in the formation of thin predominantly elastic-plastic boundary layers along the free (cold) surfaces of the spreading flows. In our model we follow Iverson [23] who assumes that the thin boundary layer behaves like an ideal plastic membrane shell enclosing the free surface. The effect of the membrane shell is then formally identical to a surface tension-like boundary condition for the normal stress at the free surface

Instabilities across the scales: simple models for shear banding, plate subduction and mantle convection in geodynamics

The Earth shows different modes of deformation in response to thermal or gravitational driving forces. The bulk mantle convects like a viscous fluid on the global scale, while the lithosphere is broken into several plates. They show little internal deformation, but change their shapes and relative positions. Oceanic plate material is generated at divergent margins and recycled into the mantle at subduction zones, on a regional scale. The buoyant continental crust resists subduction and develops meter-scale shear bands during deformation. In this article we review Eulerian finite element (FE) schemes and a particle-in-cell (PIC) FE scheme [15]. Focussing initially on models of crustal deformation at a scale of a few tens of km, we choose a Mohr-Coulomb yield criterion based upon the idea that frictional slip occurs on whichever one of many randomly oriented planes happens to be favourably oriented with respect to the stress field. As coupled crust/mantle models become more sophisticated it is important to be able to use whichever failure model is appropriate to a given part of the system. We have therefore developed a way to represent MohrCoulomb failure within a mantle-convection fluid dynamics code. With the modelling of lithosphere deformation we use an orthotropic viscous rheology (a different viscosity for pure shear to that for simple shear) to define a preferred plane for slip to occur given the local stress field. The simple-shear viscosity and the deformation can then be iterated to ensure that the yield criterion is always satisfied. We again assume the Boussinesq approximation - neglecting any effect of dilatancy on the stress field. Subduction is modelled as a Rayleigh-Taylor instability with dense oceanic lithosphere sinking into less dense sublithospheric mantle. We use a linear viscous rheology for the mantle in this case. Parts of the lithosphere are viscous, others brittle. The values of the dynamic viscosity are different for lithosphere and mantle. The brittle behaviour of parts of the lithosphere can be modelled in the continuum limit by using a viscoplastic rheology. Turning to the largest planetary scale, we present an outline of the mechanics of unified models plate-mantle models and then show how computational solutions can be obtained for such models using Escript. The consequent results for different types of convection are presented and the stability of the observed flow patterns with respect to different initial conditions and computational resolutions is discussed

Dynamic subsidence of Eastern Australia during the Cretaceous

During the Early Cretaceous Australia's eastward passage over sinking subducted slabs induced widespread dynamic subsidence and formation of a large epeiric sea in the eastern interior. Despite evidence for convergence between Australia and the paleo-Pacific, the subduction zone location has been poorly constrained. Using coupled plate tectonic–mantle convection models, we test two end-member scenarios, one with subduction directly east of Australia's reconstructed continental margin, and a second with subduction translated ~ 1000 km east, implying the existence of a back-arc basin. Our models incorporate a rheological model for the mantle and lithosphere, plate motions since 140 Ma and evolving plate boundaries. While mantle rheology affects the magnitude of surface vertical motions, timing of uplift and subsidence depends on plate boundary geometries and kinematics. Computations with a proximal subduction zone result in accelerated basin subsidence occurring 20 Myr too early compared with tectonic subsidence calculated from well data. This timing offset is reconciled when subduction is shifted eastward. Comparisons between seismic tomography and model temperature cross-sections, and an absence of subduction zone volcanism in eastern Australia in the Early Cretaceous provide support for the back-arc basin scenario

Dynamics of slab tear faults: Insights from numerical modelling

Tear resistance at the edge of a slab is an important parameter controlling the evolution of subduction zones. However, compared with other subduction parameters such as plate strength, plate viscosity, plate thickness and trench width, the dynamics of tearing are poorly understood. Here we obtain a first-order understanding of the dynamics and morphology of subduction zones to resistance during tear propagation, by developing and using a novel computational modelling technique for subducting slabs, with side boundaries described by visco-plastic weak zones, developing into tear faults. Our 3D model is based upon a visco-plastic slab that sinks into the less dense mantle, generating poloidal and toroidal flows. The asthenospheric mantle field is static and only develops flow due to the subducting slab. We use the finite element code eScript/Finley and the level set method to describe the lithosphere to solve this fluid dynamics problem. Our results show the importance of tear resistance for the speed of trench migration and for shaping the final geometry of subduction systems. We show that slab tearing along a weak layer can result in a relatively straight slab hinge shape, while increasing the strength in the weak layer results in the curvature of the hinge increasing substantially. High tear resistance at the slab edges may hinder rollback to the extent that the slab becomes stretched and recumbently folded at the base of the domain. Tear resistance also controls whether the subducting lithosphere can experience accelerating rollback velocities or a constant rollback velocity. © 2009 Elsevier B.V. All rights reserved

A model comparison study of large-scale mantle-lithosphere dynamics driven by subduction

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 feature 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 simulating 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 kin wide slab at high-resolution in three-dimensions, the successfuly achievement 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 curvature 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. (C) 2008 Elsevier B.V. All rights reserved