Direct Hydraulic Parameter and Function Estimation for Diverse Soil Types Under Infiltration and Evaporation

A new computationally efficient direct method is applied to estimating unsaturated hydraulic properties during steady-state infiltration and evaporation at soil surface. For different soil types with homogeneous and layered heterogeneity, soil hydraulic parameters and unsaturated conductivities are estimated. Unlike the traditional indirect inversion method, the direct method does not require forward simulations to assess the measurement-to-model fit; thus, the knowledge of model boundary conditions (BC) is not required. Instead, the method employs a set of local approximate solutions to impose continuity of pressure head and soil water fluxes throughout the inversion domain, while measurements act to condition these solutions. Given sufficient measurements, it yields a well-posed system of nonlinear equations that can be solved with optimization in a single step and is thus computationally efficient. For both Gardner’s and van Genuchten’s soil water models, unsaturated hydraulic conductivities and pressure heads (including the unknown BC) can be accurately recovered. When increasing measurement errors are imposed, inversion becomes less accurate, but the solution is stable, i.e., estimation errors remain bounded. Moreover, when the unsaturated conductivity model is known, inversion can recover its parameters; if it is unknown, inversion can recover a nonparametric, piecewise continuous function to which soil parameters can be obtained via fitting. Overall, inversion accuracy of the direct method is influenced by (1) measurement density and errors; (2) rate of infiltration or evaporation; (3) variation of the unsaturated conductivity; (4) flow direction; (5) the number of soil layers.     

The analysis of groundwater flows in unconfined aquifers with nonuniform hydraulic conductivity

Horizontal groundwater flows in unconfined aquifers with horizontal lower boundaries can be found exactly by the seepage analysis that allows the hydraulic conductivity to vary in the vertical direction. The exact analysis of flows when the lower boundary of the aquifer is not a horizontal plane, requires the soil-water pressure on this boundary to be known, and this is not generally the case except in the situation of a freshwater aquifer overlying saline water fed from the sea. For aquifers with spatial variations of hydraulic conductivity in both the vertical and horizontal directions, the seepage analysis can be modified to give groundwater flows in situations where the hydraulic conductivity can be represented by the product of independent functions of the three spatial coordinates. Different forms of three-dimensional variation are generated from suitably chosen functions. The use of such forms in calculations of equivalent uniform hydraulic conductivities of some groundwater flow regions demonstrates the dependence of equivalent hydraulic conductivity values on the flow boundary conditions. The exact groundwater flows calculated for particular groundwater situations by the seepage analysis provide results that are useful in validating numerical methods for solving groundwater problems in heterogeneous soils.     

New Approaches to the Propagation of Nonlinear Transients in Porous Media

Many problems in regional groundwater flow require the characterization and forecasting of variables, such as hydraulic heads, hydraulic gradients, and pore velocities. These variables describe hydraulic transients propagating in an aquifer, such as a river flood wave induced through an adjacent aquifer. The characterization of aquifer variables is usually accomplished via the solution of a transient differential equation subject to time-dependent boundary conditions. Modeling nonlinear wave propagation in porous media is traditionally approached via numerical solutions of governing differential equations. Temporal or spatial numerical discretization schemes permit a simplification of the equations. However, they may generate instability, and require a numerical linearization of true nonlinear problems. Traditional analytical solutions are continuous in space and time, and render a more stable solution, but they are usually applicable to linear problems and require regular domain shapes. The method of decomposition of Adomian is an approximate analytical series to solve linear or nonlinear differential equations. It has the advantages of both analytical and numerical procedures. An important limitation is that a decomposition expansion in a given coordinate explicitly uses the boundary conditions in such axis only, but not necessarily those on the others. In this article we present improvements of the method consisting of a combination of a partial decomposition expansion in each coordinate in conjunction with successive approximation that permits the consideration of boundary conditions imposed on all of the axes of a transient multidimensional problem; transient modeling of irregularly-shaped aquifer domains; and nonlinear transient analysis of groundwater flow equations. The method yields simple solutions of dependent variables that are continuous in space and time, which easily permit the derivation of heads, gradients, seepage velocities and fluxes, thus minimizing instability. It could be valuable in preliminary analysis prior to more elaborate numerical analysis. Verification was done by comparing decomposition solutions with exact analytical solutions when available, and with controlled experiments, with reasonable agreement. The effect of linearization of mildly nonlinear saturated groundwater equations is to underestimate the magnitude of the hydraulic heads in some portions of the aquifer. In some problems, such as unsaturated infiltration, linearization yields incorrect results.     

Modeling full-tensor anisotropy in groundwater flow via an iterative scheme for mixed finite elements

This paper presents an iterative scheme for the efficient simulation of groundwater flow in a two-dimensional, heterogeneous aquifer in which the hydraulic conductivity is anisotropic. The scheme is applicable to matrix equations arising from both mixed finite-element and cell-centered finite-difference approximations to the flow equations, and it extends readily to three space dimensions. The scheme, which generalizes an earlier technique for isotropic aquifer, admits a fast multigrid solver for hydraulic heads. Numerical experiments illustrate both the effectiveness of the scheme and the importance of accurately treating anisotropy: Small changes in the off-diagonal terms in the conductivity tensor cause relatively large changes in both the predicted heads and the Darcy velocities.     

Integrating Hydrogeochemical and Geophysical Data for Testing a Finite Volume Based Numerical Model for Saltwater Intrusion

In this paper, hydrogeological and geophysical data are used to validate a numerical model developed to predict seawater intrusion into coastal aquifers. The cell-centered finite volume method is adopted here to solve the set of coupled partial differential equations describing the motion of saltwater and freshwater separated by a sharp interface. These equations are based on the Dupuit approximation and are obtained from integration of 3D flow equations for fresh and salt water zones over the vertical dimension. In order to have flexibility upon complex configurations domain, non structured grid meshing is utilized. To approximate the diffusion fluxes, Green-Gauss type reconstruction, based on diamond-cell and least squares interpolation, is performed. The model is first validated using academic test case studies with known closed form solutions. The mathematical model has been calibrated using hydrogeochemical and geophysical data. The geophysical method applied in this study has been a frequency domain electromagnetic method. In this method the apparent electrical conductivity is measured by induction using two separate hand-held transmitter and receiver coils. During the operation the transmitter coil is energized by a low frequency alternating current that radiates an electromagnetic field and the receiver coil detects the resulting field. Taking into account the relationship between the bulk conductivity of the subsoil and the conductivity of groundwater, EM soundings have been interpreted to provide complementary information to hydrogeochemical data to outline the fresh–saltwater interface. This methodology has been applied to the case of saltwater intrusion into the Llobregat delta aquifer, near Barcelona, Spain.     

Calculation of Groundwater Flows in Coastal Pressurized Reservoirs

A model of a fresh groundwater flow through a rectangular horizontal pressurized reservoir toward a salt-water sea (basin, reservoir, trench, etc.) is examined within the framework of two-dimensional steady-state flow theory. In order to study the model, a mixed multi-parameter boundary value problem of the theory of analytic functions is formulated and solved using the Polubarinova-Kochina method. The structure and the characteristic features of the simulated process and the effect of all the determining physical parameters of the model on the nature of the flow are analyzed using the analytic dependences obtained and numerical calculations. An approximate hydraulic solution of the problem is compared with the exact solution obtained.     

Modeling the Flow Pattern at the Fractured Granites in Porto Alegre, Brazil

The studied 25 by 30 km site of Porto Alegre Southern Brazil fractured granites, accounts for surface lineaments assumed to mimic its fractures. Observations address the lineaments spatial characteristics, significant deviations between stationary groundwater levels condensed at part of the domain while sparse at most of the area, and no pumping/injection rates. Fractured media feature, e.g., dead end pathways, lack of flow interconnections, and preferential infiltration paths. These characteristics are not in line with the implementation of continuum-based macroscopic balance equations. Subject to such data and the inherent features, the objective was to verify if with lumped parameter modeling (LPM) following a flow directed graph approach we can assess a groundwater flow pattern. We at first addressed the site lineaments layout for which the evaluation of the hydraulic heads revealed the existence of isolated lineament clusters, leading to flow or no-flow zones. Aiming at better spatial distribution of flow distributions and based on the site surface lineaments, we established a virtual fracture network (VFN) for which the domain was subdivided into representative elementary area (REA) cells. Each REA was chosen so that the ratio of all its lineament lengths over that of the cell area remained practically unchanged between two consecutive subdivisions. Within each REA, lineaments with similar geometrical characteristics were considered as segment groups. The VFN was established upon elongating segments that created intersection points with other stretched segments from cells at the circumference of a considered REA. The evaluated steady state hydraulic head was compared between two LPM solutions: (1) Referring to flow along the VFN branches between intersection points, and (2) Using a flux interconnected network (FIN) for which each REA was replaced by a pole communicating flow to other poles. Computation of the FIN approach was significantly less intense. Both of these approaches resolved with hydraulic head isolines consistently similar to those obtained by interpolation between the observed groundwater levels. One reported event of a Nitrate polluted well and its plausible contaminating source within the study domain, show that it is in line with the predicted resultant flow direction following the FIN map.     

Segmented multigrid domain decomposition procedure for incompressible viscous flows

The application of grid stretching or grid adaptation is generally required in order to optimize the distribution of nodal points for fluid-dynamic simulation. This is necessitated by the presence of disjoint high gradient zones, that represent boundary or free shear layers, reversed flow or vortical flow regions, triple deck structures, etc. A domain decomposition method can be used in conjunction with an adaptive multigrid algorithm to provide an effective methodology for the development of optimal grids. In the present study, the Navier-Stokes (NS) equations are approximated with a reduced Navier-Stokes (RNS) system, that represents the lowest-order terms in an asymptotic Re expansion. This system allows for simplified boundary conditions, more generality in the location of the outflow boundary, and ensures mass conservation in all subdomain grid interfaces, as well as at the outflow boundary. The higher-order (NS) diffusion terms are included through a deferred corrector, in selected subdomains, when necessary. Adaptivity in the direction of refinement is achieved by grid splitting or domain decomposition in each level of the multigrid procedure. Normalized truncation error estimates of key derivatives are used to determine the boundaries of these subdomains. The refinement is optimized in two co-ordinate directions independently. Multidirectional adaptivity eliminates the need for grid stretching so that uniform grids are specified in each subdomain. The overall grid consists of multiple domains with different meshes and is, therefore, heavily graded. Results and computational efficiency are discussed for the laminar flow over a finite length plate and for the laminar internal flow in a backward-facing step channel.     

Discrete non‐local absorbing boundary condition for exterior problems governed by Helmholtz equation

The finite element method is employed to approximate the solutions of the Helmholtz equation for water wave radiation and scattering in an unbounded domain. A discrete, non‐local and non‐reflecting boundary condition is specified at an artificial external boundary by the DNL method, yielding an equivalent problem that is solved in a bounded domain. This procedure formulates a boundary value problem in a bounded region by imposing a relation in the discrete medium between the nodal values at the two last layers. For plane geometry, this relation can be found by straightforward eigenvalue decomposition. For circular geometry, the plane condition is applied at the external layer and this condition is condensed through a structured annular region, resulting in a condition at an inner radius. Exterior problems with a bounded internal physical obstacle are considered. It is well‐known that these kind of problems are well‐posed, and have a unique solution. Numerical studies based on standard Galerkin methodology examine the dependence of the DNL condition with respect to the circular annular region width. The DNL condition is compared with local boundary conditions of several orders. Numerical examples confirm the important improvement in accuracy obtained by the DNL method over standard conditions. Copyright © 1999 John Wiley & Sons, Ltd.     

Numerical simulation of high‐Reynolds number flow around circular cylinders by a three‐step FEM–BEM model

An innovative computational model, developed to simulate high‐Reynolds number flow past circular cylinders in two‐dimensional incompressible viscous flows in external flow fields is described in this paper. The model, based on transient Navier–Stokes equations, can solve the infinite boundary value problems by extracting the boundary effects on a specified finite computational domain, using the projection method. The pressure is assumed to be zero at infinite boundary and the external flow field is simulated using a direct boundary element method (BEM) by solving a pressure Poisson equation. A three‐step finite element method (FEM) is used to solve the momentum equations of the flow. The present model is applied to simulate high‐Reynolds number flow past a single circular cylinder and flow past two cylinders in which one acts as a control cylinder. The simulation results are compared with experimental data and other numerical models and are found to be feasible and satisfactory. Copyright © 2001 John Wiley & Sons, Ltd.     

基于区域分解算法的地下水耦合模型及其应用

根据岩体中不同区域裂隙发育规模的差异分为连续区域和离散区域,在不同区域中使用不同的地下水运动数学模型,应用区域分解算法来解决这类问题。其中连续区域采用了等效连续介质模型,离散区域采用了随机裂隙网络模型,通过区域公共边界上水位和流量连续的条件将两模型耦合求解。将基于区域分解算法的耦合模型应用于锦屏水电站坝址区三维渗流场的模拟中,通过钻孔观测水位和计算水位的对比发现,该方法是有效的,能够应用于实际工程。     

Modeling Three-Dimensional Groundwater Flows by the Body-Fitted Coordinate (BFC) Method: II. Free and Moving Boundary Problems

This paper presents a methodology and solution procedure of the time-dependent body-fitted coordinate (BFC) method for the analysis of transient, three-dimensional groundwater flow problems characterized by free and moving boundaries. The technique consists of numerical grid generation, time-dependent body-fitted coordinate transformation, and application of the finite difference method (FDM) to the transformed partial differential equations. Based on the time-dependent BFC method, a three-dimensional finite-difference computer code, BFC3DGW, was developed and used to solve two unconfined flow problems. The code was verified by comparing numerical results with analytical solutions for a steady-state seepage problem. In order to demonstrate capability of the method in dealing with flow problems with irregular and moving boundary surfaces, an unconfined well-flow problem was solved by the developed code. Difficulties associated with the free and moving irregular boundary have been successfully overcome by employing this method.     

Toward low‐noise synthetic turbulent inflow conditions for aeroacoustic calculations

The accuracy of boundary conditions for computational aeroacoustics is a well‐known challenge, due in part to the necessity of truncating the flow domain and replacing the analytical boundary conditions at infinity with numerical boundary conditions. In particular, the inflow boundary condition involving turbulent velocity or scalar fields is likely to introduce spurious waves into the domain, therefore degrading the flow behavior and deteriorating the physical acoustic waves. In this work, a method to generate low‐noise, divergence‐free, synthetic turbulence for inflow boundary conditions is proposed. It relies on the classical view of turbulence as a superposition of random eddies convected with the mean flow. Within the proposed model, the vector potential and the requirement that the individual eddies must satisfy the linearized momentum equations about the mean flow are used. The model is tested using isolated eddies convected through the inflow boundary and an experimental benchmark data for spatially decaying isotropic turbulence. Copyright © 2013 John Wiley & Sons, Ltd.     

Fluid flow in a porous medium containing partially closed fractures

The model of Snow, in which a fracture is represented by two parallel channel walls, has frequently been used to study the flow of fluid in fractured reservoirs. Although this model gives important insight into the flow in fractures, very few naturally occurring fractures have smooth parallel faces. In this paper, a simple model of partially contacting and en-echelon fractures frequently found in geological materials is presented. In this model, a fracture is viewed as a planar region where separation and contact zones both exist. To analyse the fluid flow in a porous medium containing fractures of this type, a planar array of periodically spaced fracture segments is analysed. The flow through a single fracture is deduced by taking the limit as the spacing between neighbouring fractures becomes large. The hydraulic conductivity parallel to the fractures is found to be the parallel combination of the conductivity of the porous matrix and the system of parallel fractures, the individual fracture conductance being a series combination of the hydraulic conductance of the separation and contact zones. This interpretation enables the conductance of the contact zones to be evaluated and the results to be generalised to the case in which the material in the contact regions has a hydraulic conductivity different to that of the matrix. This may arise, for example, from grain-size reduction during fracturing or may result from a partial mineralisation or cementation of the fracture.     

Rotational inviscid flow with circulation zones in a channel of variable cross section

The development of the theory of rotational motion of inviscid fluids for the purposes of describing channel flow encounters certain difficulties in connection with the appearance of viscosity effects near the walls. In the potential-rotational model [1], in which the vorticity is nonzero only in a closed circulation zone surrounded by potential flow, it is assumed that the separation and attachment points are known in advance. For example, for flow around a cavity these points coincide with the extreme corner points of the contour. The problem of determining the vorticity in a closed zone for the potential-rotational model has been investigated in a number of studies [2, 3], etc. In the case of an incompressible fluid the vorticity in the circulation zone is constant for two-dimensional flow and proportional to the distance from the axis for axisymmetric flow. The value of the constant is found from the steady-state condition for the adjoining viscous layers. If the channel walls have a smooth profile without corner points, then for determining the boundaries of the circulation zones additional conditions must be used. This study employs another scheme, in which the vorticity is formed outside the region of flow and in a particular problem is specified in the form of a boundary condition. An analytic solution describing the rotational flow of an inviscid fluid in a channel with a slightly varying cross section is obtained. Three types of entrance flow nonuniformity are considered: 1) uniform shear flow, 2) wake-type flow, and 3) potential flow with a narrow wall boundary layer. Streamline patterns with circulation zones are constructed for flows in diffuser channels with the above-mentioned types of entrance nonuniformity. A model of flow separation in a channel with a turbulent boundary layer on the walls is discussed.Translated from Izvestiya Akademii Nauk SSSR, Mekhanika Zhidkosti i Gaza, No. 2, pp. 31–37, March–April, 1985.In conclusion the author wishes to thank E. Yu. Shal'man, A. N. Kraiko, and A. B. Vatazhin for useful discussions and advice.     

Infiltration of Salt Solute in Homogeneous and Saturated Porous Media – An Analytical Solution Evaluated by Numerical Simulations

In many cases various land disposal activities (e.g. infiltration, injection wells) constitute an important potential source of groundwater contamination. Using a 2D physical model, the behaviour of the infiltration of a salt solute, locally injected in a homogeneous and saturated porous medium, has been analysed. Under various experimental conditions (density effects, injection flow rate) the salt solute penetrates the porous media and leads to a steady-state regime inside the mixing zone. By using experimental observations, the basic equations describing the flow and transport phenomena can be simplified and an analytical solution obtained. Its validity is subject to numerical verification. The numerical model, based on the development of the mass balance equation expressed by its conservative form, uses a combination of the mixed hybrid finite element (MHFE) and discontinuous finite element (DFE) methods. The efficiency of this numerical model was previously verified on standard benchmarks, for example Elder's problem and Henry's problem. In the first step, the qualitative good agreement between the experimental and numerical results enabled us to use the numerical model in order to verify some hypotheses resulting from visual observations. Thus, the numerical results reveal the existence of a steady-state regime inside the mixing zones. Nevertheless, both its vertical and longitudinal extensions are less than those observed in the physical model. In the second step, the numerical results enable to establish the validity domain as well as the accuracy of the proposed analytical solution.     

Two theorems on the steady-state conduction in anisotropic body

In this paper, steady-state plane conduction problems are considered. Two theorems are proven on the conductance of non-homogeneous and anisotropic plane conductors. Here, the conductor is a two-dimensional bounded plane domain with a conducting matter having two seperated boundary terminals and two seperated insulated boundary segments. Known theorems referring to homogeneous and isotropic conductors have been extended to non-homogeneous and anisotropic ones. The results of the paper can be directly used in following fields: heat flow, diffusion, irrotational hydraulic flow, flow of electrical current and anti-plane shear deformation. Application of formulae derived is illustrated in the examples of heat flow and anti-plane elastic shear deformation.     

Stochastic Flow Simulation and Particle Transport in a 2D Layer of Random Porous Medium

A stochastic numerical method is developed for simulation of flows and particle transport in a 2D layer of porous medium. The hydraulic conductivity is assumed to be a random field of a given statistical structure, the flow is modeled in the layer with prescribed boundary conditions. Numerical experiments are carried out by solving the Darcy equation for each sample of the hydraulic conductivity by a direct solver for sparse matrices, and tracking Lagrangian trajectories in the simulated flow. We present and analyze different Eulerian and Lagrangian statistical characteristics of the flow such as transverse and longitudinal velocity correlation functions, longitudinal dispersion coefficient, and the mean displacement of Lagrangian trajectories. We discuss the effect of long-range correlations of the longitudinal velocities which we have found in our numerical simulations. The related anomalous diffusion is also analyzed.     

Experimental Determination of Hydraulic Conductivity of Pine Samples in the Longitudinal Direction During Convective Drying

In studying the process of drying wood samples, the hydraulic conductivity of pine samples exposed to a convective air flow was measured in the direction perpendicular to the cross-sectional plane of the tree (along the wood grain). The kinetic curves of the drying process were treated with a technique based on the approximate solution of the one-dimensional diffusion equation for hydraulic conductivity of wood with a boundary condition of the third kind. The technique was tested on the basis of the known value of hydraulic conductivity in the direction tangential to annual rings of the tree. It is shown that the hydraulic conductivity in the longitudinal direction is 17 times the hydraulic conductivity in the direction tangential to the annual rings in the cross-sectional plane of the tree. By means of numerical simulation of the process, based on solving the initial-boundary problem for the two-dimensional linear equation of hydraulic conductivity, the effect of anisotropy of hydraulic conductivity coefficients on the dependence of the mean humidity on time and on the local humidity distribution is studied.     

MODELLING AND CALCULATION OF THE TURBULENT BOUNDARY LAYER ON A ROTATING CYLINDER SURFACE WITH STRONG SUCTION 1)

提出了湍流边界层的一种简单、快速计算方法, 用以求解强吸气作用下旋转圆筒表面边界层流动. 首先, 理论分析了同心圆筒间的旋转流体运动, 外筒静止、内筒旋转且为多孔吸气条件. 强吸气情况下旋转流动主要表现为内筒壁面附近的边界层流动, 基于这一事实得到了周向速度分布的解析表达式. 其次, 通过引入新参数扩展Cebeci-Smith代数湍流模型, 使其能考虑流线曲率、壁面吸气、低Reynolds数效应等因素. 针对这些因素的综合影响, 采用解析修正和经验参数对模型进行调整. 同时, 基于Reynolds应力湍流模型的仿真结果, 校准代数湍流模型中的经验参数. 最后, 给出基于广义Cebeci-Smith湍流模型的旋转壁面边界层流动的迭代算法, 该算法适用于需要特殊迭代过程的轴向及周向流动均匀情况. 计算了不同旋转速度和吸气强度组合工况下的边界层流动, 其周向速度和湍流强度分布与基于Reynolds应力湍流模型的计算结果非常接近. 并且表明, 当Reynolds应力湍流模型数值模拟预测内筒边界层为稳定层流时, 该方法也再现了相同初始条件下的层流边界层.