Numerical modelling of geoelectric structures using the galerkin and finite elements methods

Summary The Galerkin and finite elements methods are used to model numerically H- and E-polarized electromagnetic fields in inhomogeneous media. The error of approximation is studied and an example of numerical results is given.     

基于非结构网格求解三维D'Alembert介质中声波方程的并行加权Runge-Kutta间断有限元方法

间断有限元方法（Discontinuous Galerkin method,简称DGM）在求解地震波动方程时具有低数值频散、网格剖分灵活等优点,因此,为适应数值模拟对模拟精度和复杂地质结构的要求,本文提出一种新的加权Runge-Kutta间断有限元（weighted Runge-Kutta discontinuous Galerkin,简称WRKDG）方法,用于求解三维D'Alembert介质中声波方程.本文不仅详细推导了其数值格式,特别地,根据常微分方程理论给出了满足数值稳定性条件的一般经验公式,并首次对该方法的数值频散和耗散进行了深入分析,且考虑了耗散参数对结果的影响.同时,我们也对该方法进行了精度测试,并分析了3D情形下WRKDG方法的并行加速比,结果表明3D WRKDG方法具有良好的并行性.最后,我们给出了包含均匀模型、非规则几何模型以及非均匀Marmousi模型在内的数值模拟算例.结果表明,该方法不仅计算准确,能与解析解很好地吻合,且能有效模拟包含球体在内的非规则模型及非均匀Marmousi模型中的衰减声波波场.数值模拟实验进一步验证了WRKDG方法在求解三维D'Alembert介质中声波方程时的正确性和有效性,并获得了对这种强衰减介质中波传播特征的规律性新认识.

     

A discontinuous Galerkin method for two-dimensional flow and transport in shallow water

A discontinuous Galerkin (DG) finite element method is described for the two-dimensional, depth-integrated shallow water equations (SWEs). This method is based on formulating the SWEs as a system of conservation laws, or advection–diffusion equations. A weak formulation is obtained by integrating the equations over a single element, and approximating the unknowns by piecewise, possibly discontinuous, polynomials. Because of its local nature, the DG method easily allows for varying the polynomial order of approximation. It is also “locally conservative”, and incorporates upwinded numerical fluxes for modeling problems with high flow gradients. Numerical results are presented for several test cases, including supercritical flow, river inflow and standard tidal flow in complex domains, and a contaminant transport scenario where we have coupled the shallow water flow equations with a transport equation for a chemical species.     

Coupling of boundary and finite elements for soil-structure interaction problems

The investigation of complex soil-structure interaction problems is usually carried out with numerical solution procedures such as the finite element or the boundary element method. It must be noted, however, that the choice of one or the other of these approaches is not just a matter of preferences; depending on the type of the problem under consideration, either boundary or finite elements may be more advantageous. A considerable expansion in the computational power can be obtained, on the other hand, if one resorts to hybrid schemes which retain the main advantages of the two methods and eliminate their respective disadvantages. This paper presents results obtained with a boundary element-finite element coupling procedure, and discusses its applicability to some representative soil-structure interaction problems. The structures considered are elastic systems, such as foundations, tunnels and filled trenches (modelled by finite elements), which are coupled with homogeneous elastic halfspaces (modelled by boundary elements). The examples demonstrate the importance of using a model that includes wave radiation effects. The coupling approach is formulated entirely in the time domain so that an extension of the algorithm to non-linear analyses seems to present no further difficulties.     

Treatment of time derivative and calculation of flow when solving groundwater flow problems by Galerkin finite element methods

Two techniques connected with the use of the finite element Galerkin method for solving the linear parabolic differential equation describing unsteady groundwater flow in an anisotropic non-homogeneous aquifer are introduced. The first is a mode superposition technique for dealing with the time derivative which involves computing the smallest eigenvalues and associated eigenvectors of the matrices arising from the Galerkin method. It is shown how such a technique allows us to interpret the response of the groundwater level to input in terms of parallel linear reservoirs. It is further argued that if properly implemented, the technique will have computational advantages over standard finite difference methods, e.g. in the case when the input function is constant over relatively large time subintervals. The second is a technique based on so-called generalized flow formulae for calculating flow values across external or internal boundaries, posterior to obtaining the groundwater level values. The implementation of the technique in the case of linear triangular elements on an irregular grid is discussed. It is finally argued from simplified cases that, apart from guaranteeing a match with prescribed input, the technique may often be expected to give more accurate flow values than those obtained directly from the groundwater gradients.     

Solving density driven flow problems with efficient spatial discretizations and higher-order time integration methods

Modelling density driven flow problems requires an excessive computational time and/or heavy equipments due to the non-linear coupling between flow and transport equations. In this work, we develop a robust numerical model with efficient advanced approximations for both spatial and temporal discretizations in order to reduce the excessive computational requirement while maintaining accuracy.     

Positivity-preserving high order well-balanced discontinuous Galerkin methods for the shallow water equations

Shallow water equations with a non-flat bottom topography have been widely used to model flows in rivers and coastal areas. An important difficulty arising in these simulations is the appearance of dry areas where no water is present, as standard numerical methods may fail in the presence of these areas. These equations also have still water steady state solutions in which the flux gradients are nonzero but exactly balanced by the source term. In this paper we propose a high order discontinuous Galerkin method which can maintain the still water steady state exactly, and at the same time can preserve the non-negativity of the water height without loss of mass conservation. A simple positivity-preserving limiter, valid under suitable CFL condition, will be introduced in one dimension and then extended to two dimensions with rectangular meshes. Numerical tests are performed to verify the positivity-preserving property, well-balanced property, high order accuracy, and good resolution for smooth and discontinuous solutions.     

Coupling method of finite and infinite elements for strip foundation wave problems

Starting from two-dimensional wave equations and making use of Galerkin weighted residual approximations, discretized formulations for wave problems of a visco-elastic foundation have been derived. With considerations of geometrical and mechanical characteristics of a semi-infinite domain, a frequency-dependent compatible infinite element is also presented. Finally, by coupling the infinite elements with ordinary finite elements the system is used for simulation of propagating waves in a semi-infinite foundation. This model is not only suitable for simulations of complicated variations of geometrical conditions, but also for describing the unbounded behaviour of arbitrary multiple layers. Examples given indicate the model has excellent computational accuracy and feasibility for analysing the effects on foundation response due to the existence of faults or any other soft layers.     

Comparison of local grid refinement methods for MODFLOW

Many ground water modeling efforts use a finite-difference method to solve the ground water flow equation, and many of these models require a relatively fine-grid discretization to accurately represent the selected process in limited areas of interest. Use of a fine grid over the entire domain can be computationally prohibitive; using a variably spaced grid can lead to cells with a large aspect ratio and refinement in areas where detail is not needed. One solution is to use local-grid refinement (LGR) whereby the grid is only refined in the area of interest. This work reviews some LGR methods and identifies advantages and drawbacks in test cases using MODFLOW-2000. The first test case is two dimensional and heterogeneous; the second is three dimensional and includes interaction with a meandering river. Results include simulations using a uniform fine grid, a variably spaced grid, a traditional method of LGR without feedback, and a new shared node method with feedback. Discrepancies from the solution obtained with the uniform fine grid are investigated. For the models tested, the traditional one-way coupled approaches produced discrepancies in head up to 6.8% and discrepancies in cell-to-cell fluxes up to 7.1%, while the new method has head and cell-to-cell flux discrepancies of 0.089% and 0.14%, respectively. Additional results highlight the accuracy, flexibility, and CPU time trade-off of these methods and demonstrate how the new method can be successfully implemented to model surface water-ground water interactions.     

Orthogonal collocation and alternating-direction procedures for unsaturated flow problems

The alternating-direction collocation method has recently been developed for general parabolic equations. In order to test the applicability of the procedure to highly nonlinear problems, an alternating-direction collocation algorithm is developed to simulate two-dimensional flow in unsaturated porous media. The algorithm employs an alternating-direction solution procedure within the framework of a modified Picard iteration scheme. Numerical behaviour of the new procedure is compared to the behaviour of a standard two-dimensional collocation formulation. The new method is also tested on several infiltration problems of practical interest, including a layered and sloping soil. Results demonstrate the method to be accurate and highly mass conservative. The algorithm also produces significant savings in both execution time and storage.     

A new modeling approach for simulating microtopography-dominated,discontinuous overland flow on infiltrating surfaces

Realistic modeling of discontinuous overland flow on irregular topographic surfaces has been proven to be a challenge. This study is aimed to develop a new modeling framework to simulate the discontinuous puddle-to-puddle (P2P) overland flow dynamics for infiltrating surfaces with various microtopographic characteristics. In the P2P model, puddles were integrated in a well-delineated, cascaded drainage system to facilitate explicit simulation of their dynamic behaviors and interactions. Overland flow and infiltration were respectively simulated by using the diffusion wave model and a modified Green–Ampt model for the DEM-derived flow drainage network that consisted of a series of puddle-based units (PBUs). The P2P model was tested by using a series of data from laboratory overland flow experiments for various microtopography, soil, and rainfall conditions. The modeling results indicated that the hierarchical relationships and microtopographic properties of puddles significantly affected their connectivity, filling–spilling dynamics, and the associated threshold flow. Surface microtopography and rainfall characteristics also exhibited strong influences on the spatio-temporal distributions of infiltration rates, runoff fluxes, and unsaturated flow. The model tests demonstrated its applicability in simulating microtopography-dominated overland flow on infiltrating surfaces.     

Direct and inverse problems in local electromagnetic induction

Various methods of solving direct and inverse problems in local electromagnetic induction are presented. In the section dealing with direct problems some improvments are suggested in the finite difference method in the case of two-dimensional modeling. Two ways of dealing with inverse problems are presented, the first continous, the other parametric. Emphasis is laid upon algebraic aspects of dealing with one-dimensional inverse problems.     

Coupling of finite element and optimization methods for the management of groundwater systems

The paper describes an optimization method for the solution of groundwater management problems. The method consists of a combination of the computation of horizontal plane groundwater flow with a free surface (finite element method) and a linear optimization procedure (simplex algorithm). Considering the special structure of data which result form computing the groundwater flow with the finite element method, and modifying the simplex algorithm, the solution of management problems with complex groundwater flow is realized without any difficulties. Compared to a flow computation alone the additional effort of the optimization (computer time and scope for data storage) is only small.     

Nonlinear iteration methods for nonequilibrium multiphase subsurface flow

Fully implicit, fully coupled techniques are developed for simulating multiphase flow with nonequilibrium mass transfer between phases, with application to groundwater contaminant flow and transport. Numerical issues which are addressed include: use of MUSCL or Van Leer flux limiters to reduce numerical dispersion, use of full or approximate Jacobian for flux limiter methods, and variable substitution for increased Newton iteration efficiency. A comparison of the performance of equilibrium and nonequilibrium models is also presented.     

Higher order time integration methods for two-phase flow

Time integration methods that adapt in both the order of approximation and time step have been shown to provide efficient solutions to Richards' equation. In this work, we extend the same method of lines approach to solve a set of two-phase flow formulations and address some mass conservation issues from the previous work. We analyze these formulations and the nonlinear systems that result from applying the integration methods, placing particular emphasis on their index, range of applicability, and mass conservation characteristics. We conduct numerical experiments to study the behavior of the numerical models for three test problems. We demonstrate that higher order integration in time is more efficient than standard low-order methods for a variety of practical grids and integration tolerances, that the adaptive scheme successfully varies the step size in response to changing conditions, and that mass balance can be maintained efficiently using variable-order integration and an appropriately chosen numerical model formulation.     

Weakly-reflective boundary conditions for two-dimensional shallow water flow problems

In this paper weakly-reflective boundary conditions are derived for the two-dimensional shallow water equations, including bottom friction and Coriolis force. The essential aspects of the derivation are given. Zeroth and first order approximations are applied to the test problem of an initially Gaussian-shaped free surface elevation. For the numerical solution a finite element program is used and various aspects of the numerical implementation are discussed. For small scale practical problems a rather simple (one parameter) formulation might be sufficient. The influence of this parameter is discussed on the weakly-reflectiveness of the boundary condition.     

On different numerical inverse Laplace methods for solute transport problems

Numerical inversion is required when Laplace transform cannot be inverted analytically by manipulating tabled formulas of special cases. However, the numerical inverse Laplace transform is generally an ill-posed problem, and there is no universal method which works well for all problems. In this study, we selected seven commonly used numerical inverse Laplace transform methods to evaluate their performance for dealing with solute transport in the subsurface under uniform or radial flow condition. Such seven methods included the Stehfest, the de Hoog, the Honig–Hirdes, the Talbot, the Weeks, the Simon and the Zakian methods. We specifically investigated the optimal free parameters of each method, including the number of terms used in the summation and the numerical tolerance. This study revealed that some commonly recommended values of the free parameters in previous studies did not work very well, especially for the advection-dominated problems. Instead, we recommended new values of the free parameters for some methods after testing their robustness. For the radial dispersion, the de Hoog, the Talbot, and the Simon methods worked very well, regardless of the dispersion-dominated or advection-dominated situations. The Weeks method can be used to solve the dispersion-dominated problems, but not the advection-dominated problems. The Stehfest, the Honig–Hirdes, and the Zakian methods were recommended for the dispersion-dominated problems. The Zakian method was efficient, while the de Hoog method was time-consuming under radial flow condition. Under the uniform flow condition, all the methods could present somewhat similar results when the free parameters were given proper values for dispersion-dominated problems; while only the Simon method, the Weeks method, and the de Hoog method worked well for advection-dominated problems.     

波动方程反演的全局优化方法研究

复杂介质波动方程反演是地球物理研究中的重要问题,通常表述为特定目标函数最优化,难点是多参数、非线性和不适定性.局部和全局优化方法都不能实现快速全局优化.本文概述了地震波勘探反演问题的理论基础和研究进展,阐述了反演中优化问题的解决方法和面临的困难,并提出了一种确定性全局优化的新方法.通过在优化参数空间识别并划分局部优化解及其附近区域,只需有限次参数空间划分过程就能发现所有局部解（集合）;基于复杂目标函数多尺度结构分析,提出多尺度参数空间分区优化方法的研究方向.该方法收敛速度快,优化结果不依赖初始解的选取,是对非线性全局优化问题的一个新探索.     

Definition and interests of reciprocity and reciprocity gap principles for groundwater flow problems

We introduce the reciprocity and reciprocity gap principles for flow problems in hydrogeology and illustrate their interest in addressing identification problems. The reciprocity principle is derived from mechanics and establishes for flow problems a relationship between different sets of forcing terms, including sources, sinks and boundary conditions, and the resulting head fields. The reciprocity gap principle compares different head fields resulting from the same forcing terms applied to different structures. We give general 2D expressions of the reciprocity and reciprocity gap principles for transient flow problems and give two examples of applications for the identification of transmissivity values and interfaces between different transmissivities. Identification capacities of the reciprocity and reciprocity gap principles yielding direct inversion methods could be used as initial guesses for more advanced inverse problem methodologies.     

Modelling converted seismic waveforms in isotropic and anisotropic 1-D gradients: discontinuous versus continuous gradient representations

Over the past decade, there have been numerous receiver function studies directed at imaging the lithosphere-asthenosphere boundary (LAB). Although it is generally accepted that receiver function phases observed in these studies are derived from physical mode conversions at depth within the lithosphere-asthenosphere transition, it is still debatable as to whether these phases are directly indicative of the LAB. This is because interpretation of receiver function LAB signals relies on understanding the elastic characteristics of the Earth??s outer thermal boundary layer. The main issues for receiver function imaging are the sharpness of the elastic material property transition and, more importantly, what specifically are the material gradients. To test the various transition models, a forward modelling approach is required that allows accurate waveform synthetics for a range of discontinuous and continuous gradients in anisotropic, elastic media. We present a derivation of the reflection and transmission response for continuous one-dimensional (1-D) gradients in generally anisotropic elastic media. We evaluate the influence of 1-D isotropic and anisotropic elastic gradients on the seismic waveform by comparing numerical results of models for discontinuous and continuous transitions. The results indicate that discontinuous representations using layers each with uniform parameters and with thicknesses on the order of approximately 1/3 to 1/8 of the dominant seismic wavelength can be used to accurately model P-to-S and S-to-P mode conversions due to continuous transitions of both isotropic and anisotropic elastic properties. From a practical point of view, when comparing synthetic modelling with observation, this constraint can be relaxed further. The presence of signal noise and/or the result of receiver function stacking techniques will likely obscure these subtle waveform e ff ects. Hence this study suggests that accurate synthetic waveforms for LAB transitions can be modelled with discontinuous gradient representations using a reasonable number of discrete transition layers with layer thicknesses no greater than 1/2 to 1/3 the dominant seismic wavelength.