Finite element analysis of three-dimensional groundwater flow problems

A numerical method is presented for analysing either steady state or transient three-dimensional groundwater flow problems. The governing equation is formulated in terms of the finite element process using the Galerkin approach, and cubic isoparametric elements are used to simulate the flow domain as these permit accurate modelling of curved boundaries. Particular attention is paid to the time dependent movement of the phreatic surface where an iterative technique based on the replacement of the original transient problem by a discrete number of steady state problems is used to effect a solution. Furthermore, in tracing the movement of the surface use is made of the element formulation theory in order to compute the normal to the boundary.The validity of the technique is first established by analysing a radially symmetrical problem for which an alternative analytical solution is available. Finally, a general three-dimensional flow system is studied for which there is no known analytical solution. It is shown that relatively few elements are required to yield practical solutions.     

Galerkin finite element method and the groundwater flow equation: 1. Convergence of the method

This paper is concerned with the convergence of the Galerkin finite element method applied to a groundwater flow problem containing a borehole, with special reference to quadrature effects and the accuracy of the solution. It is shown that there exists an optimal quadrature rule for every choice of piecewise polynomial basis functions. Another interesting result proved here is that, in a direct application of the method the accuracy is very nearly independent of the degree of the polynomial basis functions, but strongly dependent on the distance of the borehole from the boundary if this is small.     

Coupled groundwater flow and transport: 1. Verification of variable density flow and transport models

This work examines variable density flow and corresponding solute transport in groundwater systems. Fluid dynamics of salty solutions with significant density variations are of increasing interest in many problems of subsurface hydrology. The mathematical model comprises a set of non-linear, coupled, partial differential equations to be solved for pressure/hydraulic head and mass fraction/concentration of the solute component. The governing equations and underlying assumptions are developed and discussed. The equation of solute mass conservation is formulated in terms of mass fraction and mass concentration. Different levels of the approximation of density variations in the mass balance equations are used for convection problems (e.g. the Boussinesq approximation and its extension, fully density approximation). The impact of these simplifications is studied by use of numerical modelling.Numerical models for nonlinear problems, such as density-driven convection, must be carefully verified in a particular series of tests. Standard benchmarks for proving variable density flow models are the Henry, Elder, and salt dome (HYDROCOIN level 1 case 5) problems. We studied these benchmarks using two finite element simulators - ROCKFLOW, which was developed at the Institute of Fluid Mechanics and Computer Applications in Civil Engineering and FEFLOW, which was developed at the Institute for Water Resources Planning and Systems Research Ltd. Although both simulators are based on the Galerkin finite element method, they differ in many approximation details such as temporal discretization (Crank-Nicolson vs predictor-corrector schemes), spatial discretization (triangular and quadrilateral elements), finite element basis functions (linear, bilinear, biquadratic), iteration schemes (Newton, Picard) and solvers (direct, iterative). The numerical analysis illustrates discretization effects and defects arising from the different levels of the density of approximation. We contribute new results for the salt dome problem, for which inconsistent findings exist in literature. Applications of the verified numerical models to more complex problems, such as thermohaline and three-dimensional convection systems, will be presented in the second part of this paper.     

Lanczos method for the solution of groundwater flow in discretely fractured porous media

One of the more advanced approaches for simulating groundwater flow in fractured porous media is the discrete-fracture approach. This approach is limited by the large computational overheads associated with traditional modeling methods. In this work, we apply the Lanczos reduction method to the modeling of groundwater flow in fractured porous media using the discrete-fracture approach. The Lanczos reduction method reduces a finite element equation system to a much smaller tridiagonal system of first-order differential equations. The reduced system can be solved by a standard tridiagonal algorithm with little computational effort. Because solving the reduced system is more efficient compared to solving the original system, the simulation of groundwater flow in discretely fractured media using the reduction method is very efficient. The proposed method is especially suitable for the problem of large-scale and long-term simulation. In this paper, we develop an iterative version of Lanczos algorithm, in which the preconditioned conjugate gradient solver based on ORTHOMIN acceleration is employed within the Lanczos reduction process. Additional efficiency for the Lanczos method is achieved by applying an eigenvalue shift technique. The “shift” method can improve the Lanczos system convergence, by requiring fewer modes to achieve the same level of accuracy over the unshifted case. The developed model is verified by comparison with dual-porosity approach. The efficiency and accuracy of the method are demonstrated on a field-scale problem and compared to the performance of classic time marching method using an iterative solver on the original system. In spite of the advances, more theoretical work needs to be carried out to determine the optimal value of the shift before computations are actually carried out.     

Parameter-independent model reduction of transient groundwater flow models: Application to inverse problems

A new methodology is proposed for the development of parameter-independent reduced models for transient groundwater flow models. The model reduction technique is based on Galerkin projection of a highly discretized model onto a subspace spanned by a small number of optimally chosen basis functions. We propose two greedy algorithms that iteratively select optimal parameter sets and snapshot times between the parameter space and the time domain in order to generate snapshots. The snapshots are used to build the Galerkin projection matrix, which covers the entire parameter space in the full model. We then apply the reduced subspace model to solve two inverse problems: a deterministic inverse problem and a Bayesian inverse problem with a Markov Chain Monte Carlo (MCMC) method. The proposed methodology is validated with a conceptual one-dimensional groundwater flow model. We then apply the methodology to a basin-scale, conceptual aquifer in the Oristano plain of Sardinia, Italy. Using the methodology, the full model governed by 29,197 ordinary differential equations is reduced by two to three orders of magnitude, resulting in a drastic reduction in computational requirements.     

Application of the discontinuous spectral Galerkin method to groundwater flow

The discontinuous spectral Galerkin method uses a finite-element discretization of the groundwater flow domain with basis functions of arbitrary order in each element. The independent choice of the basis functions in each element permits discontinuities in transmissivity in the flow domain. This formulation is shown to be of high order accuracy and particularly suitable for accurately calculating the flow field in porous media. Simulations are presented in terms of streamlines in a bidimensional aquifer, and compared with the solution calculated with a standard finite-element method and a mixed finite-element method. Numerical simulations show that the discontinuous spectral Galerkin approximation is more efficient than the standard finite-element method (in computing fluxes and streamlines/pathlines) for a given accuracy, and it is more accurate on a given grid. On the other hand the mixed finite-element method ensures the continuity of the fluxes at the cell boundaries and it is particular efficient in representing complicated flow fields with few mesh points. Simulations show that the mixed finite-element method is superior to the discontinuous spectral Galerkin method producing accurate streamlines even if few computational nodes are used. The application of the discontinuous Galerkin method is thus of interest in groundwater problems only when high order and extremely accurate solutions are needed.     

Finite element solutions for mass oscillations in surge tanks

This paper describes a finite element technique using the method of weighted residuals for the solution of mass oscillations in surge tanks. Three weighting functions, uniform, linear and Galerkin, are applied and the results are compared with those from alternative techniques. The relatively simple case of surge analysis with flow rate change in the penstock, but neglecting tunnel friction, is first considered as a direct analytical solution is available. Finally friction is included for comparison with a graphical and analogue solution.     

Recharge process and aquifer models of a small watershed

Abstract

Kanchanapally watershed covering an area of about 11 km² in Nalgonda district, Andhra Pradesh, India is located in granitic terrain. Groundwater recharge has been estimated from a water balance model using hydrometeorological data from 1978–1994. The monthly recharge estimates obtained from the water balance model formed input for the groundwater flow model during transient model testing. The groundwater flow model has been prepared to simulate steady state groundwater conditions of 1977 using the nested squares finite difference method. The transient groundwater flow model has been tested during 1977–1994 using the estimated recharge values. The present study helped verify the usefulness of monthly recharge estimates for accounting dynamic variations in recharge as reflected in water level fluctuations in hydrographs.     

Coupling local discontinuous and continuous galerkin methods for flow problems

In this paper, we discuss the local discontinuous Galerkin (LDG) method applied to elliptic flow problems and give details on its implementation, focusing specifically on the case of piecewise linear approximating functions. The LDG method is one a family of discontinuous Galerkin (DG) methods proposed for diffusion models. These DG methods allow for very general hp finite element meshes, and produce locally conservative fluxes which can be used in coupling flow with transport. The drawback to DG methods, when compared to their continuous counterparts, is the number of degrees of freedom required to compute the solution. This motivates a coupled approach, discussed herein, where the solution is allowed to be continuous or discontinuous on a node-by-node basis. This coupled approximation is locally conservative in regions where the numerical solution is discontinuous. Numerical results for fully discontinuous, continuous and coupled discontinuous/continuous solutions are given, where we compare solution accuracy, matrix condition numbers and mass balance errors for the various approaches.     

Comparison of Galerkin-type discretization techniques for two-phase flow in heterogeneous porous media

A comparison of Standard Galerkin, Petrov-Galerkin, and Fully-Upwind Galerkin methods for the simulation of two-phase flow in heterogeneous porous media is presented. On the basis of the coupled pressure-saturation equations, a generalized formulation for all three finite element methods is derived and analysed. For flow in homogeneous media, the Petrov-Galerkin method gives excellent results. But this method fails miserably for problems with heterogeneous media. This is because it is not able to capture correctly processes that take place at interfaces when, for instance, the capillary pressure-saturation relationship after Brooks and Corey is assumed. The Fully-Upwind Galerkin method is superior to the Petrov-Galerkin approach because it is able to give correct results for flow in homogeneous and heterogeneous media for the two models of van Genuchten and Brooks-Corey. The widely used formulation which is correct for the homogeneous case cannot be used for heterogeneous media. Instead the straightforward approach of gradp_c in combination with a chord-slope technique must be utilized.     

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.     

Nodal domain integration model of two-dimensional advection-diffusion processes

The nodal domain integration method is applied to a two-dimensional advection-diffusion process in an anisotropic inhomogeneous medium. The domain is discretised into the union of irregular triangle finite elements with vertex-located nodal points and a linear trial function is used to approximate the governing flow equation's state variable in each element. Non-linear parameters are assumed quasi-constant for small durations in time in each element. The resulting numerical model represents the Galerkin and subdomain integration weighted residual methods and the integrated finite difference method as special cases. Both Dirichlet and Neumann boundary conditions are accommodated in a manner similar to the Galerkin finite element approach.     

A 3-D finite element conjugate gradient model of subsurface flow with automatic mesh generation

The 3-D flow modelling of groundwater systems of realistic size generally requires a big effort for the preparation of the input data as well as large computational costs. A numerical finite element model (MAITHREE) is developed for the efficient analysis of the steady and unsteady behaviour of natural confined 3-D basins. Starting from an initial triangular grid the code automatically generates a set of tetrahedral elements in each of the geologic units or subunits specified by the user in the vertical profile. The original element incidences list is then rearranged in order to provide conforming 3-D elements throughout the domain. The model is designed with a view to saving much of the labour involved in setting a 3-D grid and to providing flexibility as well as economical convenience through a high computational efficiency. The latter task is achieved by the aid of a solver based on the modified conjugate gradient (MCG) method which has proved to be an excellent technique for the solution of large linear finite element sets of sparse 3-D subsurface equations. Some examples derived from both hypothetical and real-world situations are discussed to illustrate the innovative features of MAITHREE and its computational performance.     

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.     

Estimating Aquifer Properties Using Derivative Analysis of Water Level Time Series from Active Well Fields

Water level time series from groundwater production wells offer a transient dataset that can be used to estimate aquifer properties in areas with active groundwater development. This article describes a new parameter estimation method to infer aquifer properties from such datasets. Specifically, the method analyzes long‐term water level measurements from multiple, interacting groundwater production wells and relies on temporal water level derivatives to estimate the aquifer transmissivity and storativity. Analytically modeled derivatives are compared to derivatives calculated directly from the observed water level data; an optimization technique is used to identify best‐fitting transmissivity and storativity values that minimize the difference between modeled and observed derivatives. We demonstrate how the consideration of derivative (slope) behavior eliminates uncertainty associated with static water levels and well‐loss coefficients, enabling effective use of water level data from groundwater production wells. The method is applied to time‐series data collected over a period of 6 years from a municipal well field operating in the Denver Basin, Colorado (USA). The estimated aquifer properties are shown to be consistent with previously published values. The parameter estimation method is further tested using synthetic water level time series generated with a numerical model that incorporates the style of heterogeneity that occurs in the Denver Basin sandstone aquifers.     

Nonlinear dynamics of long-wave perturbations of the Kolmogorov flow for large Reynolds numbers

The nonlinear dynamics of long-wave perturbations of the inviscid Kolmogorov flow, which models periodically varying in the horizontal direction oceanic currents, is studied. To describe this dynamics, the Galerkin method with basis functions representing the first three terms in the expansion of spatially periodic perturbations in the trigonometric series is used. The orthogonality conditions for these functions formulate a nonlinear system of partial differential equations for the expansion coefficients. Based on the asymptotic solutions of this system, a linear, quasilinear, and nonlinear stage of perturbation dynamics is identified. It is shown that the time-dependent growth of perturbations during the first two stages is succeeded by the stage of stable nonlinear oscillations. The corresponding oscillations are described by the oscillator equation containing a cubic nonlinearity, which is integrated in terms of elliptic functions. An analytical formula for the period of oscillations is obtained, which determines its dependence on the amplitude of the initial perturbation. Structural features of the field of the stream function of the perturbed flow are described, associated with the formation of closed vortex cells and meandering flow between them. As a supplement, an asymptotic analysis of nonlinear dynamics of long-wave perturbations superimposed on a damped by small viscosity Kolmogorov flow (very large, but finite Reynolds numbers) is made. It is strictly shown that all velocity components of the perturbed flow remain bounded in this case.     

Groundwater dynamics in a till hillslope: flow directions,gradients and delay

Knowledge of groundwater dynamics is important for the understanding of hydrological controls on chemical processes along the water flow pathways. To increase our knowledge of groundwater dynamics in areas with shallow groundwater, the groundwater dynamics along a hillslope were studied in a boreal catchment in Southern Sweden. The forested hillslope had a 1‐ to 2‐m deep layer of sandy till above bedrock. The groundwater flow direction and slope were calculated under the assumption that the flow followed the slope of the groundwater table, which was computed for different triangles, each defined by three groundwater wells. The flow direction showed considerable variations over time, with a maximum variation of 75°. During periods of high groundwater levels the flow was almost perpendicular to the stream, but as the groundwater level fell, the flow direction became gradually more parallel to the stream, directed in the downstream direction. These findings are of importance for the interpretation of results from hillslope transects, where the flow direction usually is assumed to be invariable and always in the direction of the hillslope. The variations in the groundwater flow direction may also cause an apparent dispersion for groundwater‐based transport. In contrast to findings in several other studies, the groundwater level was most responsive to rainfall and snowmelt in the upper part of the hillslope, while the lower parts of the slope reached their highest groundwater level up to 40 h after the upper parts. This can be explained by the topography with a wetter hollow area in the upper part. Copyright © 2011 John Wiley & Sons, Ltd.     

Finite element analysis of water flow in variably saturated soil

A two-dimensional Galerkin finite element model for water flow in variably saturated soil is presented. A fourth-order Runge-Kutta time integration method is employed which allows use of time steps at least 2 times greater than for a traditional finite difference approximation of time derivatives. For short total simulation times computer execution costs for the Runge-Kutta method are greater than for the finite difference approximation due to the start up cost of the Runge-Kutta method, but for longer simulation times the Runge-Kutta method requires considerably less computational effort even when automatic time-step adjustment is used with the finite difference procedure. A comparison of the method of influence coefficients and 2 × 2 Gaussian integration to compute element matrices indicates that the influence coefficient method reduces total execution time to 60% of that required for numerical quadrature. Computed pressure heads using the influence coefficient method and numerical integration are found to be in close agreement with each other even under conditions of highly non-linear soil properties in a heterogeneous domain. Fluxes computed by the two methods are also generally in close agreement except under extremely non-linear conditions when some deviations were observed at short simulation times.     

Cyclic hysteretic flow in porous medium column: model, experiment, and simulations

A periodic vertical movement of the groundwater table results in a subsequent cyclic response of the water content and pressure profiles in the vadose zone. The sequence of periodic wetting and drying processes can be affected by hysteresis effects in this zone. A one-dimensional saturated/unsaturated flow model based on Richards’ equation and the Mualem (Soil Sci. 137 (1984) 283) hysteresis model is formulated which can take into account multi-cycle hysteresis effects in the relation between capillary pressure and water content. The numerical integration of the unsaturated flow equation is based on a Galerkin-type finite element method. The flow domain is discretised by finite elements with linear shape functions. Simulations start with static water content and pressure profiles, which correspond to either a boundary drying or wetting retention curve. To facilitate the numerical solution of the hysteretic case an implicit non-iterative procedure was chosen for the solution of the nonlinear differential equation. Laboratory experiments were performed with a vertical sand column by imposing a high frequency periodic pressure head at the lower end of the column. The total water volume in the column, and the periodic water content profile averaged over time were measured. The boundary drying and wetting curves of the relation between water content and capillary pressure were determined by independent experiments. The simulations of the experimental conditions show a clear effect of the hysteresis phenomenon on the water content profile. The simulations with hysteresis agree well with the measurements. Computed dimensionless water content profiles are presented for different oscillation frequencies with and without consideration of hysteresis.     

Finite difference modeling of contaminant transport using analytic element flow solutions

The two-dimensional implementation of the analytic element method (AEM) is commonly used to simulate steady-state saturated groundwater flow phenomena at regional and local scales. However, unlike alternative groundwater flow simulation methods, AEM results are not ordinarily used as the basis for simulation of reactive solute transport. The use of AEM-simulated flow fields is impeded by the discrepancy between a continuous representation of flow and a typically discrete representation of transport, and requires translation of the flow solution to a discrete analog. This paper presents a variety of methods for analytically calculating conservative discrete water fluxes and integrated components of the dispersion tensor across cell interfaces. An Eulerian finite difference method based on these AEM-derived parameters is implemented for use in simulation of 2D (vertically averaged) solute transport. This implementation is first benchmarked against existing methods that use standard finite difference flow solutions, then used to investigate the effects of an inaccurate discrete water balance. It is shown that improper translation of AEM fluxes leads to significant water balance errors and inaccurate simulation of contaminant transport.