DI21A-0337 

A Comparison of Finite Difference Formulations for the Stokes Equations in the Presence of Strongly Varying Viscosity

* Deubelbeiss, Y (yolanda.deubelbeiss@erdw.ethz.ch AF: AF: AF: AF:

Numerical modeling of geodynamic processes typically requires the solution of the Stokes equations for creeping, highly viscous flows. Since material properties such as effective viscosity of rocks can vary many orders of magnitudes over small spatial scales, the Stokes solver needs to be robust even in the case of highly variable viscosity. Currently, a number of different techniques (e.g. finite element, finite difference and spectral methods) are in use by different authors. Benchmark studies indicate that the accuracy of the velocity solution is satisfying for most methods. The accuracy of deviatoric stresses and pressures, however, is typically less than that of velocities. In the case of highly varying viscosity, some methods even result in oscillating pressures. Over recent years there has been an increased demand for accurate pressures. E.g. melt migration through compacting, two- phase flow materials requires solving equations for the fluid and the solid matrix. Shear-localization in partially molten rocks couples moving fluid within a deforming solid. Therefore it is necessary to have accurate knowledge of pressures, which feed back to the solution. It becomes increasingly important to understand the accuracy of numerical methods for Stokes flow in the presence of large variations in material properties. The objective of this study is therefore to evaluate the accuracy of the pressure solution for a number of numerical techniques. Thereby, we make use of a 2-D analytical solution for the stress distribution inside and around a viscous inclusion in a matrix of different viscosity subjected to pure-shear boundary conditions. Furthermore, numerical simulations have been compared with the analytical solution for density-driven flow (Rayleigh-Taylor instabilities). Results are presented for a staggered grid velocity- pressure finite difference method, a stream function finite difference method and a rotated staggered grid velocity- pressure finite difference method. The staggered grid and stream function formulations require viscosities to be defined both at center and at corner points of control volumes, while the rotated staggered grid finite difference method only requires viscosity defined at center points. We demonstrate that the manner in which viscosities are defined at these locations is of extreme importance for the accuracy of the overall solution. The problem is investigated by studying a simple physical quasi 1-D model with a contact of two media representing the contact between an inclusion embedded in a matrix (2-D case). Analytically and numerically, it is demonstrated that viscosity interpolation using harmonic averaging yields the best results. 2-D numerical results for the above mentioned setups show that for different interpolation methods the errors can vary one order of magnitude. Accuracy of velocity solutions are more than half an order of magnitude better than pressure solutions. The Rayleigh-Taylor instability test, on the other hand, has a weaker sensitivity to viscosity interpolation methods. Results are mainly dependent on the manner in which density is interpolated, which is the driving force in this system. Differences between the three numerical schemes for both setups are secondary compared to the effect of the viscosity interpolation. The best averaging method, for the setups studied here, is a geometric- harmonic averaging of viscosity and an arithmetic averaging of density.

DI21A-0338 

Guaranteed Convergence Rates for Iterative Solutions of Elliptic PDEs with Strongly Varying Coefficients

* Elgersma, M R (michael_elgersma@juno.com), Honeywell, 3660 Technology Drive, Minneapolis, MN 55418, United States Yuen, D A (daveyuen@gmail.com), University of Minnesota, Department of Geology and Geophysics and University of Minnesota Supercomputing Institute, Minneapolis, MN 55455, United States

The convergence rate of iterative solvers is determined by the condition number of the matrix associated with the linearized system of equations. The number of iterations required by Krylov solvers or by Jacobi iteration for the coarsest grid in multigrid methods, both depend on the condition number. A matrix representing an elliptic PDE on a discretized region can be considered as a graph with the diagonal matrix elements representing the vertices, and the off-diagonal matrix elements representing the edges of the graph. Spectral graph theory gives upper and lower bounds for the condition number of any graph in terms of the vertex degree and edge weights. Given the largest ball that fits inside the computational domain, the condition number of a matrix is directly related to the ratio of the sum of the vertex degrees inside the ball, and the sum of the edge weights at the boundary of the ball. This knowledge can be used to construct a grid that has the lowest possible condition number. The lowest condition number and lowest number of cells is obtained by gridding or triangulating the region as coarsely as possible while still preserving the topology of boundaries and surfaces of discontinuity. This is achieved by using small cells near boundaries and larger cells further from boundaries. The resulting iteration number bounds are verified using a system of two coupled Laplacians to solve a mantle convection problem with dozens of sinking spheres whose velocity is several orders of magnitude greater than the surrounding mantle.

DI21A-0339 

Joint Inversion of Mantle Viscosity and Thermal Structure: Applications of the Adjoint of Mantle Convection with Observational Constraints

* Liu, L (lijun@gps.caltech.edu), California Institute of Technology, 1200 E California Blvd, Pasadena, CA 91125, United States Gurnis, M (gurnis@gps.caltech.edu), California Institute of Technology, 1200 E California Blvd, Pasadena, CA 91125, United States

The adjoint method widely used in meteorology and oceanography was introduced into mantle convection by Bunge et al (2003) and Ismail et al (2004). We implemented the adjoint method in CitcomS, a finite-element code that solves for thermal convection within a spherical shell. This method constrains the initial condition by minimizing the mismatch of prediction to observation. Since the present day mantle thermal structure is inferred from seismic tomography, we converted seismic velocity to temperature, an uncertain conversion. Moreover, since mantle viscosity is also uncertain, the inference of mantle initial conditions from tomography is not unique. We have developed a method that incorporates dynamic topography as an additional constraint and are able to jointly invert for mantle viscosity and the seismic to thermal scaling. We assume the thermal structure of present day mantle has the same ¡°pattern" as inferred from tomography, but leave the scaling to temperature as an unknown. The other constraint is the evolving dynamic topography recorded at specific points on earth's surface. From the governing equations of mantle convection, we derive the relations between dynamic topography, thermal anomaly and mantle viscosities. These relations allow a two- layer looping algorithm that inverts for viscosity and thermal anomaly: the inner loop takes the tomographic image as a constraint and the outer loop takes dynamic topography and its rate of change. Starting with incorrect values of thermal anomaly and viscosities, we show with synthetic experiments that all variables converge to their correct values after a finite number of iterations. Our method is examined both in a uniformly viscous mantle and a mantle with depth- and temperature-dependent viscosity. The method has been applied to the descent of the Farallon slab beneath North America.

DI21A-0340 

2D/3D Numerical Modelling of Lithosphere-Mantle Interaction.

* Kaus, B J (boris.kaus@erdw.ethz.ch), Geophysical Fluid Dynamics, ETH Zurich, Schafmattstrasse 30, Zurich, 8093, Switzerland * Kaus, B J (boris.kaus@erdw.ethz.ch), Department of Earth Sciences, University of Southern California, 3651 Trousdale Parkway, Los Angeles, CA 90089-740, United States

Whereas 3D numerical modelling of mantle convection is a standard procedure these days, numerical modelling of lithospheric-scale processes have so far mainly been performed in two spatial dimensions. Part of this is probably due to the somewhat more complicated rheologies that are thought to be important for lithospheric- scale processes. Rather than just behaving viscously, surface-near rocks may also deform in an elastic or plastic (brittle) manner, which manifests itself in faults or shear zones. Moreover, rocks on a lithospheric scale are fairly heterogeneous, undergo complex phase transitions, partial melting and have strain-dependent material properties. All of these complexities result in numerical simulations with a rather large number of adjustable ("free") parameters. In order to get an insight in the physics of geodynamic processes, various workers have therefore concentrated on smaller scale problems (e.g. crustal deformation) in which the "bigger" picture was introduced by kinematically prescribed boundary conditions. Whereas such studies undoubtedly increased our understanding of geodynamic processes, it is also important to understand the lithosphere-mantle system in a more self-consistent manner (in which mantle flow drives lithospheric deformation). From a theoretical point-of-view it is clear how this can be done: take a mantle flow code, add a free-surface, elasticity and plasticity and work on ways to track heterogeneous phases (with phase transitions). From a practical point of view, things appear to be slightly more complicated. The self-consistent free-surface requires a very small time step for stability reasons; elasticity requires tracking of stress tensors and plasticity and material heterogeneities may result in large variations in viscosity over small length scales (which deteriorates the efficiency of multigrid solvers). Questions arise whether the low-order (inadmissible) elements typically used in mantle flow code are still sufficient or whether higher order (admissible) elements are necessary. Here we discuss how we addressed some of these issues in a 2D thermo-mechanical finite element code (SloMo) that has been employed to study coupled lithosphere-mantle interaction processes. A semi-implicit free surface algorithm (SIFSA) was developed to overcome some of the time step restrictions related to the presence of a free surface. Tracers are used to track material properties, stress tensors and temperature and the governing equations are solved on a Lagrangian background solver (with remeshing to allow large deformations). Furthermore we discuss a recently developed 3D parallel finite element code that was written in the PETSc framework, has direct, iterative and multigrid solvers and both linear and higher order elements (as well as tracers to track material properties). Finally we give examples of self-consistent modelling of lithosphere-mantle interaction (a case study on Taiwan and preliminary results on deformation of the Western US), as well as address questions such as: does elasticity modify large-scale geodynamic processes? What is the effect if surface-processes on mantle and crustal deformation? Do we need a free-surface or is a free-slip upper boundary condition sufficient?

DI21A-0341 

Investigating Lithosphere Strength With Thin-Shell Tectonic Modeling

* Moder, C (moder@geophysik.uni-muenchen.de), Dept. of Earth and Env. Sci., Theresienstr. 41, Muenchen, 80333, Germany Carena, S (scarena@geophysik.uni-muenchen.de), Dept. of Earth and Env. Sci., Luisenstr. 37, Muenchen, NJ 80333, Germany

The behavior of many major faults on Earth can only be explained if they are assumed to be much weaker than expected from Byerlee's Law alone. However, there is no agreement over what is a realistic range of friction parameters for faults, or its possible dependency on fault network geometry. Both can be studied with numerical forward modeling, but this requires knowledge of the detailed 3-D geometry of the faults. The latter is now available for most of California, thanks to the SCEC Community Fault Model (southern California) and to the USGS program "3-D Geologic Maps and Visualization" (San Francisco Bay and surrounding region). We model the behavior of the California fault network with the finite-element code SHELLS. We use as input a coarse global grid, with local high-resolution representation of actual faults based on the existing 3-D fault maps. By comparing the simulation results with data on fault-slip rates, we can determine how the faults in this network interact, the role of small faults, and we can quantify the typical fault strength in a continental transform plate boundary setting.

DI21A-0342 

Dynamics of Double Subduction: Numerical Predictions and Critical Observations.

* Mishin, Y (yury.mishin@erdw.ethz.ch), Institute of Geophysics, ETH Zurich, Schafmattstr. 30/HPP, Zurich, 8093, Switzerland Gerya, T (taras.gerya@erdw.ethz.ch), Institute of Geophysics, ETH Zurich, Schafmattstr. 30/HPP, Zurich, 8093, Switzerland Burg, J (jean-pierre.burg@erdw.ethz.ch), Institute of Geology, ETH Zurich, Leonhardstr. 19/LEB, Zurich, 8092, Switzerland

Double subduction is a complex geodynamic process in which two plates following each other are synchronously subducted in the same direction. Double subduction episodes are characteristic for both modern and ancient plate tectonics and are, in particular, inferred in the history of the Himalayan collision zone. However, our knowledge about this process is limited by conceptual schemes and double subduction remains unexplained in terms of physical factors controlling its initiation, duration, dynamics as well as its relation to magmatic activity. We present first results on numerical simulation of double subduction process. Our current high-resolution 2D coupled geochemical-petrological-thermomechanical numerical model employs visco-plastic rheology of rocks and allows simultaneous treatment of heat, mass and water transport, metamorphic phase transformation, partial melting and melt extraction. We studied the influence of different physical factors on initiation, duration and dynamics of the process. In particular we explored the effect of varying convergence rate (0.0 - 7.0 cm/yr), age of the slab (10 - 100 Myrs), water propagation velocity (0 - 10 cm/yr) and dislocation creep activation volume (0.6 - 1.0 J/bar). Depending on these physical parameters (primarily on dimensionless ratio between convergence rate and water percolation velocity) numerical experiments show large variations in double subduction dynamics characterized by (i) different amount of shortening/extension in two simultaneously developing subduction zones, (ii) strong spatial and temporal oscillations of magmatic productivity (separated magmatic episodes) within two parallel volcanic arcs, and (iii) different modes of interaction of two subducting slabs (penetrating/non- penetrating) with the 660 km discontinuity. We compare numerical results with two geologically and geophysically investigated examples of double subduction: the past Karakoram and Kohistan Arcs and the active Izu-Bonin- Marianas and Ryukyu Arcs. Numerical predictions show important similarities with geological information and shed new light in interpreting the natural case stories. In particular, they allow deciphering magmatic and structural/kinematic interplays during double subduction tectonics.

DI21A-0343 

Gale: Large Scale Tectonics Modelling With Free Software

* Landry, W (walter@geodynamics.org), CIG/Caltech, 2750 E. Washington Blvd Suite 210, Pasadena, CA 91107, Hodkinson, L (luke@vpac.org), Victorian Partnership for Advanced Computing, PO Box 201, Carlton South, VIC 3053, Australia

In response to requests from the long timescale tectonics community, we have developed Gale, a parallel 2D and 3D finite element code. Gale's focus is on orogenesis, rifting, and subduction, although it is flexible enough to be applied to such diverse problems as coronae formation on Venus and 3D evolution of crustal fault systems. Gale solves the Stokes and heat transport equations with a large selection of viscous and plastic rheologies. Material properties are tracked using particles, allowing Gale to accurately track interfaces and simulate large deformations. In addition, Gale has a true free surface and a simple programming interface that allows you to plug in your own surface process model. Gale supports a wide variety of boundary conditions, including inflow/outflow, fixed, stress, and static and dynamic friction. Gale has been extensively tested and validated and is exhaustively documented with a 100+ page manual. Gale has been run on everything from laptops to 1000+ processor clusters. Source and prebuilt binaries are freely available at the CIG website. We will discuss Gale's capabilities, present benchmark results, and demonstrate solutions to realistic problems. http://geodynamics.org/cig/software/packages/long/gale/

DI21A-0344 

Present-day Three-Dimensional Temperature Distribution From a Mantle Flow Model.

* Wang, X (wangxin@umich.edu), Dept. of Geological Sciences,University of Michigan, 2534 C. C. Little Building, 1100 North University Ave, Ann Arbor, MI 48109-1005, United States Lithgow-Bertelloni, C (crlb@umich.edu), Dept. of Earth Sciences, University College London, Gower St, London, WC1E 6BT, United Kingdom

In an attempt to understand the global temperature distribution in the mantle and its consequences for Earth structure we construct a model for the instantaneous temperature field of the mantle, assuming downwellings slabs to be the most important source of density heterogeneity and temperature variations. We neglect the contributions due to active upwellings. We use a model for the history of subduction derived from tectonic reconstructions and compute the 3-D velocity field for an incompressible Newtonian fluid. We solve the advection- diffusion equation in steady state for a spherical shell using the finite element package ABAQUS. We choose free-slip, 3000 K velocity-temperature boundary conditions at the core-mantle boundary, and at the surface we constrain velocities to be plate velocities and temperatures to be 300 K. The vertical resolution is on the order of ~5 km at the top and bottom boundary layers, and ~200 x 200 km horizontally. We recover the half-space cooling behavior in the lithosphere and obtain reasonable values of the heat flow, indicating that our predicted temperature field behaves as expected. Not surprisingly, the 3-D variations in the entire mantle are dominated by the presence of slabs in regions of long-lived subduction We use our predicted temperature fields to compute the expected phase assemblage for a mantle of constant bulk composition. We will focus our discussion on the expected topography on major seismic discontinuities (410 and 660) and comparisons to three-dimensional seismological observations.

Author(s) (2007), Title, Eos Trans. AGU, 88(52), Fall Meet. Suppl., Abstract #####-##.