We use a finite element (FEM) formulation of the level set method to model geological fluid flow problems involving interface propagation. Interface problems are ubiquitous in geophysics. Here we focus on a Rayleigh-Taylor instability, namely mantel plumes evolution, and the growth of lava domes. Both problems require the accurate description of the propagation of an interface between heavy and light materials (plume) or between high viscous lava and low viscous air (lava dome), respectively. The implementation of the models is based on Escript which is a Python module for the solution of partial differential equations (PDEs) using spatial discretization techniques such as FEM. It is designed to describe numerical models in the language of PDEs while using computational components implemented in C and C++ to achieve high performance for time-intensive, numerical calculations. A critical step in the solution geological flow problems is the solution of the velocitypressure problem. We describe how the Escript module can be used for a high-level implementation of an efficient variant of the well-known Uzawa scheme (Arrow et al., 1958). We begin with a brief outline of the Escript modules and then present illustrations of its usage for the numerical solutions of the problems mentioned above.

No abstract

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.

[1] Modeling volcanic phenomena is complicated by free surfaces often supporting large rheological gradients. Analytical solutions and analogue models provide explanations for fundamental characteristics of lava flows. However, more sophisticated models are needed, incorporating improved physics and rheology to capture realistic events. To advance our understanding of the flow dynamics of highly viscous lava in Peléean lava dome formation, axisymmetrical finite element method models of generic endogenous dome growth have been developed. We use a novel technique, the level set method, which tracks a moving interface, leaving the mesh unaltered. The model equations are formulated in an Eulerian framework. In this paper we test the quality of this technique in our numerical scheme by considering existing analytical and experimental models of lava dome growth which assume a constant Newtonian viscosity. We then compare our model against analytical solutions for real lava domes extruded on Soufrière, St. Vincent, West Indies in 1979 and Mount St. Helens in October 1980 using an effective viscosity. The level set method is found to be computationally light and robust enough to model the free surface of a growing lava dome. Also, by modeling the extruded lava with a constant pressure head this naturally results in a drop in extrusion rate with increasing dome height, which can explain lava dome growth observables more appropriately than when using a fixed extrusion rate. From the modeling point of view, the level set method will ultimately provide an opportunity to capture more of the physics while benefiting from the numerical robustness of regular grids.

S U M M A R YA 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.