Tectonic episodicity and convective feedback mechanisms

The evolution of the Earth is characterized by irreversible processes: radioactive decay of the major heat-producing elements, thermal convection and chemical segregation. The prevailing heating from within and the temperature dependence of the viscosity are essential for thermal convection. In the present paper, the chemical and thermal evolution of the mantle and the generation of the continent material are represented by a two-dimensional and finite-difference Boussinesq convection model. We have incorporated the above-mentioned principal features in this model, a geochemical paper by Hofmann (1988, Earth Planet. Sci. Lett., 90: 297–314) constituting our starting point for the distribution of the radionuclides. The concentration of the radionuclides and the viscosity are functions of the location and time developing according to our system of differential equations. Although the real Earth is a much more complex system, we have dared to make a comparison with observed geophysical and geological data; we obtain a depleted upper mantle and acceptable values for the heat flow on the surface of the Earth as well as for the distribution of temperature, viscosity and of the velocity of creep in the mantle. The ups and downs of the convective vigour of the model roughly resemble the supercontinental cycles, the world-wide distribution of mineral dates in time, the sea-level variations and the variations of a number of geochemical parameters.     

一种时空域优化的高精度交错网格差分算子与正演模拟

如何有效压制数值频散是有限差分正演模拟研究中的关键问题之一.近年来,许多学者对二阶声波方程的差分算子开展了大量的优化工作,在压制频散方面取得不错的效果.一阶压强-速度方程广泛用于研究地震波在地下变密度模型中传播规律,目前针对一阶方程的优化工作大多只是在空间差分算子上展开.本文在前人研究的基础上,推导出一阶声波方程中压强场与偏振速度场之间的解析关系,据此在传统交错网格基础上给出一种高精度的显式时间递推格式,该递推格式将时间差分与空间差分算子结合在一起,并采用共轭梯度法得到精确时间递推匹配系数,实现时空差分算子的同时优化.在编程实现算法的基础上,通过频散分析与三个典型模型测试表明:本文方法能够较为有效地压制时间频散与空间频散,提高数值计算精度;同时对复杂模型也有很好适用性.     

一种优化的频率域三维声波有限差分模拟方法

为提高频率域有限差分（FD,finite-difference）正演模拟技术的计算精度和效率,基于旋转坐标系统的优化差分格式被广泛应用,但是只应用于正方形网格的情况.基于平均导数法（ADM）的优化差分格式,应用于正方形和长方形网格模拟.这些频率域有限差分算子,各自具有不同的差分格式和对应的优化系数求解表达式.本文基于三维声波方程发展了一种新的优化方法,只要给定FD模板形式,可直接构造频散方程,求取FD模板上各节点的优化系数.此方法的优点在于频率域FD算子的优化系数对应各个节点,可扩展优化其他格式.运用此优化方法,计算得到了不同空间采样间距比情况下27点和7点格式的优化系数.数值实验表明,优化27点格式与ADM 27点格式具有相同的精度,优化7点格式比经典的7点格式具有更小的数值频散.     

利用低阶偏微分方程组的大倾角差分偏移

利用Claerbout方程进行地震资料偏移,只适用于小倾角的情况。为了克服这一限制,R.Stolt和A.Berkhout等人导出了高阶近似的单程波方程,它们是比较复杂的高阶偏微分方程,在数值求解上存在一定的困难。本文讨论了低阶方程组形式的高阶近似,对它们构造了一些合适的差分格式。提出了求解这些差分方程的具体算法,并与15°差分偏移算法相比较,分析了此算法的计算工作量。本文提出的大倾角差分偏移方法十分有效且容易实现。     

A high order Godunov scheme with constrained transport and adaptive mesh refinement for astrophysical and geophysical MHD

We propose to extend the well-known MUSCL-Hancock scheme for Euler equations to the induction equation modeling the magnetic field evolution in kinematic dynamo problems. The scheme is based on an integral form of the underlying conservation law which, in our formulation, results in a “finite-surface” scheme for the induction equation. This naturally leads to the well-known “constrained transport” method, with additional continuity requirement on the magnetic field representation. The second ingredient in the MUSCL scheme is the predictor step that ensures second order accuracy both in space and time. We explore specific constraints that the mathematical properties of the induction equations place on this predictor step, showing that three possible variants can be considered. We show that the most aggressive formulations reach the same level of accuracy than the other ones, at a lower computational cost. More interestingly, these schemes are compatible with the Adaptive Mesh Refinement (AMR) framework. It has been implemented in the AMR code RAMSES. It offers a novel and efficient implementation of a second order scheme for the induction equation. The scheme is then adaptated to solve for the full MHD equations using the same methodology. Through a series of test problems, we illustrate the performances of this new code using two different MHD Riemann solvers (Lax–Friedrich and Roe) and the need of the Adaptive Mesh Refinement capabilities in some cases. Finally, we show its versatility by applying it to the ABC dynamo problem and to the collapse of a magnetized cloud core.     

Well-balanced high-order centred schemes for non-conservative hyperbolic systems. Applications to shallow water equations with fixed and mobile bed

This paper concerns the development of high-order accurate centred schemes for the numerical solution of one-dimensional hyperbolic systems containing non-conservative products and source terms. Combining the PRICE-T method developed in [Toro E, Siviglia A. PRICE: primitive centred schemes for hyperbolic system of equations. Int J Numer Methods Fluids 2003;42:1263–91] with the theoretical insights gained by the recently developed path-conservative schemes [Castro M, Gallardo J, Parés C. High-order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products applications to shallow-water systems. Math Comput 2006;75:1103–34; Parés C. Numerical methods for nonconservative hyperbolic systems: a theoretical framework. SIAM J Numer Anal 2006;44:300–21], we propose the new PRICE-C scheme that automatically reduces to a modified conservative FORCE scheme if the underlying PDE system is a conservation law. The resulting first-order accurate centred method is then extended to high order of accuracy in space and time via the ADER approach together with a WENO reconstruction technique. The well-balanced properties of the PRICE-C method are investigated for the shallow water equations. Finally, we apply the new scheme to the shallow water equations with fix bottom topography and with variable bottom solving an additional sediment transport equation.     

The unsteady plane flow of ice-sheets: A parabolic problem with two moving boundaries

Abstract

Finite difference algorithms have been developed to solve a one-dimensional non-linear parabolic equation with one or two moving boundaries and to analyse the unsteady plane flow of ice-sheets. They are designed to investigate the response of an ice-sheet to changes in climate, and to reconstruct climatic changes implied by past ice-sheet variations inferred from glacial geological data. Two algorithms are presented and compared. The first, a fixed domain method, replaces time as an independent variable with span. The grid interval in real space is kept constant, and thus the number of grid points changes with span. The second, a moving mesh method, retains time as one of the independent variables, but normalises the spatial variable relative to the span, which now enters the diffusion and advection coeficients in the parabolic equation for the surface profile.

Crank-Nicholson schemes for the solution of the equations are constructed, and iterative schemes for the solution of the resulting non-linear equations are considered.

Boundary (margin) motion is governed by the surface slope at the margin. Differentiation of the evolution equations results in an evolution equation for the margin slopes. It is shown that incorporation of this evolution equation, while not formally increasing the accuracy of the finite difference schemes, in practice increases accuracy of the solution.     

高阶方程偏移的分裂算法

对于用高阶偏微分方程进行地震偏移,本文提出一种分裂算法。有限差分法目前只能解波动方程的二阶和三阶的近似方程的偏移问题,更高阶的近似方程还没有一种有效的解法。利用高阶方程进行偏移会提高偏移的精度和效果。本文就是为了解决这个问题所进行的探索。 本文阐述了分裂算法的原理,并给出了应用该算法解决偏移问题的实例。     

Fourier finite-difference migration for steeply dipping reflectors with complex overburden1

The improvement in accuracy and efficiency of wave-equation migration techniques is an ongoing topic of research. The main problem is the correct imaging of steeply dipping reflectors in media with strong lateral velocity variations. We propose an improved migration method which is based on cascading phase-shift and finite-difference operators for downward continuation. Due to these cascaded operators we call this method‘Fourier finite-difference migration’(FFD migration). In our approach we try to generalize and improve the split-step Fourier migration method for strong lateral velocity variations using an additional finite-difference correction term. Like most of the current migration methods in use today, our method is based on the one-way wave equation. It is solved by first applying the square-root operator but using a constant velocity at each depth step which has to be the minimum velocity. In a second step, the approximate difference between the correct square-root operator and this constant-velocity squareroot operator (the error made in the first step) is implemented as an implicit FD migration scheme, part of which is the split-step Fourier correction term. Some practical aspects of the new FFD method are discussed. Its performance is compared with that of split-step and standard FD migration schemes. First applications to synthetic and real data sets are presented. They show that the superiority of FFD migration becomes evident by migrating steeply dipping reflectors with complex overburden having strong lateral velocity variations. If velocity is laterally constant, FFD migration has the accuracy of the phase-shift method. The maximum migration angle is velocity adaptive, in contrast to conventional FD migration schemes. It varies laterally depending on the local level of velocity variation. FFD migration is more efficient than higher-order implicit FD schemes. These schemes use two cascaded downward-continuation steps in order to attain comparable migration performance.     

Finite strain theories and comparisons with seismological data

All the finite strain equations that we are aware of that are worth considering in connection with the interior of the Earth are given, with the assumptions on which they are based and corresponding relationships for incompressibility and its pressure derivatives in terms of density. In several cases, equations which have been presented as new or independent are shown to be particular examples of more general equations that are already familiar. Relationships for deriving finite strain equations from atomic potential functions or vice versa are given and, in particular it is pointed out that the Birch-Murnaghan formulation implies a sum of power law potentials with even powers. All the equations that survive simple plausibility tests are fitted to the lower mantle and outer core data for the PEM earth model. For this purpose the model data are extrapolated to zero temperature, using the Mie-Grüneisen equation to subtract the thermal pressure (at fixed density) and the pressure derivative of this equation to substract the thermal component of incompressibility. Fitting of finite strain equations to such zero temperature data is less ambiguous than fitting raw earth model data and leads immediately to estimates of the low temperature zero pressure parameters of earth materials. On this basis, using the best fitting equations and constraining core temperature to give an extrapolated incompressibilityK ₀=1.6×10¹¹Pa, compatible with a plausible iron alloy, the following numerical data are obtained: Core-mantle boundary temperature 3770 K Zero pressure, zero temperature densities: lower mantle 4190 kg m^–3 outer core (solidified) 7500 kg m^–3 Zero pressure, zero temperature incompressibility of the lower mantle 2.36×10¹¹PaHowever, an inconsistency is apparent betweenP() andK() data, indicating that, even in the PEM model, in which the lower mantle is represented by a single set of parameters, it is not perfectly homogeneous with respect to composition and phase.     

The Homogeneous Finite-Difference Formulation of the P-SV-Wave Equation of Motion

Two different approaches to finite-difference modeling of the elastodynamic equations have been used: the heterogeneous and the homogeneous. In the heterogeneous approach, boundary conditions at interfaces are treated implicitly; in the homogeneous, they are explicitly discretized. We present a homogeneous finite-difference scheme for the 2-D P-SV-wave case. This scheme represents a generalization of earlier such schemes, being able to model media with arbitrary non-uniformities, provided only that all interfaces are aligned with the numerical grid. We perform a detailed comparison of the generalized homogeneous scheme with the analogous heterogeneous scheme, and show the two schemes to be identical for media with a spatially constynt Poisson's ratio. For media where Poisson's ratio is spatially varying, the schemes differ by terms first-order in the spatial step size. However, a comparison of the numerical results produced by the two schemes shows that the resulting differences are negligible for a wide range of values of the Poisson's ratio contrast.     

全波形时间域航空电磁响应三维有限差分数值计算

全波形时域航空电磁测量具有解释精度高、分辨率高等优点,已成为新一代时域航空电磁探测发展的趋势.为此,本文基于无源有限差分法,将发射电流离散为若干梯形脉冲元,将脉冲元在时间上产生的电磁场转换为空间初始场,在差分迭代过程中分时引入差分方程,逐步将有源有限差分方程转换为无源有限差分方程,最终实现了发射电流为任意波形时,全波形时间域航空电磁响应的三维数值计算.计算结果表明:在梯形波发射时,对于典型地电模型的航空时间域电磁响应,计算稳定时间大于10 ms;与积分方程法相比,发射线圈中心点二次场响应平均误差小于1%.任意源三维全波形数值计算的实现为全波形三维反演和仪器设计奠定了基础.     

Viscosity of peridotite liquid up to 13 GPa: Implications for magma ocean viscosities

The viscosity of synthetic peridotite liquid has been investigated at high pressures using in-situ falling sphere viscometry by combining a multi-anvil technique with synchrotron radiation. We used a newly designed capsule containing a small recessed reservoir outside of the hot spot of the heater, in which a viscosity marker sphere is embedded in a forsterite + enstatite mixture having a higher solidus temperature than the peridotite. This experimental setup prevents spheres from falling before a stable temperature above the liquidus is established and thus avoids difficulties in evaluating viscosities from velocities of spheres falling through a partially molten sample.

Experiments have been performed between 2.8 and 13 GPa at temperatures ranging from 2043 to 2523 K. Measured viscosities range from 0.019 (± 0.004) to 0.13 (± 0.02) Pa s. At constant temperature, viscosity increases with increasing pressure up to  8.5 GPa but then decreases between  8.5 and 13 GPa. The change in the pressure dependence of viscosity is likely associated with structural changes of the liquid that occur upon compression. By combining our results with recently published 0.1 MPa peridotite liquid viscosities [D.B. Dingwell, C. Courtial, D. Giordano, A. Nichols, Viscosity of peridotite liquid, Earth Planet. Sci. Lett. 226 (2004) 127–138.], the experimental data can be described by a non-Arrhenian, empirical Vogel-Fulcher-Tamman equation, which has been modified by adding a term to account for the observed pressure dependence of viscosity. This equation reproduces measured viscosities to within 0.08 log₁₀-units on average. We use this model to calculate viscosities of a peridotitic magma ocean along a liquid adiabat to a depth of  400 km and discuss possible effects on viscosity at greater pressures and temperatures than experimentally investigated.     

一种全局优化的隐式交错网格有限差分算法及其在弹性波数值模拟中的应用

在数值模拟中,隐式有限差分具有较高的精度和稳定性.然而,传统隐式有限差分算法大多由于需要求解大型矩阵方程而存在计算效率偏低的局限性.本文针对一阶速度-应力弹性波方程,构建了一种优化隐式交错网格有限差分格式,然后将改进格式由时间-空间域转换为时间-波数域,利用二范数原理建立目标函数,再利用模拟退火法求取优化系数.通过对均匀模型以及复杂介质模型进行一阶速度-应力弹性波方程数值模拟所得单炮记录、波场快照分析表明:这种优化隐式交错网格差分算法与传统的几种显式和隐式交错网格有限差分算法相比不但降低了计算量,而且能有效的压制网格频散,使弹性波数值模拟的精度得到有效的提高.     

Study on the rheology of subducting slabs

We calculate thermal and phase structures of subducting slabs for different subducting velocities by a modified coupling code of the kinetic phase-transformation equations and the heat-diffusion equation with latent-heat release. Whereafter, we estimate their rheology structures based on the thermal and phase structures from the mineral physical point of view. At shallow depth, the upper layer has a high effective viscosity greater than 10³⁴Pa · s; while the lower layer has a relatively low effective viscosity, which is greater than 10²⁶Pa · s nevertheless. The effective viscosities below the kinetic phase boundary of olivine to wadsleyite decrease obviously, and reach a minimum of 10²²Pa · s. Small areas with higher effective viscosities exist above the depth of about 700 km in subducting slabs, which are produced by lower temperatures that are related with endothermic phase transformation of spinel to perovskite and magnesiowustite. The 1% and 99% isograds of spinel proportion delineate tortuous belts with low effective viscosities, which would affect the geodynamic behavior of subducting slabs.     

基于谱元法的频率域三维海洋可控源电磁正演模拟

高精度、快速有效的正演模拟算法是三维电磁正反演的前提.为了提高海洋电磁三维数值模拟的精度和效率,本文提出利用基于Gauss-Lobatto-Chebyshev（GLC）基函数的谱元法进行海洋可控源三维电磁正演模拟.谱元法结合有限元法和谱方法的优点.我们通过应用伽辽金加权残差法离散二次电场矢量亥姆赫兹方程,在单元内选择混合阶GLC多项式的张量积作为高阶矢量插值基函数,在求解大型稀疏线性方程组时利用直接求解器进行快速求解,从而实现了三维海洋可控源电磁快速高精度正演模拟.一维和三维模型正演结果验证了本文算法的有效性和准确性.典型模型的数值结果表明谱元法是一种有效的三维海洋可控源电磁正演数值方法,能在稀疏网格剖分情况下获得精确的海洋电磁正演模拟响应.     

Convection of a fluid with strongly temperature and pressure dependent viscosity

Plate tectonics on the Earth is a surface manifestation of convection within the Earth’s mantle, a subject which is as yet improperly understood, and it has motivated the study of various forms of buoyancy-driven thermal convection. The early success of the high Rayleigh number constant viscosity theory was later tempered by the absence of plate motion when the viscosity is more realistically strongly temperature dependent, and the process of subduction represents a continuing principal conundrum in the application of convection theory to the Earth. A similar problem appears to arise if the equally strong pressure dependence of viscosity is considered, since the classical isothermal core convection theory would then imply a strongly variable viscosity in the convective core, which is inconsistent with results from post-glacial rebound studies. In this paper we address the problem of determining the asymptotic structure of high Rayleigh number convection when the viscosity is strongly temperature and pressure dependent, i.e. thermobaroviscous. By a method akin to lid-stripping, we are able to extend numerical computations to extremely high viscosity contrasts, and we show that the convective cells take the form of narrow, vertically-oriented fingers. We are then able to determine the asymptotic structure of the solution, and it agrees well with the numerical results. Beneath a stagnant lid, there is a vigorous convection in the upper part of the cell, and a more sluggish, higher viscosity flow in the lower part of the cell. We then offer some comments on the possible meaning and interpretation of these results for planetary mantle convection.     

On the interaction between convection and crystallization in cooling magma chambers

We investigate the interaction of thermal convection and crystallization in large aspect-ratio magma chambers. Because nucleation requires a finite amount of undercooling, crystallization is not instantaneous. For typical values of the rates of nucleation and crystal growth, the characteristic time-scale of crystallization is about 10³–10⁴ s. Roof convection is characterized by the quasi-periodic formation and instability of a cold boundary layer. Its characteristic time-scale depends on viscosity and ranges from about 10² s for basaltic magmas to about 10⁷ s for granitic magmas. Hence, depending on magma viscosity, convective instability occurs at different stages of crystallization. A single non-dimensional number is defined to characterize the different modes of interaction between convection and crystallization.Using realistic functions for the rates of nucleation and crystal growth, we integrate numerically the heat equation until the onset of convective instability. We determine both temperature and crystal content in the thermal boundary layer. Crystallization leads to a dramatic increase of viscosity which acts to stabilize part of the boundary layer against instability. We compute the effective temperature contrast driving thermal convection and show that it varies as a function of magma viscosity and hence composition.In magmas with viscosities higher than 10⁵ poise, the temperature contrast driving convection is very small, hence thermal convection is weak. In low-viscosity magmas, convective breakdown occurs before the completion of crystallization, and involves partially crystallized magma. The convective regime is thus characterized by descending crystal-bearing plumes, and bottom crystallization proceeds both by in-situ nucleation and deposition from the plumes. We suggest that this is the origin of intermittent layering, a form of rhythmic layering described in the Skaergaard and other complexes. We show that this regime occurs in basic magmas only at temperatures close to the liquidus and never occurs in viscous magmas. This may explain why intermittent layering is observed only in a few specific cases.     

3D finite volume groundwater and heat transport modeling with non-orthogonal grids,using a coordinate transformation method

Many popular groundwater modeling codes are based on the finite differences or finite volume method for orthogonal grids. In cases of complex subsurface geometries this type of grid either leads to coarse geometric representations or to extremely fine meshes. We use a coordinate transformation method (CTM) to circumvent this shortcoming. In computational fluid dynamics (CFD), this method has been applied successfully to the general Navier–Stokes equation. The method is based on tensor analysis and performs a transformation of a curvilinear into a rectangular unit grid, on which a modified formulation of the differential equations is applied. Therefore, it is not necessary to reformulate the code in total. We applied the CTM to an existing three-dimensional code (SHEMAT), a simulator for heat conduction and advection in porous media. The finite volume discretization scheme for the non-orthogonal, structured, hexahedral grid leads to a 19-point stencil and a correspondingly banded system matrix. The implementation is straightforward and it is possible to use some existing routines without modification. The accuracy of the modified code is demonstrated for single phase flow on a two-dimensional analytical solution for flow and heat transport. Additionally, a simple case of potential flow is shown for a two-dimensional grid which is increasingly deformed. The result reveals that the corresponding error increases only slightly. Finally, a thermal free-convection benchmark is discussed. The result shows, that the solution obtained with the new code is in good agreement with the ones obtained by other codes.     

Using MODFLOW code to approach transient hydraulic head with a sharp‐interface solution

Saltwater intrusion problems have been usually tackled through analytical models because of its simplicity, easy implementation and low computational cost. Most of these models are based on the sharp‐interface approximation and the Ghyben–Herzberg relation, which neglects mixing of fresh water and seawater and implicitly assumes that salt water remains static. This paper provides insight into the validity of a sharp‐interface approximation defined from a steady state solution when applied to transient seawater intrusion problems. The validation tests have been performed on a 3D unconfined synthetic aquifer, which include spatial and temporal distribution of recharge and pumping wells. Using a change of variable, the governing equation of the steady state sharp‐interface problem can be written with the same structure of the steady confined groundwater flow equation as a function of a single potential variable (?). We propose to approach also the transient problem solving a single potential equation (using also the ? variable) with the same structure of the confined groundwater flow equation. It will allow solving the problem by using the classical MODFLOW code. We have used the parameter estimation model PEST to calibrate the parameters of the transient sharp‐interface equation. We show how after the calibration process, the sharp‐interface approach may provide accurate enough results when applied to transient problems and improve the steady state results, thus avoiding the need of implementing a density‐dependent model and reducing the computational cost. This has been proved by comparing results with those obtained using the finite difference numerical code SEAWAT for solving the coupled partial differential equations of flow and density‐dependent transport. The comparison was performed in terms of piezometric heads, seawater penetration, transition zone width and critical pumping rates. Copyright © 2014 John Wiley & Sons, Ltd.