17 Groundwater Hydrology IV (Coupled Flow and Transport)
Prabhas K Yadav
Objectives:
This module will present introductory topics relevant to flow and transport problems in groundwater. Emphasis will be on d processes, coupling them and eventually establishing the mathematical formulation of the problem. The target groups are higher level undergraduate students and the first year PG students. The module, as it is introductory, will only briefly introduce methods for solving transport problems. The specific objectives of the module are:
1. Derivation of the flow model
2. Recognizing factors and processes affecting transport problems
3. Derivation of transport models
4. Developing a systematic approach to solve transport problems.
17. Introduction
In this module we will get introduced to fundamentals of the flow and transport equations or models. To distinguish, flow deals with the quantity aspect of groundwater whereas transport deals with the quality. Hence the flow problem is governed by the energy gradient where as mass gradient is significant in the transport problem. The complexities in both problems arises from the properties of the porous media, such as porosity, conductivity etc, which are mostly varying in space (heterogeneity versus homogeneity) and direction (isotropic versus anisotropic). In this module we will first introduce the flow model and subsequently learn about different transport processes and develop a transport model. Emphasising on the transport problem we will introduce methods to solve it. Our approach will remain mostly mathematical as such we will attempt to solve and visualize few examples of transport problems. In this module we will focus on saturated (all voids filled), homogeneous (aquifer properties such as conductivity do not change along the extend) and isotropic (aquifer properties is uniform also in all directions) aquifers.
17.1 Flow Model
17.1.1 The Representative Control Volume (RCV)
Aquifers as we know already are extensive, at least in the horizontal and lateral directions, whereas its properties such as conductivity, transmissivity often changes in a very small scale, in many cases also at a pore-scale. As such to understand aquifer properties and underlying processes, we define a Representative Control Volume (RCV, see Fig. 1). Ideally a RCV is very much smaller than the aquifer that is investigated, and at the same it is very much larger than individual grain or pore. The essential condition is that Darcy’s law has to be applicable in the RCV. A big advantage of RCV is that in the aquifer processes analysis it can be oriented as per the coordinate system, e.g., rectangular for Cartesian coordinate.
17.1.2 The derivation of flow model
The groundwater flow problem (and hence the model) can be classified according to: the problem dimensionality (1D, 2D, or 3D), the time-dependent state (steady and transient), the variation in aquifer properties with respect to space (homogeneous and heterogeneous), the variation in aquifers properties with respect to direction (isotropic and anisotropic), the aquifer condition (confined and unconfined) and based on weather aquifer is with or without source and sink. Each of the above combination corresponds to a certain equation governing groundwater flow. The derivation of the equation requires two fundamental principles: the conservation of volume and the Darcy’s law, which is given as:
With q = Q/A i.e., the discharge Q [L3T-1] and Area A[L2]. vf is often referred to as Darcy’s velocity or specific discharge. K [LT-1] refers to hydraulic conductivity and gradh [-] is gradient vector of hydraulic head h [L]. The negative sign is indicator of the flow from larger to smaller head. With Darcy’s law already defined, we now develop a groundwater flow model for several scenarios. We begin with the derivation of 3D groundwater model.
We consider a unit 3D RCV (Fig. 2) in which vfx, vfy and vfz represents the specific discharge entering the RCV from the left face of the RCV. Since the RCV is very small, the change of the specific discharge across the RCV (at the right face) can be regarded as linear, i.e., we add the product of first-order derivative and the corresponding distance between faces. The volume budget is given by:
Next we find the relationship between the change in head and the corresponding change in water volume. This relation is expressed as
Where q [T-1] represents the volumetric rate of the source sink per unit aquifer volume. Few examples of q can be water injection/extraction through wells, water transferred to aquifer from/to rivers/streams. For anisotropic case the conductivity (K) becomes the function of spatial direction and the eq. (11) takes the following form:
With T and S defined, the 2D groundwater flow equation for a confined anisotropic aquifer with source/sink is
With source/sink term included, eq. (17) is also termed as Boussinesq equation. Further, eq. (17) is also valid for the unconfined aquifer, in which the top of the aquifer is the groundwater table, provided that Dupuit’s assumptions, that groundwater flow is horizontal, the flow velocity do not vary along the aquifer thickness, and that Darcy’s law (eq. 1) is valid at water table, are observed.
17.1.3 Complete formulation of flow model
The complete formulation of the flow model requires that we define the following in addition to an appropriate flow equation (eqs. 10 – 14 and 17):
1. Specify the geometric properties of the region of interest
2. Specify values of aquifer parameters (hydraulic conductivity, storage coefficient) by considering spatial variability and anisotropy, if necessary
3. Specify the initial condition (IC): Head values at time t=0 for the transient problems only.
4. Specify boundary conditions (BCs): Along the complete boundary, and which may be time-dependent. The three type of BCs used in groundwater flow problems are defined next.
4.1 Boundary condition of the first kind or Dirichlet condition: The head value is given.
4.2 Boundary condition of the second kind or Neumann condition: The component of the head gradient, which is perpendicular to the boundary, is provided.
4.3 Boundary condition of the third kind or Cauchy condition or Robin boundary condition: A relationship between the head value and the component of the head gradient, which is perpendicular to the boundary, is given.
17.2 Transport Processes
In the transport model we consider mass (solute, chemical) that is flowing along with the groundwater. Before we develop the transport model, we will first briefly learn about the different transport processes that affect the concentration of mass in the aquifer.
17.2.1. Advection
Advection (also called convection) is the transport of mass with the movement of a moving medium. The moving medium can be termed carrier, and in the groundwater case it is usually fluid. To quantify advection process, let us consider a laboratory column of a constant cross-sectional area, A [L2] and the steady-state discharge, Q [L3T-1]. Further we will assume ne [ ] the effective porosity of the column packing and v as average linear velocity of water flow. Note that v = q/ne, where q is the Darcy velocity (eq. 1). With these information, we quantify the advective mass flow rate Jadv[MT-1] as
J adv = QC (18)
We may want to modify eq. (18) by considering the continuity equation in the column
Q = ne Av = constant
Substituting eq. (19) in eq. (18), and defining advective flux, jadv[MT-1L2] = Jadv/A , we get
17.2.2 Mechanical Dispersion
The spread of mass causes dispersion. Spread of mass is not entirely a physical processes but physiochemical processes also causes spread. We refer mechanical dispersion to the spread due to physical processes only. Next, we quantify dispersive mass flow rate, Jdis[MT-1]. Standard hydrogeology texts, such as Domenico and Schwartz, (1998), Chahar (2014), suggest that dispersion can be explained as a random phenomena. Based on observations, the dispersive mass flow due to mechanical dispersion in porous media is found to be
with α [L] a constant of proportionality called dispersivity , ne [ -] is effective porosity, A [L2] is surface area of the flow, v [LT-1] is the groundwater flow velocity, L [L] is the flow length and C [ML-3]. The negative sign is used in equality relation to identify that transport is from a higher concentration towards the lower one. Dispersivity is a property of porous media alone and is not dependent on fluid type or the flow characteristics. In a more formal manner the ratio ΔC/L [ML-4] is termed as concentration gradient.
As mass transport is influenced by both porous medium and flow properties, a product of dispersivity and linear velocity called the mechanical dispersion coefficient, Dmech [L2T -1], is more often used in analysing dispersive transport. This changes eq. (21) to
with Dmech = αv. Instead of mass flow rate one may use mass flux jdis [MT-1L2] , which is
Since α characterizes spread length, its value is dependent on the flow coordinate direction under consideration. Advection and the mechanical dispersion are two processes that are generally used in analysing groundwater transport problems. Nevertheless, it is important for us to introduce a very important transport mechanism called diffusion that is not flow or mechanically driven but rather it is concentration gradient driven.
17.2.3 Diffusion
The quantification of diffusive mass transport, Jdif [MT-1] is based on the work in Fick, (1855). In that work the mass flow due to diffusion (1D) is described to be
with A [L2] is surface area of the flow, L [L] is the flow length and C [ML-3] and Ddif [L2T-1], is the constant of proportionality called diffusion coefficient. The negative sign in eq. (24) represents transport of mass from higher concentrations toward the lower one. From eq. (24) we can get the diffusive flux jdis [MT-1L2] from
It the analysis of groundwater transport problems the dispersion coefficient and the diffusion coefficient are summed and a new coefficient called hydrodynamic dispersion coefficient Dhyd is defined as in
This summation is based on the fact that both dispersion and diffusion lead to spread of mass, and can be justified only on practical basis. It is to be noted that dispersion is physically driven process whereas difference is chemically driven process.
17.3. Joint Action of advection and Dispersion
The transport of conservative (non-reacting) solute in an aquifer can be understood as a superposition of advection, mechanical dispersion and pore diffusion. To quantify the transported mass over the time, we sum all the mass flow rate, J [MT-1] that we have considered. i.e., eqs. (18, 22 and 24). Thus we get a combined equation for the mass flow
Where Jdis,h = Jdis + Jdif refer to mass flow due to hydrodynamic dispersion. Eq (27) can alternatively written as
The spreading of mass due to advection and dispersion can be quantified by combining a mass budget and the corresponding laws of motion. This result in a transport equation called advection-dispersion or convection-dispersion equation, which we derive in the next section.
17.4. Derivation of transport model
Transport equations are derived by using a representative control volume (or RCV see Fig 1). For the transport problems (2D), the RCV extends from the aquifer bottom to the aquifer top in confined aquifers. Furthermore average concentration along the thickness is used. The fundamental transport equation is based on mass balance (similar to eq 2, volume replaced by mass) relation and laws of motion used in advection and dispersion. For the derivation we will consider a 2D transport in the horizontal (xy-plane) with advection along the -axis. The groundwater flow is steady, i.e., vx is constant over time. we will consider 2D dispersion. The mass balance equation in the RCV is then (see Fig. 3)
Where ΔM is change in mass in the RCV over time Δt. For deriving the transport equation, we will have to find an expression for each term of eq. (29). Let ne be the porosity of RCV then we have:
I. The volume of the RCV, V = Δx Δy m
II. The dissolved mass, = ne V C
III. Thus total mass, M, in RCV is: M = ne VC
in which ne and V are fixed quantities, then the change of mass ΔM over time ΔM is given as
with ΔC being the change in solute concentration. Since the RCV is very small compared to the actual study area, the components of mass flux (j =J/A) can be assumed to change linearly across the extension of a RCV. These changes are obtained by multiplying first-order derivatives with corresponding distances, i.e., using Taylor series expansion (see Fig. 3) . We thus get
Note that we have used Dhyd = αv in the above equations. We will insert eqs. (36 and 37) to complete our transport equation (eq. 37), thus we get
Eq. (39), a 2D transport equation, can be converted to 1D transport equation by removing the y-components and likewise it can be converted to 3D transport equation by adding the z-components. Next we learn about other essential requirements for transport problems.
17.4.1 Initial and boundary conditions
A complete set of transport problem includes the transport equation, e.g., eq. (39), the initial conditions (time dependent) and the boundary conditions (space dependent). The initial condition describes the distribution of mass or specifies concentration in the entire area of investigation at some starting time, usually t=0 is considered. The initial conditions are required for the transport problems that are time dependent (e.g., eq. 39). These are so called transient problems and they comprise of general transport problems. A mathematical statement of initial condition is
C(t = 0) = C0 (40)
i.e., at t=0, C0 is the concentration in the entire investigated area.
For time-independent problems, i.e, ∂C/∂C = 0, also called steady-state problems boundary conditions are required to be specified. The boundary conditions quantify the impacts of condition at the surrounding just adjacent to the area under investigation. As the name suggests, the boundary conditions have to be specified for the entire boundary of the investigation area, even for the cases when the investigated area in unlimited or extend are infinite. The boundary conditions may or may not be time dependent. For transport problems three types of boundary conditions have been defined.
A first type of boundary condition, also called Dirichlet type, specifies the value of mass or concentration at the boundary. The mathematical statement of the first type boundary condition is
C(x, y, z, t) = C0 (41)
i.e., at any t, the concentration at specified space (x,y,z) is C0. Tracer injection, effluent concentration from polluted area can be few examples of first type boundary conditions. In a second type boundary condition, also called Neumann type condition, the component of the concentration gradient perpendicular to the boundary is specified. The mathematical statement of the second type boundary condition along x-axis is
Since the gradient of concentration is proportional to diffusive flux, the second type boundary condition is also called a flux specified boundary. Second type boundary conditions are also used for specifying the no-flow condition, i.e., no concentration gradient exist across the boundary. However it has to be noted that absolute no-flow can only exist when there is no concentration gradient and the velocity vector is zero. Rock formation at the bottom of aquifer will result to a no flow type of boundary condition.
The third type, Cauchy or Robin type, conditions combines the first two conditions, i.e. both a specific value and gradient is specified. This condition generalizes the other two boundary conditions. Mathematical statement of the third type boundary condition (along x-axis) is
17.5 Solving Transport Problems
In the previous sections we learned about the transport problems and developed the transport equation, also called Advection-Dispersion (AD) equation. In this section we will attempt to solve the AD model, which has been one of our objective from this module. The techniques for solving the transport problems are extremely varied and therefore we will restrict to simpler problems (1D, and few reaction components etc.). We will very briefly get introduced to complex problems (higher dimensions, multiple reaction).
Direct integrating of the AD equation (e.g. eq. 39) is only possible in very simplified cases (e.g., only advection case). However, transformation methods (Laplace, Fourier etc.) can be applied to directly solve ADR equation under some conditions. Thus obtained solutions are called analytical models, which are exact and often a closed-form. The conditions very generally include that aquifer is of regular geometry, is homogeneous and isotropic and the chemical reactions are linear models based. Numerical methods, e.g. Finite-Difference Method (FDM), Finitie-Element Method (FEM), that provide approximate solution of the ADR equation are very commonly used to solve transport problems. In this module we will restrict to an analytical model and very briefly discuss numerical methods.
17.5.1 Analytical Solution
The solution provided in the succeeding subsections generally assumes the following:
1. Fluid of constant density and viscosity
2. Flow in x-direction only, and velocity is constant.
3. The longitudinal dispersion coefficient Dx is constant.
In addition we assume aquifer is homogeneous and isotropic. Additional assumptions is provided when required.
17.5.2 The Adevetion-Dispersion (AD) equation
The simplest ADR equation can be a transient 1D equation without the reaction term, i.e. AD equation, which is
Readers should check the wikipedia site (https://en.wikipedia.org/wiki/Error_function) to find methods to approximate the value of erfc.
In order to visualize eq. (45) we perform a column (3 m long) experiment with C0= 0 mg/L, v =1 m/h and D = 1m2/h. Let our Cin = 1, 3 and 5 mg/L. We will check the concentration in the column at different times. Fig. 4 provides the results of our experiment. The result suggest that it will require 5 hours for our homogeneous distribution of mass in the column. For further visualisation of the Ogata and Banks solution, readers may check the GNU-Octave https://www.gnu.org/software/octave/) code (ogata1D.m) provided at https://github.com/prabhasyadav/UGC-Transport .
17.5.3 An introduction to numerical method for transport problems
Numerical modelling is generally used for solving transport problems. The main reason behind this is that analytical solution are possible for only few cases, often limited to regular domain geometry and linear (differential) equations. These limitations are easily overcome by numerical method. This topic is a very extensive. Here, we intend to only introduce the topic very briefly. Among few numerical methods that are available, Finite Difference Method (FDM) and Finite Element Method (FEM) are the two very extensively used in solving transport problems. We will focus on this two methods. Towards the end of the section we will use an example problem to illustrate FDM, the more common of the two methods mentioned above.
In the FDM the governing partial differential equation is replaced by a set of difference equations applicable to the system of nodes, i.e., the area under investigation is discretized into structured meshes. Taylor series expansion or polynomial fitting techniques are used to approximate all space/time derivative in terms of concentration. For example the first and second order derivatives of the transport equation is approximated using
where i refers to the node i and i+1 and i–1 are two adjacent nodes. Δx represents the width of the node, which may be a constant. A system of algebraic equation results when each node in the domain in considered. The system of equation is then solved using matrix algebra techniques. The most important aspect of FDM is the ease of formulating difference equations. Anderson and Woessner, (1992) provides an excellent introduction to the method. Readers are suggested to start with that literature if interested in learning the FDM for groundwater studies. For modelling, MODFLOW (an open-source) developed and maintained by United States Geological Survey (USGS) can be considered academic, research and industry standard. Further details can be obtained from http://water.usgs.gov/ogw/modflow/.
The mathematics of FEM is not as straightforward as that of FDM. In this method the area under investigation is subdivided into elements that are defined by nodes. The element can be of different shapes, although triangular and quadrilateral shapes are the most common one used. Generation of FEM mesh is very tedious, as such FEM codes usually include mesh generation software. The solution of the differential equation using FEM is found as a combination of shape functions. The shape function (also calledLagrangian function) for the triangular element are linear within each element. Thus in Cartesian coordinates the shape function f is
within element α. All coefficients aαj for all elements are computed, which is derived from the integral equivalent of differential equation, also called weak formulation (Huyakorn and Pinder, 1983).
Summary:
This sub-module provided:
1. A systematic approach to understand the flow and transport problems
2. Introduced several processes that influence the coupled (flow and transport) problems..
3. The derivation of governing equation to solve for the flow and transport problem.
4. Approaches to solve the flow and transport problems.
| you can view video on Groundwater Hydrology IV (Coupled Flow and Transport) | |
References:
- Anderson, M. P., and W. W. Woessner (1992), Applied groundwater modeling: simulation of flow and advective transport, Academic Press, San Diego.
- Chahar, B. R. (2014), Groundwater Hydrology, McGraw hill education, New Delhi, India
- Domenico, P. A., and F. W. Schwartz (1998), Physical and chemical hydrogeology, Wiley, New York.
- Fick, A. (1855), Ueber Diffusion, Ann. Phys. Chem., 170(1), 59–86
- Freeze, R. A., and J. A. Cherry (1979), Groundwater, Prentice-Hall, Englewood Cliffs, N.J.
- Huyakorn, P. S., and H. R. Pinder (1983), Computational methods in subsurface flow, Elsevier Science, Oxford.
- Ogata, A., and R. B. Banks (1961), A solution of the differential equation of longitudinal dispersion in porous media, U.S. Geological Survey.