DI31A-1770 INVITED
Scalable Adaptive Mantle Convection Simulation on Petascale Supercomputers
Our goal is to conduct global mantle convection simulations that can resolve faulted plate boundaries, down to 1 km scales. Uniform resolution at these scales would result in meshes with a trillion finite elements, which would elude even petaflops supercomputers. Thus parallel adaptive mesh refinement and coarsening (AMR) is essential. In this talk, we present RHEA, a new generation mantle convection code designed to scale to hundreds of thousands of cores. RHEA is built on ALPS, a parallel octree-based adaptive mesh finite element library that provides new distributed data structures and parallel algorithms for dynamic coarsening, refinement, rebalancing, and repartitioning of the mesh. Using TACC's 579 teraflops Ranger supercomputer, we demonstrate excellent weak and strong scalability of parallel AMR on up to 62,464 cores for the advection- diffusion component of mantle convection and up to 16,384 cores for the full mantle convection problem. With RHEA's adaptive capabilities, we have been able to reduce the number of elements by over three orders of magnitude, thus enabling us to simulate large-scale mantle convection with a finest local resolution of 1.5 km.
DI31A-1771
Multigrid-based Simulation Code For Mantle Convection In Spherical Shell Using Yin-Yang Grid
A new simulation code of mantle convection in a three-dimensional spherical shell is presented. Major innovation of the code comes from an combination of two numerical techniques which we had developed for large-scale simulations of solid earth sciences. First ingredient of the code is the Yin-Yang grid, which is a kind of the overset grid on a sphere. In the Yin-Yang grid, a spherical surface is decomposed into two identical component grids, each of which is a part of low-latitude region defined by the spherical polar coordinate system. The most cruicial benefit of the Yin-Yang grid is its quasi-uniform grid spacings over the entire surface. This significantly relaxes a severe restriction of the time increments coming from the CFL condition. Second ingredient is the "ACuTE" smoothing algorithm, particularly designed for the multigrid solution the flow field of mantle convection. In the ACuTE algorithm, the equations for conservation of mass and momentum for highly viscous and incompressible fluid with strongly variable viscosity are iteratively solved for the velocity and pressure fields in a combination of artificial compressibility and local time-stepping methods. In addition, a careful optimization of the multigrid computations helps to enhance the overall parallel efficiency for a massively parallel environment. The new simulation code is based on the finite-volume discretization using the staggered grid in a spherical shell. Benchmark comparisons for the steady convection for low Rayleigh numbers with previous calculations revealed that accurate results are successfully reproduced not only for isoviscous cases but also for those with moderately temperature-dependent viscosity. We also demonstrated that our code can reproduce the change in convective flow patterns into the "Sluggish-Lid" regime with increasing the viscosity variation up to 104.
DI31A-1772
A Community Benchmark for 2D Cartesian Compressible Convection
The results from five Cartesian, compressible mantle convection codes are reported for constant viscosity
and temperature-dependent viscosity calculations in a unit aspect ratio domain from Rayleigh numbers
104 to 106 and dissipation numbers from 0 to 2.0 are compared.
Four codes are finite element codes and one is a finite volume code.
The finite element codes include both Uzawa and direct solution techniques and use either a penalty or
integrated method (Taylor Hood elements).
All five codes agree to better than 2% when comparing surface heat flux and rms velocity for constant
viscosity, steady, Bousinessq convecion (Blankenbach, 1989 1a, 1b, 1c) with approximately 60 by 60 grid
points (elements).
All five codes agree to better than 2% when comparing surface heat flux and rms velocity for constant
viscosity, steady, extended Bousinessq and Truncated Anelastic liquid convecion with approximately 60 by 60
grid points up to Rayleigh numbers of 2 times 105.
Some codes settle on one-cell solutions while other codes settle on two-cell solutions, and this has a
significant effect on the global flow diagnostics.
We do not yet know if the one-cell versus two-cell solutions result from differences in grids, initial conditions
or reflects differences in the time-evolution algorithm.
For temperature-dependent problems, the difference between one-cell solutions and two-cell solutions is
more pronounced.
As has been previously noted, the balance between work and viscous dissipation is better when the pressure
effect of density is included in the buoyancy (Anelastic Liquid versus Truncated Anelastic Liquid) although
the effect on the surface heat flow and rms-velocity is less than 1-2%.
http://conman.geos.vt.edu/~sdk/benchmarks/compr
DI31A-1773
Empirical error estimation for CitcomCU Stokes solver
An estimate of local error of the finite element solution is invaluable in assessing the quality of the mesh and solution. If errors are concentrated in small regions of the mesh, a local grid refinement can achieve high accuracy without resorting to globally refining the mesh. With the recent progress of Adaptive Mesh Refinement in codes of mantle convection, there is a need for an error estimator that is of good quality while cheap to compute. An empirical error estimator is constructed for the Stokes solver in CitcomCU. Trilinear elements for velocity and constant element for pressure are used in the Stokes solver. Solutions of the Stokes problem are computed on two grids of different resolutions in a 3D Cartesian box. The difference of the solutions is treated as the empirical error of the coarse grid solution. The empirical error is spatially correlated with various fields on the coarse grid, including temperature, temperature gradient, viscosity, viscosity gradient, stress, strain rate, and strain rate energy. The correlation are obtained from various temperature distributions and viscosity laws. The result suggests that the second invariant of the strain rate tensor is a good estimator for the local error of the Stokes solution.
DI31A-1774
Assessing the performance of a parallel MATLAB-based 3D convection code
We are currently building 2D and 3D MATLAB-based parallel finite element codes for mantle convection and melting. The codes use the MATLAB implementation of core MPI commands (eg. Send, Receive, Broadcast) for message passing between computational subdomains. We have found that code development and algorithm testing are much faster in MATLAB than in our previous work coding in C or FORTRAN, this code was built from scratch with only 12 man-months of effort. The one extra cost w.r.t. C coding on a Beowulf cluster is the cost of the parallel MATLAB license for a >4core cluster. Here we present some preliminary results on the efficiency of MPI messaging in MATLAB on a small 4 machine, 16core, 32Gb RAM Intel Q6600 processor-based cluster. Our code implements fully parallelized preconditioned conjugate gradients with a multigrid preconditioner. Our parallel viscous flow solver is currently 20% slower for a 1,000,000 DOF problem on a single core in 2D as the direct solve MILAMIN MATLAB viscous flow solver. We have tested both continuous and discontinuous pressure formulations. We test with various configurations of network hardware, CPU speeds, and memory using our own and MATLAB's built in cluster profiler. So far we have only explored relatively small (up to 1.6GB RAM) test problems. We find that with our current code and Intel memory controller bandwidth limitations we can only get ~2.3 times performance out of 4 cores than 1 core per machine. Even for these small problems the code runs faster with message passing between 4 machines with one core each than 1 machine with 4 cores and internal messaging (1.29x slower), or 1 core (2.15x slower). It surprised us that for 2D ~1GB-sized problems with only 3 multigrid levels, the direct- solve on the coarsest mesh consumes comparable time to the iterative solve on the finest mesh — a penalty that is greatly reduced either by using a 4th multigrid level or by using an iterative solve at the coarsest grid level. We plan to measure how these results scale with problems up to 32GB in size in both 2D and 3D. We will also compare the speed of our code per iteration with a community scalar 2D finite code (CONMAN) and parallel 3D code (CITCOM).
DI31A-1775
Thermal-chemical Mantle Convection Models With Adaptive Mesh Refinement
In numerical modeling of mantle convection, resolution is often crucial for resolving small-scale features. New techniques, adaptive mesh refinement (AMR), allow local mesh refinement wherever high resolution is needed, while leaving other regions with relatively low resolution. Both computational efficiency for large- scale simulation and accuracy for small-scale features can thus be achieved with AMR. Based on the octree data structure [Tu et al. 2005], we implement the AMR techniques into the 2-D mantle convection models. For pure thermal convection models, benchmark tests show that our code can achieve high accuracy with relatively small number of elements both for isoviscous cases (i.e. 7492 AMR elements v.s. 65536 uniform elements) and for temperature-dependent viscosity cases (i.e. 14620 AMR elements v.s. 65536 uniform elements). We further implement tracer-method into the models for simulating thermal-chemical convection. By appropriately adding and removing tracers according to the refinement of the meshes, our code successfully reproduces the benchmark results in van Keken et al. [1997] with much fewer elements and tracers compared with uniform-mesh models (i.e. 7552 AMR elements v.s. 16384 uniform elements, and ~83000 tracers v.s. ~410000 tracers). The boundaries of the chemical piles in our AMR code can be easily refined to the scales of a few kilometers for the Earth's mantle and the tracers are concentrated near the chemical boundaries to precisely trace the evolvement of the boundaries. It is thus very suitable for our AMR code to study the thermal-chemical convection problems which need high resolution to resolve the evolvement of chemical boundaries, such as the entrainment problems [Sleep, 1988].
DI31A-1776
3-D Spherical Mantle Convection with Radial Basis Functions
In the past 25 years a wide variety of numerical methods, such as finite-difference, finite-volume , finite- elements, and pseudospectral methods have been employed to study the problem of 3-D mantle convection. All have specialized strengths but also serious weaknesses. The first three methods are generally considered low-order and can involve high algorithmic complexity (as in triangular elements). Spectrally accurate methods do not practically allow for local mesh refinement and often involve cumbersome algebra. Here, we introduce a new grid/mesh-free approach using radial basis functions (RBFs). It has the advantage of being spectrally accurate for arbitrary node layouts in multi-dimensions with extreme algorithmic simplicity, and naturally permits local node refinement. It has been shown for shallow-water equations and vortex flows that RBFs outperform other numerical methods in the sense that they obtain a much higher accuracy for the same spatial resolution while being able to take unusually large time steps. One virtue of the RBF scheme is the ability to use a simple Cartesian geometry while implementing the required boundary conditions for the temperature, velocity and stresses on a spherical surface of both the outer( planetary surface ) and inner shell ( core-mantle boundary ). The velocity and stress components are expressed in terms of the scalar potential approach (Zebib and Schubert, 1982) and the other remaining variable is the perturbed temperature field. We have studied the problem from the onset of convection to a modest nonlinear regime.
DI31A-1777 INVITED
Plate generation and damage theory
The formation of narrow, rapidly deforming plate boundaries separating strong plate interiors are integral components of the generation of plate tectonics from mantle convection. The development of narrow plate boundaries requires the interaction of a nonlinear rheology and convection. One such nonlinear rheology is two-phase damage theory which employs a non-equilibrium relation between interfacial surface energy, pressure, and viscous deformation, thereby forming a theoretical model for void generation. Two-phase damage theory was recently extended to allow for deformational work to increase the fineness (reduce the grain size) of the matrix phase. We present results testing two-phase damage theory in both a two- dimensional convectively driven system and two-dimensional simple shear calculations where we allow for (1) pure void-generating damage, (2) pure fineness-generating damage, and (3) combined void- and fineness- generating damage. Pure void-generating damage is found to be unsuccessful at producing plate-like features. Fineness-generating damage is successful at inducing plate-like behavior in certain circumstances, including increasing viscosity sensitivity to fineness and certain regimes of damage input and healing rate. Cases with combined void- and fineness-generating damage lead to a spectrum of localized and distributed deformation results. The interaction of micro-cracks and grain size reduction in two-phase damage theory suggests a rheological model for shear localization necessary for the generation of plate tectonic boundaries. We also find that the interplay of damage and the thermal state of the lithosphere as influenced by surface temperatures may have a profound influence on the development of plate tectonics.
DI31A-1778
How to make the most out of level sets for geodynamical modeling
A key challenge in geodynamical modeling is how to represent the interface between different fluids, such that (1) no restrictions on the deformability of the interface in the flow field are imposed, (2) mass is conserved in both fluids, and (3) the discontinuous jumps in fluid properties at the interface are maintained. In three dimensions, these challenges become even more intricate, and the appropriate choice of numerical scheme has a profound influence on the tractability, accuracy, robustness and efficiency of the computational simulation. We present a new approach for modeling multi-phase flows based on level sets. The level set approach has a number of attractive computational aspects, including efficiency, robustness, accuracy, and the ability to handle topological change in the evolving interface in both two and three dimensional flows. A key aspect of our algorithm is the utilization of the extension velocity approach, which accurately couples the flow solver to the interface solver, and avoids spurious mass loss and artificial front repositioning that plague reinitialization techniques. We also present a new numerical scheme to solve the Stokes equation in a ghost fluid setup, which obviates the need for smoothing of sharp discontinuities. Instead, interfaces are represented sharply with sub-grid interpolation. This enables us to compute the correct flow field even in the presence of huge and discontinuous jumps in material parameters, such as viscosity contrasts of several orders of magnitude. We validate our code by reproducing the benchmark results for the isothermal Rayleigh-Taylor instability by van Keken et al. 1997. Our analysis compares our approach with other existing techniques, including particle, marker chain and color function methods, detailing the advantages of this approach and demonstrating its ability to track complex structures in evolving flows. We present several examples of geodynamical problems, where the method presented is particularly promising.
DI31A-1779
Modelling localization during simple shear deformation of the lithosphere
We study the formation of shear bands in 3D numerical experiments of welded viscoplastic layers
undergoing simple shear with a small component of extension or compression. We first examine the spacing
of shear bands in viscoplastic-over-viscous layers with a 2.5D geometry and compare the results with the
equivalent spacing in 2D extensional models. The relationship between spacing and integrated strength of
the layers is found to be similar in each case.
We find 3 distinctive modes of deformation. 1) For a weak lower (viscous) layer, deformation in the upper
layer quickly localizes to a single shear band. 2) As the lower layer becomes stronger - comparable with the
upper layer, multiple shear bands in the upper layer co-exist due to the resistance of the lower layer to strong
localization in the upper. 3) When the strength of the lower layer exceeds that of the upper, it is common for
a single shear band to form in the upper layer with a near-horizontal detachment at the interface of the two
layers.
We also examine how the shear bands evolve and interact in full 3D simulations --- in particular how the initial
orientation of incipient shear bands is inherited as they link into larger, through-going structures. Inherited
structure can be amplified into significant topography as strain accumulates. We examine how the
deformation evolves for each of the observed deformation modes.
http://www.underworldproject.org
DI31A-1780
Links Between Displacement Rates and Erosion in Experimental Tectonic Wedges
Erosional redistribution of mass along Earth's surface modifies the near-surface lithostatic stresses, altering displacement rates and the kinematics within orogens. In this study we use analogue experiments of a deforming sand wedge to systematically examine the impact that erosion may have had on the kinematics of the Argentine Precordilleran fold-and-thrust belt at ~32.5°S. Here, the history of deformation has been superbly documented by others, and that work resolves changes in shortening rates over time throughout the range. Specifically, total shortening rates across the fold-and-thrust belt may have changed over time, and out-of-sequence thrusting may have played an important role accommodating deformation at various times in the history of the fold-and-thrust belt. We hypothesize that such changes may be the response of the fold-and-thrust belt to changing erosion of these ranges. To this end, we have constructed an analogue sandbox experiments whose specific layered rheology is akin to that documented in the Precordillera fold-and-thrust belt in central Argentina. Our contractional experimental apparatus (sandbox) includes a servo-controlled feedback system that allows for a variety of boundary conditions to be applied to the moving wall, including constant displacement rate, time-varying displacement rate, constant loading, and time-varying loading. The application of a loading rate allows us to explicitly investigate feedbacks between topographic construction, erosion, strain softening within the dry sand, and temporal changes in total shortening rates that would be difficult to examine using the constant velocity conditions that are usually applied to the analogue models. We also apply Particle Image Velocimetry (PIV) techniques to digital images from the experimental model to derive high-resolution kinematics and calculate strain, uplift and exhumation rates. Preliminary results indicate that changes in the erosional efficiency in the experimental wedge produce instantaneous changes in displacement rates as topography is built or destroyed. Future analogue experiments will erode material from the simulated Precordillera fold-and-thrust belt according to specific geomorphic transport rules to investigate the impact that climate changes and/or the progressive unroofing of erosionally resistant rocks may have had on the Precordillera fold-and-thrust belt.
DI31A-1781
Postseismic stress relaxation in a layered spherical Earth with linear viscoelastic rheologies
We present the new capabilities offered by the application of the Post-Widder Laplace inversion formula to semianalytical postseismic stress modeling in a spherical, viscoelastic self-gravitating Earth. Within the normal-mode solution framework, the Post-Widder algorithm allows to recover the temporal dependence of the physical observables without the explicit computation of Bromwhich path integrals; in this way, many of the limitations of normal-mode based models are completely bypassed. With this approach, the detail of the rheological profile is limited only by the available computing resources, and any linear viscoelastic rheology can be implemented. At the same time, the numerical codes are much simplified and can be easily parallelized in order to take advantage of multiprocessor systems.
DI31A-1782
Modeling and Testing Yielding Rheologies, Friction, and Magmatic Dikes in a Stokes Flow Simulation
For numerical models of the lithosphere, there are a number of
features that can be quite crucial to understanding the underlying
phenomena. We will focus on how to model and, more importantly, test
three key features in the context of the Stokes flow code Gale.
Yielding rheologies are often play a determinative role in the
evolution of lithospheric systems. While there is an established
literature on implementing yielding rheologies, there is some
confusion about expected fault angles. We will clarify this confusion
and provide analytic, high accuracy benchmarks.
When looking at slab subduction, it can be convenient to only simulate
the subducting slab. The supporting material is then treated as
fixed, with the interface described with friction. One way to
implement this is to have a boundary layer of material. This
additional layer can have unfortunate side effects. We will describe
our implementation of friction that sidesteps that issue.
In the context of Stokes flow simulations, magmatic dikes can be
approximated as regions where material is created and flows out. If
this condition is imposed by specifying the divergence of the velocity
inside the dike, it turns out to be equivalent to specifying a
pressure drop across the boundary of the dike. We will describe our
implementation, benchmarks, and sample simulations of this condition.
http://geodynamics.org/cig/software/packages/long/gale/
DI31A-1783
Parameterized Thermal Convection With a Multi-Layered Fluid and Plates
Thermal convection parameterization has provided a useful tool for investigating the thermal evolution of the Earth and other planets. However, many such parameterizations assume homogenous fluid properties and use averages to account for spatial variations in material properties. While reasonable approximations, these models fail to produce the appropriate wavelengths found in systems with depth dependent viscosities, such as the Earth, and are difficult to couple to surface plates. The purpose of this study was to develop a parameterized model that properly accounts for thermal convection through a fluid layer with a depth- dependent viscosity and plates on the surface. We have carried out a study of the dynamics of thermal convection in a fluid layer bounded by plates and with an internal viscosity change. The problem was decomposed into two separate and less complex flows; a buoyantly driven flow and a plate driven flow, with the superposition of the two yielding the total flow. Simple flow models were constructed for both the buoyantly driven and plate driven flows. The two flows were then coupled together using (1) the combined velocity field, (2) simple thermal boundary layer theory to approximate the magnitude of the thermal anomalies, (3) the requirement that the shear stress applied to the plate by the buoyantly driven flow match the shear stress required to force the plate driven flow, and (4) the energy required for subduction in the energy balance through a simple parameterization of plate bending. The combined model provides a parameterization for convective heat transport and yields an approximate velocity field from which the mixing rate across the viscosity change can be estimated. This mixing rate can be used to calculate chemical mixing between the different layers and can provide a means of further constraining thermal histories.
DI31A-1784
Multipathing in Models of Subducting Slabs
The geodynamic controls on the geometry of subducting slabs are still poorly understood. Previous mantle convection studies have shown that for instance different viscosity structures can result in different slab shape, width, and sharpness of slab edges. No direct information is available to study in particular the latter two, and so far seismic tomography has not been able to provide us with models with high enough resolution to determine the exact shape of slabs in the mantle. A tool that could yield more information on slab geometries uses the multipathing principle. This is based on seismic energy arriving in a different manner by presence of sharp edges along the ray path, causing waveforms to be distorted in the seismogram. In this study, several test cases have been produced in which the mantle convection code CitcomT is used to create various slab geometries in 3-D. The models are time-dependent and contain the evolution of the slab from subduction initiation onward. The rheology in these models allows for yielding, Newtonian and non- Newtonian viscosity. The developed thermal structure of the mantle is converted to a velocity structure for which synthetic seismograms are constructed, and the resulting multipathing is analysed.
DI31A-1785
Dynamical Modeling for Generally Shaped, Layered Lithospheric Geometries Using Continuous Field Variables
Modeling lithospheric deformation presents many numerical challenges, e.g., large horizontal to vertical aspect ratio of the area geometry, discontinuous material properties, and layers that pinch out. We modify the standard Euler-Lagrange variational approach for constructing the governing equations to improve the stability and accuracy of numerical calculations. The key innovation is that, for the geophysical flow and elastic deformation problems of interest, the calculations are performed by separating all spatial dependences into terms that can be related to coordinates on a reference surface and variations perpendicular to that surface. This results in considerable computational efficiency, as the problems are reduced to finite element calculations only in terms of the coordinates on the reference surface. The second key ingredient is that the computations are performed using only physical field quantities that are continuous across discontinuities in the medium, i.e., velocities/displacements and tractions. The advantage in this case is computational accuracy for laterally varying layers, including layers that pinch out. To formulate the problem in this way requires an additional modification of the standard Euler-Lagrange variational approach for constructing the governing equations. We introduce the key elements of solving the force balance equations using COMSOL Multiphysics. First, we will outline the geometry including the transformed vertical coordinates, the forms of tangential derivative and surface normals. We then explain the variational approach, including the forms of the functional minimized and the boundary constraints. Last, we discuss how medium properties are incorporated, to show that the methodology can handle rheologies that are linear and non linear, spatially varying, and isotropic or anisotropic, followed by detailed expression for the isotropic case. We present selected examples for 2D and 3D geometries in Cartesian and curvilinear coordinates. The examples include layers that pinch out and discontinuities due to slip on fault surfaces.
DI31A-1786 INVITED
A Boundary Integral Formulation for Multiparticle Flow
This presentation outlines a Boundary Integral Method (BIM) to calculate particle shapes and boundary velocities in an aggregate of viscous particles concentrated in a melt of a different viscosity. Interaction among the particles is controlled by midrange forces and the velocity fields of the neighboring particles. The calculated grain shapes are used to calculate the average elastic moduli of the aggregate, from the intergranular contact area. The results from a series of BEM experiments indicate that the average area of intergranular contact depends on the viscosity contrast between the grains and the melt and the ratio between surface tension and viscosity.
DI31A-1787
Polarity Reversals in Chaotic Geodynamo Models
Polarity reversals in numerical dynamos driven by thermo-chemical convection are analyzed in terms of their magnetic field intensity, transition field symmetry, and other observable characteristics. We compile magnetic field statistics and reversal properties of dynamos with Ekman numbers E=5-10× 10-4, Prandtl numbers Pr=1 and Pm=10, chemical and thermal Rayleigh numbers in the ratio of 1.4, and diffusivity ratio r=1, over 1-2 Myr model time periods. The most earth-like dynamos are characterized by long stable polarity chrons with dipole-dominent time-averaged fields punctuated by occasional polarity reversals, and are found within a transitional range of Rayleigh numbers Ra, between non-reversing dynamos at lower Ra and chaotic multipolar-type dynamos at higher Ra. Dynamos in this transitional regime have time-averaged Lowes magnetic field spectra with elevated dipole terms, reduced quadrupole terms, and nearly flat spectra at larger spherical harmonic degree l-values. They also have broadband dipole intensity spectra, with variance that decreases at low frequencies like f- 2, and more negative slopes at frequencies greater than f=1 Kyr-1. Histograms (PDFs) of the axial dipole intensity are quasi-gaussian. The more earth-like reversing dynamos have trimodal PDFs, with large positive and negative modes representing the stable polarity states, separated by a smaller intermediate mode representing the transition weak field state. Polarity reversals in the transitional regime are preceded by a ~ 5-fold reduction in dipole intensity and a ~ 2-fold reduction in the rms field intensity on the core-mantle boundary, typically lasting ~ 3 times longer than the dipole axis instability. The rate of dipole intensity reduction is variable, sometimes reaching -4% per century sustained decrease. The polarity reversals are accompanied by a transient change in magnetic field symmetry. Whereas the antisymmetric part of the field (the dipole family) exceeds the symmetric part of the field (the quadrupole family) during stable polarity times, these two parts are comparable to each other prior to and also during the reversal.
DI31A-1788 INVITED
Two-Phase Dynamics Simulations of the Growth and Instability of Earth's Inner Core
When the center of Earth's core began to freeze from a homogeneous liquid 1-2 billion years ago, its constitution was very likely that of a mushy region. As this incipient inner core grew by further crystallization of the outer core, an increase in gravity force allowed for the solid grains to compress against one another, undergo viscous compaction, and begin to expel remnant fluid out of the inner core by percolation. Meanwhile, inside the inner core the residual fluid and solid remained in equilibrium, and any perturbations that resulted in upwelling of the deformable mush would also be accompanied by decompression melting. Upwelling and melting regions might then increase in liquid fraction, become less dense, and hence buoyant in a way that would propel them upward at a faster rate, setting up a runaway instability and partial Rayleigh-Taylor-like overturn of Earth's inner core. Structures inherited from this event possibly include the distinct innermost inner core posited by seismologists to exist at Earth's centermost 300-600 km. We use a new two-phase dynamics code to model this scenario in axi-symmetric geometry in order to understand whether and when such an instability occurred, what size the core will have been at the onset of instability, and the degree and style of deformation that would have accompanied this episode. We have found that the growth of instability competes with the rate of background melt percolation, such that the instability would only have occurred after the inner core reaches a critical size and expelled a certain amount of liquid from its interior. A linear stability analysis confirms that there is a critical Rayleigh number for the onset of instability at a given radius. The combined constraints show that the inner core is guaranteed to have undergone this kind of instability, at a time and strength governed solely by physical properties such as grain size, density differences between liquid and solid, and viscosities of the phases.
DI31A-1789
Disequilibrium melting of a two phase multicomponent mantle
There are a number of geochemical and petrological arguments that suggest magma is not always in chemical equilibrium with the surrounding matrix as it rises to the Earth's surface. The most commonly used equations for modelling melt transport are the two phase flow equations of McKenzie (1984), and we have extended these equations to encompass multicomponent mantle melting with possible chemical disequilibrium between the two phases. The governing equations consist of conservation of mass, momentum, and energy; equations of state; and phenomenological laws relating fluxes of mass and heat to differences in temperature and chemical potential. The theory of non-equilibrium thermodynamics constrains the form of these phenomenological laws, through the necessity of positive entropy production (the second law) and the Onsager reciprocal relations. Our particular focus has been on the phenomenological laws for interphase mass transfer, which control melting and crystallisation. In different limits, fractional melting, fractional crystallisation and equilibrium melting can be recovered, but the laws also describe more general behaviour. We have performed some simple 1D melting column calculations to explore these phenomenological laws. The theory should prove useful in a wide range of two phase problems where kinetics plays a key role, and particularly in modelling melt-rock reactions on the continuum scale.
DI31A-1790
Mantle Upwelling and Melting Beneath Oblique and Segmented Ridge Axes
A mid-ocean ridge plate boundary is often idealized as a straight ridge axis with spreading orthogonal to the axis. However, observed ridges often have obliquely spreading segments and/or en-echelon transform faults interrupting much straighter spreading segments. Here we use our new 3D algorithms for viscous mantle flow and melting to compare and contrast the flow and melting patterns arising from these offset patterns. 3D flow is modeled using a MATLAB-based parallel finite element solver for variable viscosity flow. To ensure high numerical resolution beneath displaced ridge axes we use unstructured 3D meshes built up of 10-node tetrahedra. Temperature advection is performed using a semi-Lagrange algorithm on the same unstructured mesh, while diffusion is solved using standard operator- splitting methods. Our first experiments have assumed the mantle melts as a two-component basalt-eclogite + peridotite lithology. Melt is assumed to rise vertically from where it is generated. We also assume the temperature and depletion-dependent rheology of Phipps Morgan (EPSL, 1997), and run a suite of numerical experiments where melt and depletion-buoyancy may or may not occur. Preliminary results suggest that the reduction in upwelling and melting may significantly depend on the style of the ridge offset i.e. oblique ridge vs. staircase of transforms vs. single transform offset. A strong reduction in upwelling is seen beneath even small-offset oblique segments, consistent with the crustal thickness variations seen along the segments of the Southern M.A.R. near Ascension Island.