Updated: Jan 10, 2001.
Copyright © 2000 by Jacob Bear, Haifa, Israel. 
All Rights Reserved.

This is an incomplete version for demo purposes only

Computer-Mediated Distance Learning

Course on

MODELING GROUNDWATER FLOW AND CONTAMINANT TRANSPORT

INSTRUCTOR:

JACOB BEAR

Professor Emeritus

Faculty of Civil Engineering

Technion-Israel Institute of Technology

Haifa 32000, Israel

TOPIC D:

MODELING FLOW IN AQUIFERS

TOPIC B was devoted to the MOTION EQUATION. We have presented this equation for the two basic flow models:

Flow in a three-dimensional domain.

Flow in a domain whose geometry justifies the "essentially horizontal flow approximation", i.e., the flow domain is approximated as a horizontal two-dimensional one.

In both flow "modes", we have considered extensions of the original (experimentally derived) Darcy's law, viz., extensions to

Three-dimensional spatial domains.

Inhomogeneous domains.

Anisotropic domains.

Compressible fluids.

Inhomogeneous fluids (i.e., concentration, and/or temperature dependent density).

Deformable, and\or non-stationary solid matrix.

In all cases, we have shown how the specific discharge is related to the piezometric head, or to the fluid's pressure.

The first natural question, especially if at some earlier studies of basic groundwater flow you have used Darcy's law to solve some simple flow problems, is:

IS DARCY'S LAW SUFFICIENT FOR SOLVING GROUND WATER FLOW PROBLEMS, OR IS THERE A NEED FOR ADDITIONAL INFORMATION, OR EQUATION(S)?

Typically, in three-dimensional flow, the specific discharge vector is expressed by:

q = - K Ñh, (1)

in which we have TWO state variables:

The specific discharge q = q(x, y, z, t) and

The piezometric head h = h(x, y, z, t).

Equation (1) is written as a vector equation (instead of writing three scalar ones). The dot indicates a dot product. Just to make sure that you understand this nomenclature,

WRITE Eqn (1) as scalar equations. For the sake of simplicity, assume that the soil is anisotropic, but that x, y, z are principal directions.

ANSWER

Recall that the main objective of modeling flow in a given porous medium, or aquifer domain, is to predict the spatial and temporal distribution of the piezometric head, h = h (x, y, z, t), within the domain, in response to excitation of the domain (= the system).

Later we shall discuss in detail the meaning of "excitation". Here, think of "excitation" as what we do to the system (e.g., pump from it or inject into it) in order to achieve goals.

To predict the distribution of h (x, y, z, t), we have to solve a mathematical model, i.e., a set of equations for the TWO unknown variables.

The specific discharge vector q = q(x, y, z, t) and

The piezometric head h = h(x, y, z, t).

However, the only available vector equation seems to be Darcy's law,

q = - K Ñh, (1)

given above. How can we solve a single equation for two variables?

Note that we could say that we have FOUR unknown scalar variables: q_x, q_y, q_z, and h. For these variables, we have three scalar equations:

q_x = - K , q_y = - K , q_z = - K . (2)

We have written these equations for the case of an isotropic porous medium domain. However, you should have no difficulty in writing the three scalar equations for an anisotropic one. TRY. If you do not remember how to express Darcy's law in a scalar form for an anisotropic porous medium, Click here.

Again, UNFORTUNATELY, we have only THREE SCALAR EQUATIONS . We know that it is

Accordingly, we have to discuss the requirement that the balance of fluid mass be maintained within a considered porous medium domain. We shall express this requirement in the form of a 'mass balance equation'. Later we shall see that the description of transport of a chemical species in a fluid phase requires that a mass balance of that species be also satisfied.

Our objective is to develop the mass balance equation. We intend to do it for one, two, and three fluid phases in a porous medium domain. Deformation of the solid matrix will also be taken into account. For a given fluid phase, by combining the mass balance equation with an appropriate motion equation, a flow equation for that phase is obtained. The complete flow model requires, in addition, specification of appropriate initial and boundary conditions. Also, the numerical values of all the coefficients which appear in the model equations must be specified. All this supplemental information must be specified for the considered domain and fluids, in order to obtain solutions that pertain to the particular flow problem under consideration.

LECTURE 1: Fundamental Balance Equation

Following the above introductory remarks, I suggest that we proceed in the following steps:

Although our direct objective is to develop a balance equation, for the mass of a fluid phase, since the same development and principles apply also to the balance of any extensive quantity (e.g., energy, or mass of a component of a fluid phase), we shall start by considering the most general balance equation for ANY extensive quantity. This will allow us later to apply this general equation to a number of extensive quantities of interest.

In TOPIC B, LECTURE 3, we have discussed a regional water balance, i.e., a balance performed over a finite, well-defined aquifer domain. We could not use the balance to predict the behavior (piezometric head, specific discharge, etc.) at points within the considered domain, or as they change with time. We have assumed a constant density and hence, we have actually written a volume balance.

To provide such information, we shall express the balance equation in the form of a partial differential equation. Think of a partial differential equation as a "statement about what happens at a specified point within a considered domain at a specified instant of time".

Let me remind you:

At the microscopic level, we consider what happens at a point within a phase (including the case when this phase, say a fluid phase, occupies part of the void space of a porous medium domain).

At the macroscopic level, we consider what happens at a (macroscopic) point within a porous medium domain. The behavior of the phase at such point is obtained by averaging the behavior of the considered phase within a representative elementary volume (REV), centered at the point, and assigning the average value to the considered point.

At the macroscopic level, the porous medium domain is considered as a continuum.

Remember that we are interested in developing the mass balance equation, or the balance equation for any extensive quantities, at the macroscopic level. In principle, such equation can be obtained in two ways:

Start from a mass balance equation written at the microscopic level. Then we average this equation, to obtain the macroscopic equation.

(In parenthesis, we need to discuss the meaning of "averaging" as a tool for passing from the microscopic level to the macroscopic one. However, at this stage, just think of the intuitive meaning of the word. We may come back to this topic at a later stage.). Or,

Develop the mass balance equation at the macroscopic level, already considered as a continuum.

In fact, the same approach, of developing the balance equation directly, is valid also when developing the balance equation at the microscopic level, also considered as a continuum.

Accordingly, let us start by developing the balance equation at the microscopic level. This means that we are considering what happens at every point within a fluid phase continuum that occupies a specified spatial domain.

However, as suggested above, before discussing the specific subject of "mass balance equation", we shall, first, discuss the more general case of a balance equation for any extensive quantity. It will then be a straightforward procedure to apply the results to the particular case in which the mass of a phase is the considered extensive quantity.

The basic idea of a balance is very simple and straightforward:

A balance can be written for any extensive quantity, E, within a specified spatial domain, and for a specified period of time.

Let us denote the time interval by Dt, the volume of the porous medium domain for which the balance is written by U, and the surface bounding it, by S. We shall denote the density of E, i.e., the quantity of E per unit volume of the considered phase, by e.

The balance is then stated as

Does this make sense?

Let us consider these BALANCE COMPONENTS in detail:

The components of J^tE are: j^tE|_x , j^tE|_y , and j^tE|_z . For the sake of simplicity, we shall use J^tE and J, i.e., without the superscripts, interchangeably.

DIVERGENCE OF A FLUX

Let us introduce at this time the mathematical term 'divergence of a flux', of any extensive quantity, and its physical interpretation. It will play an essential role in any balance equation that takes the form of a partial differential equation. It is very important that you grasp the physical meaning of this term.

Consider the parallelepiped control volume with dimensions Dx, Dy, Dz, shown in the figure below.

Figure 1: A control volume

During a period Dt, the excess of inflow over outflow of the considered extensive quantity, E, through the surfaces that bound the control volume, is expressed (using the symbol j for total flux of E) by

The excess of inflow over outflow of the considered extensive quantity, E, per unit time, through the surfaces that bound the control volume, is obtained by dividing by Dt:

Next, we divide by DxDyDz, in order to obtain the excess of inflow over outflow of the considered extensive quantity, E, per unit time, through the surfaces that bound the control volume, per unit volume:

Next we wish to learn "what is this excess at a point. To achieve this goal, we go to the limit of the above expression , as Dx,Dy,Dz ֶ 0. We obtain:

Recalling the definition of a derivative:

We may now write:

In mathematics, given a vector, j, the negative of the sum of the three partial derivatives: the derivative with respect to x of the x-component of vector, j, plus the derivative with respect to y of the y-component of j, plus the derivative with respect to z of the z-component of the vector j, is a scalar called the divergence of the vector. Symbolically, this is written in the form:

The two expressions are equivalent!

In words, we say "divergence of the vector j ".

In defining the divergence of any vector, J, we have used here Cartesian coordinates, x, y, z. However, as emphasized earlier in recommending the use of vector notation, once we have the term "div J" in an equation, we can rewrite it in any other convenient coordinates. For example:

In polar coordinates: r, q.

To summarize:

The excess of inflow over outflow of an extensive quantity, per unit volume and per unit time, is expressed by div (flux vector), (or Ñ(flux vector)). This is a consequence of the spatial distribution of the flux. In this way, we have presented the physical interpretation of the mathematical symbol, or expressed the physical statement by a mathematical symbol. The latter is more compact, but the content is the same.

Using the above expression, and recalling the general structure of the partial differential form of the balance equation, we can immediately write the balance of any extensive quantity without the need to develop it from first principles. Let me suggest the (strange?!) form

This does not look like a "normal" equation. Doesn't it? NEVERTHELESS, ...

Whenever you need to write a balance equation, you have to determine:

The quantity of E per unit volume, (at the microscopic level considered here, this is e),
The total flux of E, and
The rate of production of E, per unit volume,

And using the above "verbal equation", you'll be able to write the differential form of the balance equation IMMEDIATELY.

Note, and never forget to put, the dot after the nabla symbol. This will mean that the "nabla dot" symbol indicates divergence of a vector, while without the dot, we have the gradient of a scalar (e.g., Ñh).

DIFFUSIVE FLUX

We still need to discuss the meaning of "velocity V^E of an extensive quantity, E."

For the case E = m, i.e., the extensive quantity is the mass of a phase,

V^E º V^m º V . (13)

This is rather simple. Thus, V^m, also denoted by the symbol V, is the velocity of the mass continuum of a phase.

However, consider the example:
E = m^g, e = r^g., i.e., the extensive quantity is the mass of a g-component in a phase. What is then the meaning of V^E= V^g?

This g-component is composed of many molecules which move continuously within the phase (random, Brownian motion). As discussed in TOPIC A, LECTURE 2, at the molecular level, we "see" the individual molecules running around. However, we are not interested in the detailed behavior at such level. Instead, we average the behavior of the many individual molecules that run around, within small representative elementary volumes, called micro-REV (mREV), and assign the result to the centroid of each such mREV. In this way, the individual molecules "disappear", and we obtain at every point only the average behavior. For example, we get the average density, r^g, and the average velocity, V^g. In other words, we replace the domain with molecules by a domain with a FLUID that has properties derived by averaging properties of the molecules.

Altogether, at every point within a phase, now considered as a continuum, we have two velocities:

V of the mass of the phase as a whole, regarded as a continuum, and

V^gof the mass of the g-component also regarded as a continuum.

The two continua occupy the same spatial domain. At every point within the latter, we have properties of both continua.

Also, we have to consider the fluxes:

The mass flux, rV, and

The flux of the mass of the g-component, r^gV^g.

However, how do we determine the velocity V^g ?

We may now return to the fundamental balance equation and rewrite it, making use of the above decomposition:

= -Ñ (eV + J^E) + rG ^E, J^E = e(V^E - V). (16)

This is the most general differential balance equation for any extensive quantity, E. Again, it is a balance at a point, in the sense that it describes what happens on the average within a small volume surrounding the point, with all values assigned to the point.

We have, thus, achieved the goal of stating the general differential balance equation for any extensive quantity, E. In the next lecture, we'll return to the balance of mass of a fluid phase.

Before turning to the next lecture, let's rewrite the general balance equation, written above in a vector form, in terms of the more commonly used coordinate systems

A. Cartesian coordinates:

At a later stage, we'll discuss how to express the velocity and the diffusive flux.

B. Cylindrical coordinates:

DIRECT APPROACH

Based on the discussion at the microscopic level, the E balance equation at the macroscopic level, can be stated in the form:

In this equation:

The quantity of E contained in a unit volume of porous medium is expressed by eq :

(E per unit volume of porous medium) =
(E per unit volume of fluid) ¥ (volume of fluid per unit volume of porous medium).

The rate of increase in this quantity is expressed by . This term appears on the l.h.s.

The total E flux at the macroscopic level is expressed by , where the symbol represents the intrinsic phase

However, if there exists also transfer of E from the considered a-fluid phase to any adjacent (fluid or solid) b-phase, through their common (microscopic interface), and we wish to express this transfer per unit volume of porous medium, we add a source denoted by the symbol .

Altogether, we may now write the macroscopic mass balance equation in the form:

I hope you can now interpret each term in the above equation as a symbolic expression of the rate at which the quantity of E is added per unit volume of porous medium.

Now to the second approach:

BY AVERAGING

Let's start by writing the microscopic balance equation for an extensive quantity E of a phase:

= -Ñ (eV + J^E) + rG ^E. (20)

Again, each term expresses the rate at which the quantity of E is added per unit volume, but this time, it is per unit volume of the fluid.

Let's consider the case in which the considered fluid phase is occupying only part of the void space, at volumetric fraction, q.

vanishes. If it is material with respect to mass, e(V - u) vanishes, and the flux through the S_ab-surface is due to diffusion only.

What is a "material surface"?

A "material surface with respect to an E-quantity" is surface such that no E-quantity crosses it. One may also have a

The figure below shows an interface, its velocity, the unit normal vector and the advective flux.

Denoting this interface transfer term by , we may rewrite Eqn. (21) in the form

Physical interpretation: The net added quantity, represented by the term on the l.h.s, is due to the r.h.s. terms:

A net influx of E, per unit time and per unit volume of porous medium, due to spatial variations in the total E -flux.

A quantity of E that leaves the phase through the S_ab-surface that bounds the phase within the REV (of volume U₀), per unit time and per unit volume of porous medium.

deviations of the density and the velocity (at the microscopic level) from their respective averages.

A diffusive flux, obtained as the average of the microscopic one (and, again, multiplied by q_{a
.}

Note that, obviously, we have arrived at the same equation by the two approaches.

We have now the general differential form of the macroscopic balance equation for any E

within an a-phase, in a (possibly) multiphase system.

Question

To check if you have understood the above discussion, consider the case in which we consider the transport of the mass of a g-component in an a-fluid phase (of density r). The g-component is produced at the rate G^g. Let c denote the mass concentration of the g-component in the a-phase. Use the symbol to denote the rate of transfer of the g-component from the a-fluid phase to all other b-(fluid and solid) phases present in the system

Write the equation and compare with my answer.

You should have no difficulty in using this fundamental balance equation for any extensive quantity of interest. In the next lecture, we'll apply this balance equation to the mass of a fluid phase.

The lecture, TOPIC D, LECTURE 2
deals with the mass balance equation for a phase.