Semi-empirical Gibbs free energy formulations for minerals and fluids for use in thermodynamic databases of petrological interest

The P-T partition function in statistical thermodynamics can be used to derive semi-empirical formulations of the Gibbs free energy G for minerals and fluids. Parameterization of these equations includes simultaneous regression of experimental heat capacity and molar volume data, allowing fitting, appraisal and optimization of various data sources, as required in the construction of internally consistent petrological data bases. This approach can also be extended to minerals with λ-transitions and to fluids by considering the Gibbs free energy as a function of pressure P, temperature T and an ordering parameter Xα, so that accurate modelled representation and extrapolation of the thermodynamic properties of large numbers of petrologically significant minerals and coexisting fluids can be attained. The ordering parameter is chosen to denote the equilibrium mole fraction (thermodynamic probability) of ordered clusters (structural units) in a substance when G(T,P, Xα)=min. The procedure is tested on existing experimental data for the system MgO-SiO2-H2O. The proposed Gibbs free energy formulation permits thermodynamic properties of minerals, fluids and phase equilibria to be described and extrapolated over a wide range of pressure (0-800 kbar) and temperature (20-3000 K), thus allowing effective use in thermodynamic data bases of petrological interes

Numerical modelling of post-seismic rupture propagation after the Sumatra 26.12.2004 earthquake constrained by GRACE gravity data

In the last decades, the development of the surface and satellite geodetic and geophysical observations brought a new insights into the seismic cycle, documenting new features of inter-, co-, and post-seismic processes. In particular since 2002 satellite mission GRACE provides monthly models of the global gravity field with unprecedented accuracy showing temporal variations of the Earth's gravity field, including those caused by mass redistribution associated with earthquake processes. When combined with GPS measurements, these new data have allowed to assess the relative importance of afterslip and viscoelastic relaxation after the Sumatra 26.12.2004 earthquake. Indeed the observed post-seismic crustal displacements were fitted well by a viscoelastic relaxation model assuming Burgers body rheology for the asthenosphere (60-220 km deep) with a transient viscosity as low as 4× 1017 Pas and constant∼1019 Pas steady state viscosity in the 60-660-km depth range. However, even the low-viscosity asthenosphere provides the amplitude of strain which gravity effect does not exceed 50 per cent of the GRACE gravity variations, thus additional localized slip of about 1 m was suggested at downdip extension of the coseismic rupture. Post-seismic slip at coseismic rupture or its downdip extension has been suggested by several authors but the mechanism of the post-seismic fault propagation has never been investigated numerically. Depth and size of localized slip area as well as rate and time decay during the post-seismic stage were either assigned a priory or estimated by fitting real geodesy or gravity data. In this paper we investigate post-seismic rupture propagation by modelling two consequent stages. First, we run a long-term, geodynamic simulation to self-consistently produce the initial stress and temperature distribution. At the second stage, we simulate a seismic cycle using results of the first step as initial conditions. The second short-term simulation involves three substeps, including additional stress accumulation after part of the subduction channel was locked; spontaneous coseismic slip; formation and development of damage zones producing afterslip. During the last substep post-seismic stress leads to gradual∼1 m slip localized at three faults around∼100-km downdip extension of the coseismic rupture. We used the displacement field caused by the slip to calculate pressure and density variations and to simulate gravity field variations. Wavelength of calculated gravity anomaly fits well to that of the real data and its amplitude provides about 60 per cent of the observed GRACE anomaly. Importantly, the surface displacements caused by the estimated afterslip are much smaller than those registered by GPS networks. As a result cumulative effect of Burgers rheology viscoelastic relaxation (which explains measured GPS displacements and about a half of gravity variations) plus post-seismic slip predicted by damage rheology model (which causes much smaller surface displacements but provides another half of the GRACE gravity variations) fits well to both sets of the real data. Hence, the presented numerical modelling based on damage rheology supports the process of post-seismic downdip rupture propagation previously hypothesized from the GRACE gravity dat

Practical analytical solutions for benchmarking of 2-D and 3-D geodynamic Stokes problems with variable viscosity

Geodynamic modeling is often related with challenging computations involving solution of the Stokes and continuity equations under the condition of highly variable viscosity. Based on a new analytical approach we have developed particular analytical solutions for 2-D and 3-D incompressible Stokes flows with both linearly and exponentially variable viscosity. We demonstrate how these particular solutions can be converted into 2-D and 3-D test problems suitable for benchmarking numerical codes aimed at modeling various mantle convection and lithospheric dynamics problems. The Main advantage of this new generalized approach is that a large variety of benchmark solutions can be generated, including relatively complex cases with open model boundaries, non-vertical gravity and variable gradients of the viscosity and density fields, which are not parallel to the Cartesian axes. Examples of respective 2-D and 3-D MatLab codes are provided with this paper

Numerical modelling of post-seismic rupture propagation after the Sumatra 26.12.2004 earthquake constrained by GRACE gravity data

International audienceIn the last decades, the development of the surface and satellite geodetic and geophysical observations brought a new insights into the seismic cycle, documenting new features of inter-, co-, and post-seismic processes. In particular since 2002 satellite mission GRACE provides monthly models of the global gravity field with unprecedented accuracy showing temporal variations of the Earth's gravity field, including those caused by mass redistribution associated with earthquake processes. When combined with GPS measurements, these new data have allowed to assess the relative importance of afterslip and viscoelastic relaxation after the Sumatra 26.12.2004 earthquake. Indeed the observed post-seismic crustal displacements were fitted well by a viscoelastic relaxation model assuming Burgers body rheology for the asthenosphere (60–220 km deep) with a transient viscosity as low as 4 × 10^17 Pas and constant ~ 10^19 Pas steady state viscosity in the 60–660-km depth range. However, even the low-viscosity asthenosphere provides the amplitude of strain which gravity effect does not exceed 50 per cent of the GRACE gravity variations, thus additional localized slip of about 1 m was suggested at downdip extension of the coseismic rupture. Post-seismic slip at coseismic rupture or its downdip extension has been suggested by several authors but the mechanism of the post-seismic fault propagation has never been investigated numerically. Depth and size of localized slip area as well as rate and time decay during the post-seismic stage were either assigned a priory or estimated by fitting real geodesy or gravity data. In this paper we investigate post-seismic rupture propagation by modelling two consequent stages. First, we run a long-term, geodynamic simulation to self-consistently produce the initial stress and temperature distribution. At the second stage, we simulate a seismic cycle using results of the first step as initial conditions. The second short-term simulation involves three substeps, including additional stress accumulation after part of the subduction channel was locked; spontaneous coseismic slip; formation and development of damage zones producing afterslip. During the last substep post-seismic stress leads to gradual ~1 m slip localized at three faults around ~100-km downdip extension of the coseismic rupture. We used the displacement field caused by the slip to calculate pressure and density variations and to simulate gravity field variations. Wavelength of calculated gravity anomaly fits well to that of the real data and its amplitude provides about 60 per cent of the observed GRACE anomaly. Importantly, the surface displacements caused by the estimated afterslip are much smaller than those registered by GPS networks. As a result cumulative effect of Burgers rheology viscoelastic relaxation (which explains measured GPS displacements and about a half of gravity variations) plus post-seismic slip predicted by damage rheology model (which causes much smaller surface displacements but provides another half of the GRACE gravity variations) fits well to both sets of the real data. Hence, the presented numerical modelling based on damage rheology supports the process of post-seismic downdip rupture propagation previously hypothesized from the GRACE gravity data

ON THE STOKES FLOW COMPUTATION ALGORITHM BASED ON WOODBURY FORMULA

The Stokes approximation is used for the description of flow in nanostructures. An algorithm for Stokes flow computation in cases when there is great variation in the viscosity over a small spatial region is described. This method allows us to overcome computational difficulties of the finite-difference method. The background of the approach is using the Woodbury formula -a discrete analog of the Krein resolvent formula. The particular example of a rectangular domain is considered in detail. The inversion of the discrete Stokes operator is made in analytic form for the case of constant viscosity

A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the ‘sticky air' method

Calculating surface topography in geodynamic models is a common numerical problem. Besides other approaches, the so-called ‘sticky air' approach has gained interest as a free-surface proxy at the top boundary. The often used free slip condition is thereby vertically extended by introducing a low density, low viscosity fluid layer. This allows the air/crust interface to behave in a similar manner to a true free surface. We present here a theoretical analysis that provides the physical conditions under which the sticky air approach is a valid approximation of a true free surface. Two cases are evaluated that characterize the evolution of topography on different timescales: (1) isostatic relaxation of a cosine perturbation and (2) topography changes above a rising plume. We quantitatively compare topographies calculated by six different numerical codes (using finite difference and finite element techniques) using three different topography calculation methods: (i) direct calculation of topography from normal stress, (ii) body-fitting methods allowing for meshing the topography and (iii) Lagrangian tracking of the topography on an Eulerian grid. It is found that the sticky air approach works well as long as the term (ηst/ηch)/(hst/L)3 is sufficiently small, where ηst and hst are the viscosity and thickness of the sticky air layer, and ηch and L are the characteristic viscosity and length scale of the model, respectively. Spurious lateral fluctuations of topography, as observed in some marker-based sticky air approaches, may effectively be damped by an anisotropic distribution of markers with a higher number of markers per element in the vertical than in the horizontal directio

Translithospheric Mantle Diapirism: Geological Evidence and Numerical Modelling of the Kondyor Zoned Ultramafic Complex (Russian Far-East)

We report new structural, microstructural, petrological, and major- and trace-element data on ultramafic rocks from the Kondyor zoned ultramafic complex in Far-East Russia. The ultramafic rocks are subdivided into three subconcentric lithologies, from core to rim: (1) a metasomatic domain where generally phlogopite-rich dykes pervasively intrude dunite; (2) a main dunite core; (3) a pyroxenite rim. The ultramafic rocks have nearly vertical contacts with the surrounding Archaean basement (gneisses, quartzites and marbles) and hornfelsed Riphean sediments. The hornfelsed sediments show a relatively steep (> 60°), outward dipping layering, which rapidly flattens to horizontal away from the inner contact. Although the Riphean sediments define a dome-like structure, the inward, shallow dipping foliation of the dunites indicates a synformal structure. Detailed petro-structural investigations indicate that the Kondyor dunites were deformed by solid-state flow under asthenospheric mantle conditions. The outward textural change from coarse- to fine-grained equigranular dunite and the outward-increasing abundance of subgrains and recrystallized olivine grains suggest dynamic recrystallization while fluid circulation was channelized within the core metasomatic zone, with a decreasing melt fraction from core to rim, and also suggest that solid-state deformation induced grain-size reduction towards the cooling border of the Kondyor massif. Based on their geochemistry, the dunites are interpreted as mantle rocks strongly affected by reaction with melts similar to the Jurassic-Cretaceous Aldan Shield lamproites. Rim pyroxenites were formed by a melt-consuming peritectic reaction, implying the existence of at least a small, conductive thermal gradient around the dunite body while the latter was still at near-solidus temperature conditions. This suggests that the zoned structure of Kondyor was initiated at mantle depths, most probably within the subcontinental lithosphere. Upon cooling, the lamproitic melts were progressively focused in the central part of the massif and drained into vein conduits where they reacted with the wall-rock dunite. Two-dimensional numerical modelling based on finite-differences with a marker-in-cell technique incorporates temperature-dependent rheologies for both molten and non-molten host rocks. The modelling consolidates the structural, petrological and geochemical interpretations, which show that the dunites represent the synformal, flat-lying apex of an asthenospheric mantle diapir, triggered by fluid pressure channelized in the core, which nearly reached the Earth's surface. We conclude that translithospheric mantle diapirism is an important mode of mass transfer in theEart

Modelling gravitational instabilities: slab break-off and Rayleigh-Taylor diapirism

A non-standard new code to solve multiphase viscous thermo–mechanical problems applied to geophysics is presented. Two numerical methodologies employed in the code are described: A level set technique to track the position of the materials and an enrichment of the solution to allow the strain rate to be discontinuous across the interface. These techniques have low computational cost and can be used in standard desktop PCs. Examples of phase tracking with level set are presented in two and three dimensions to study slab detachment in subduction processes and Rayleigh–Taylor instabilities, respectively. The modelling of slab detachment processes includes realistic rheology with viscosity depending on temperature, pressure and strain rate; shear and adiabatic heating mechanisms; density including mineral phase changes and varying thermal conductivity. Detachment models show a first prolonged period of thermal diffusion until a fast necking of the subducting slab results in the break–off. The influence of several numerical and physical parameters on the detachment process is analyzed: The shear heating exerts a major influence accelerating the detachment process, reducing the onset time to one half and lubricating the sinking of the detached slab. The adiabatic heating term acts as a thermal stabilizer. If the mantle temperature follows an adiabatic gradient, neglecting this heating term must be included, otherwise all temperature contrasts are overestimated. As expected, the phase change at 410 km depth (olivine–spinel transition) facilitates the detachment process due to the increase in negative buoyancy. Finally, simple plume simulations are used to show how the presented numerical methodologies can be extended to three dimensions.Peer ReviewedPostprint (author’s final draft

Arc-continent collision: The making of an orogen

Peer Reviewe

Order/disorder phase transition in cordierite and its possible relationship to the development of symplectite reaction textures in granulites

Based on a consistent set of empirical interatomic potentials, static structure energy calculations of various Al/Si configurations in the supercell of Mg-cordierite and Monte Carlo simulations the phase transition between the orthorhombic and hexagonal modifications of cordierite (Crd) is predicted at 1623 K. The temperature dependences of the enthalpy, entropy, and free energy of the Al/Si disorder were calculated using the method of thermodynamic integration. The simulations suggest that the commonly observed crystallization of cordierite in the disordered hexagonal form could be related to a tendency of Al to occupy T1 site, which is driven by local charge balance. The increase in the Al fraction in the T1 site over the ratio of 2/3(T1): 1/3(T2), that characterizes the ordered state, precludes formation of the domains of the orthorhombic phase. This intrinsic tendency to the crystallization of the metastable hexagonal phase could have significantly postponed the formation of the association of orthorhombic cordierite and orthopyroxene over the association of quartz and garnet in metapelites subjected to granulite facies metamorphism. The textures of local metasomatic replacement (the formation of Crd + Opx or Spr + Crd symplectites between the grains of garnet and quartz) indicate the thermodynamic instability of the association of Qtz + Grt at the moment of the metasomatic reaction. This instability could have been caused by the difficulty of equilibrium nucleation of orthorhombic cordierite