Solution of the nonlinear transport equation using modified Picard iteration

The transport and fate of reactive chemicals in groundwater is governed by equations which are often difficult to solve due to the nonlinear relationship between the solute concentrations for the liquid and solid phases. The nonlinearity may cause mass balance errors during the numerical simulation in addition to numerical errors for linear transport system. We have generalized the modified Picard iteration algorithm of Celia et al.⁵ for unsaturated flow to solve the nonlinear transport equation. Written in a ‘mixed-form’ formulation, the total solute concentration is expanded in a Taylor series with respect to the solution concentration to linearize the transport equation, which is then solved with a conventional finite element method. Numerical results of this mixed-form algorithm are compared with those obtained with the concentration-based scheme using conventional Picard iteration. In general, the new solver resulted in negligible mass balance errors (< ∥10⁻⁸∥%) and required less computational time than the conventional iteration scheme for the test examples, including transport involving highly nonlinear adsorption under steady-state as well as transient flow conditions. In contrast, mass balance errors resulting from the conventional Picard iteration method were higher than 10% for some highly nonlinear problems. Application of the modified Picard iteration scheme to solve the nonlinear transport equation may greatly reduce the mass balance errors and increase computational efficiency.     

ABAQUS中动力问题边界条件的选取

数值模拟是研究半无限域中结构地震反应的重要手段,而边界条件的选取则直接影响到数值模拟结果的可靠性.本文应用有限元软件ABAQUS,分别采用固定边界条件、黏性边界条件和远置边界条件对简单弹性半无限介质的波源问题和散射问题进行了数值模拟分析,结果显示:使用固定边界时,波动不能穿过边界面进入半无限介质,这与波动在半无限介质中...     

A two-grid method for coupled free flow with porous media flow

This paper presents a two-grid method for solving systems of partial differential equations modelling incompressible free flow coupled with porous media flow. This work considers both the coupled Stokes and Darcy as well as the coupled Navier-Stokes and Darcy problems. The numerical schemes proposed are based on combinations of the continuous finite element method and the discontinuous Galerkin method. Numerical errors and convergence rates for solutions obtained from the two-grid method are presented. CPU times for the two-grid algorithm are shown to be significantly less than those obtained by solving the fully coupled problem.     

一种基于全局优化的交错网格有限差分法

在地震波场数值模拟中, 交错网格有限差分技术得到了广泛的应用, 但是在弹性模量变化较大时, 通常会因插值而导致模拟误差增大. 旋转交错网格可以很好地克服这个缺点, 因而适合于各向异性介质正演模拟. 但是对于同样大小的网格单元, 旋转交错网格需要的步长比常规交错网格要大, 这会使梯度和散度算子的误差增大因而更易产生空间数值频散. 针对这些问题, 本文提出了旋转交错网格与紧致有限差分相结合的方法, 并基于模拟退火算法进行全局优化, 压制数值频散, 拓宽波数范围. 数值模拟结果表明, 此方法可以有效地压制数值频散, 且具有较高的模拟精度.      

Numerical simulation of seismic wavefields in TTI media using the rotated staggered-grid compact finite-difference scheme

Numerical simulation in transverse isotropic media with tilted symmetry axis(TTI) using the standard staggered-grid finite-difference scheme(SSG)results in errors caused by averaging or interpolation. In order to eliminate the errors, a method of rotated staggered-grid finite-difference scheme(RSG) is proposed. However, the RSG brings serious numerical dispersion. The compact staggered-grid finite-difference scheme(CSG) is an implicit difference scheme, which use fewer grid points to suppress dispersion more effectively than the SSG. This paper combines the CSG with the RSG to derive a rotated staggered-grid compact finite-difference scheme(RSGC). The numerical experiments indicate that the RSGC has weaker numerical dispersion and better accuracy than the RSG.     

Numerical analysis of one-dimensional unsaturated flow in layered soils

The paper deals with numerical solutions to the Richards equation to simulate one-dimensional flow processes in the unsaturated zone of layered soil profiles. The equation is expressed in the pressure-based form and a finite-difference algorithm is developed for accurately estimating the values of the hydraulic conductivity between two neighboring nodes positioned in different soil layers, often referred to as the interlayer hydraulic conductivity. The algorithm is based upon flux conservation and continuity of pressure potential at the interface between two consecutive layers, and does not add significantly to simulation run time. The validity of the model is established for a number of test problems by comparing numerical results with the analytical solutions developed by Srivastava and Yeh²⁹ which hold for vertical infiltration towards the water table in a two-layer soil profile. The results show a significant reduction in relative mass balance errors when using the proposed model. Some specific insights into its numerical performance are also gained by comparisons with a numerical model in which the more common geometric averaging operator acts on the interlayer conductivities.     

Impact of Local Groundwater Flow Model Errors on Transport and a Practical Solution for the Issue

A groundwater flow model is typically used to provide the flow field for conducting groundwater solute transport simulations. The advection term of the mass conserved formulation for groundwater transport assumes that the flow field is perfectly balanced and that all water flowing into a numerical grid cell is exactly balanced by outflows after accounting for sources/sinks or internal storage. However, in many complicated regional or site‐scale models, there may be localized flow balance errors that may be difficult to eliminate through tighter flow convergence tolerances due to simulation time constraints or numerical limits on convergence tolerances. Thus, if water is erroneously gained or lost within a grid cell during the flow computation, the solutes within it will also be numerically affected in the associated transport simulation. Transport solutions neglect this error in groundwater flow as the transport equations that are solved assume no error in flow. This flow imbalance error can however have consequences on the transport solution ranging from unnoticeable errors in the resulting concentrations to spurious oscillations that can grow in time and hinder further solution. An approach has been suggested here, to explicitly handle these flow imbalances during mass conserved advective transport computations and report them in the corresponding transport mass balance output, as corrections that are needed to handle errors originating in the flow solution. Example problems are provided to explain the concepts and demonstrate the impacts.     

Numerical versus field studies of delayed yield in response to a moving water-table

As pumping proceeds in an unconfined aquifer the water-table is lowered and slowly drains out of the initially saturated material giving a delayed yield from storage. Although considerable work has been carried out on delayed yield, a definitive solution is still lacking. The unsaturated zone was treated here as a significant factor in delayed yield analysis. The one-dimensional unsaturated flow equation was solved numerically to simulate the flow in the zone located between successive water levels in the unconfined aquifer for the case with zero flux at the upper boundary. The moisture content distribution curve was integrated numerically to determine the amount of water drained with time. The study shows that the numerical results are in close agreement with the monitored field data, which implies the important role of the unsaturated zone and the reliability of the numerical simulation in describing this physical phenomenon. The results show that the relationship between the flux and time is sensibly in exponential form, which agrees with the general applicability of Boulton's assumption.     

Numerical simulation of an unsaturated flow equation

A numerical model for an unsaturated flow problem by using the finite element method is established in order to simulate liquid moisture flow In an unsaturated zone with homogeneous soil and deep subsurface water, and with different initial and boundary conditions. For infiltration or evaporation problems, a traditional method usually yields oscillatory non-physics profiles. However, nonoscillatory solutions are obtained and non-physics solutions for these problems are evaded by using the mass-lumped finite element method. Moreover, the kind of boundary condition is handled very well. Project supported by the National Key Project of Fundamental Research ”Climate Dynamics and Climate Prediction Theory“ and China Postdoctoral Science Foundation.     

Flow in unsaturated random porous media,nonlinear numerical analysis and comparison to analytical stochastic models

This work presents a rigorous numerical validation of analytical stochastic models of steady state unsaturated flow in heterogeneous porous media. It also provides a crucial link between stochastic theory based on simplifying assumptions and empirical field and simulation evidence of variably saturated flow in actual or realistic hypothetical heterogeneous porous media. Statistical properties of unsaturated hydraulic conductivity, soil water tension, and soil water flux in heterogeneous soils are investigated through high resolution Monte Carlo simulations of a wide range of steady state flow problems in a quasi-unbounded domain. In agreement with assumptions in analytical stochastic models of unsaturated flow, hydraulic conductivity and soil water tension are found to be lognormally and normally distributed, respectively. In contrast, simulations indicate that in moderate to strong variable conductivity fields, longitudinal flux is highly skewed. Transverse flux distributions are leptokurtic. the moments of the probability distributions obtained from Monte Carlo simulations are compared to modified first-order analytical models. Under moderate to strong heterogeneous soil flux conditions (σ²_y≥1), analytical solutions overestimate variability in soil water tension by up to 40% as soil heterogeneity increases, and underestimate variability of both flux components by up to a factor 5. Theoretically predicted model (cross-)covariance agree well with the numerical sample (cross-)covarianaces. Statistical moments are shown to be consistent with observed physical characteristics of unsaturated flow in heterogeneous soils.©1998 Elsevier Science Limited. All rights reserved     

A semi-analytical solution based on a numerical solution of the solute transport of a conservative and nonreactive tracer

In the evaluation of potential risk from ingestion of groundwater near an impacted site, numerical simulation of fate and transport processes of chemicals of concern is often required. If there is potential concern about multiple chemicals, numerical simulation of each chemical separately is often needed. In this paper, a semi-analytical solution is presented based on a numerical solution of the transport of a conservative and nonreactive tracer. When multiple chemicals undergoing sorption and first-order degradation need to be modeled, we can avoid performing individual numerical simulations for each chemical by applying the semi-analytical solution. Numerical test runs were conducted to verify the semi-analytical solution; simulation results reveal that the concentrations derived from the semi-analytical solution are identical to those derived from the individual numerical fate and transport model simulations. The semi-analytical solution requires steady-state flow conditions, no continuing contaminant source, and similar initial source concentration distributions.     

Finite-element analysis of the transport of water and solutes in title-drained soils

The finite-element method based on a Galerkin technique was used to formulate the problem of simulating the two-dimensional (cross-sectional) transient movement of water and solute in saturated or partially saturated nonuniform porous media. The numerical model utilizes linear triangular elements. Nonreactive, as well as reactive solutes whose behaviour can be described by a distribution coefficient or first-order reaction term were considered. The flow portion of the model was tested by comparison of the model results with experimental and finite-difference results for transient flow in an unsaturated sand column and the solute transport portion of the model was tested by comparison with analytical solution results. The model was applied to a hypothetical case involving movement of water and solutes in tile-drained soils. The simulation results showed the development of distinct solute leaching patterns in the soil as drainage proceeded. Although applied to a tile drainage problem in this study, the model should be equally useful in the study of a wide range of two-dimensional water and solute migration problems.     

基于LS-RSGFD方法优化的横向各向同性(TI)介质一阶qP波高精度数值模拟

为克服各向异性介质弹性波数值模拟中存在着计算量大和波场分离困难等局限,研究了声学近似的VTI介质和TTI介质一阶qP波数值模拟方法.首先对VTI介质弹性波方程进行声学近似,推导了VTI介质一阶qP波方程;然后基于精确的TTI介质频散关系,引入一个包含各向异性控制参数σ的新辅助波场,推导了稳定的TTI介质二阶耦合qP波波动方程,并通过引入波场的伪速度分量,推导了等价的一阶应力-速度形式.结合旋转交错网格有限差分(RSGFD)和基于最小二乘优化的有限差分(LS-FD)两种各具优势的方法,研究了最小二乘旋转交错网格有限差分(LS-RSGFD)方法,并用其数值求解VTI和TTI介质一阶qP波方程,然后通过构造其LS-RSGFD格式,实现了高精度的各向异性介质qP波波场数值模拟.数值模拟结果表明:TI介质一阶qP波方程能够准确地模拟各向异性介质中qP波的运动学特征,引入控制参数σ能够有效地减弱不稳定性问题,保证非均匀TTI介质中qP波场的稳定传播;利用优化的LS-RSGFD方法可以得到高精度的合成地震记录,同时还可以相对地提高计算效率.     

Vulnerability of three spatially variable soil groups for solute leaching

Solute leaching in unsaturated soil is influenced by the variability in hydraulic functions (water retention and conductivity) that govern the flow process. Variability in measured soil hydraulic functions of a coarse-, medium- and fine-textured soil group was quantified with the scaling theory of similar media. Solute leaching in these soils was calculated with Monte Carlo simulation assuming, successively, hydraulic conductivity, K, volumetric water content, 0, and pressure head, h, to be constant. In addition to variability in hydraulic functions, variability in the solute retardation factor was also taken into account. To examine this effect five solutes were considered: a conservative solute (chloride), a non-retarded solute subject to decay (nitrate), a retarded solute that does not decay (cadmium) and two organic solutes which are retarded but have different sorption and decay parameters (the pesticide atrazine and a chlorinated hydrocarbon). The numerical results obtained with Monte Carlo simulation were in a number of instances verified with analytical solutions. The three soil groups distinguished showed considerable differences in vulnerability for leaching of the five solutes, emphasizing the importance of the effect of variability in soil hydraulic functions when studying solute leaching. Numerical and analytical results showed good agreement. Therefore, in relatively simple situations analytical solutions are attractive. However, in complicated situations, analytical solutions are cumbersome and numerical solutions are the only realistic alternative.     

Analytical Solutions for Water Flow and Solute Transport in the Unsaturated Zone

Analytical solutions for the water flow and solute transport equations in the unsaturated zone are presented. We use the Broadbridge and White nonlinear model to solve the Richards’ equation for vertical flow under a constant infiltration rate. Then we extend the water flow solution and develop an exact parametric solution for the advection-dispersion equation. The method of characteristics is adopted to determine the location of a solute front in the unsaturated zone. The dispersion component is incorporated into the final solution using a singular perturbation method. The formulation of the analytical solutions is simple, and a complete solution is generated without resorting to computationally demanding numerical schemes. Indeed, the simple analytical solutions can be used as tools to verify the accuracy of numerical models of water flow and solute transport. Comparison with a finite-element numerical solution indicates that a good match for the predicted water content is achieved when the mesh grid is one-fourth the capillary length scale of the porous medium. However, when numerically solving the solute transport equation at this level of discretization, numerical dispersion and spatial oscillations were significant.     

On approximation for fractional stochastic partial differential equations on the sphere

This paper gives the exact solution in terms of the Karhunen–Loève expansion to a fractional stochastic partial differential equation on the unit sphere \({\mathbb {S}}^{2} \subset {\mathbb {R}}^{3}\) with fractional Brownian motion as driving noise and with random initial condition given by a fractional stochastic Cauchy problem. A numerical approximation to the solution is given by truncating the Karhunen–Loève expansion. We show the convergence rates of the truncation errors in degree and the mean square approximation errors in time. Numerical examples using an isotropic Gaussian random field as initial condition and simulations of evolution of cosmic microwave background are given to illustrate the theoretical results.     

利用PETSc和FEPG编制海啸数值模拟程序的研究

海啸的数值模拟是海啸研究的一个重要领域, 它对于帮助理解海啸的基本物理特性和预防减灾具有重要意义。 海啸数值模拟程序的编制是一项繁杂的工作, 该文介绍了利用两种功能强大的通用软件(PETSc和FEPG)来进行海啸数值模拟程序编写的研究。 PETSc和FEPG采用有限差分、 有限元、 有限体积等多种离散方法, 可以对数值问题给出稳定的求解。 该文采用了直角坐标系下的非线性浅水波方程作为海啸波传播的控制方程进行离散求解, 并将其结果与TUNAMI N1模型进行了比较, 表明这两种方法方便而且有效。     

Grid Refinement in Cartesian Coordinates for Groundwater Flow Models Using the Divergence Theorem and Taylor's Series

Grid refinement is introduced in a numerical groundwater model to increase the accuracy of the solution over local areas without compromising the run time of the model. Numerical methods developed for grid refinement suffered certain drawbacks, for example, deficiencies in the implemented interpolation technique; the non‐reciprocity in head calculations or flow calculations; lack of accuracy resulting from high truncation errors, and numerical problems resulting from the construction of elongated meshes. A refinement scheme based on the divergence theorem and Taylor's expansions is presented in this article. This scheme is based on the work of De Marsily (1986) but includes more terms of the Taylor's series to improve the numerical solution. In this scheme, flow reciprocity is maintained and high order of refinement was achievable. The new numerical method is applied to simulate groundwater flows in homogeneous and heterogeneous confined aquifers. It produced results with acceptable degrees of accuracy. This method shows the potential for its application to solving groundwater heads over nested meshes with irregular shapes.     

A quasi three-dimensional model for flow and transport in unsaturated and saturated zones: 1. Implementation of the quasi two-dimensional case

A quasi three-dimensional (QUASI 3-D) model is presented for simulating the subsurface water flow and solute transport in the unsaturated and in the saturated zones of soil. The model is based on the assumptions of vertical flow in the unsaturated zone and essentially horizontal groundwater flow. The 1-D Richards equation for the unsaturated zone is coupled at the phreatic surface with the 2-D flow equation for the saturated zone. The latter was obtained by averaging 3-D flow equation in the saturated zone over the aquifer thickness. Unlike the Boussinesq equation for a leaky-phreatic aquifer, the developed model does not contain a storage term with specific yield and a source term for natural replenishment. Instead it includes a water flux term at the phreatic surface through which the Richards equation is linked with the groundwater flow equation. The vertical water flux in the saturated zone is evaluated on the basis of the fluid mass balance equation while the horizontal fluxes, in that equation, are prescribed by Darcy law. A 3-D transport equation is used to simulate the solute migration. A numerical algorithm to solve the problem for the general quasi 3-D case was developed. The developed methodology was exemplified for the quasi 2-D cross-sectional case (QUASI2D). Simulations for three synthetic problems demonstrate good agreement between the results obtained by QUASI2D and two fully 2-D flow and transport codes (SUTRA and 2DSOIL). Yet, simulations with the QUASI2D code were several times faster than those by the SUTRA and the 2DSOIL codes.