<< . .

. 2
( : 59)

. . >>

0.4 we focus on reaction and transport processes in porous media. These
sections are independent of the remaining parts and may be skipped by
the reader. Section 0.5, however, should be consulted because it ¬xes some
notation to be used later on.

0.1 The Basic Partial Di¬erential Equation Models
Partial di¬erential equations are equations involving some partial deriva-
tives of an unknown function u in several independent variables. Partial
di¬erential equations which arise from the modelling of spatial (and tempo-
ral) processes in nature or technology are of particular interest. Therefore,
we assume that the variables of u are x = (x1 , . . . , xd )T ∈ Rd for d ≥ 1,
representing a spatial point, and possibly t ∈ R, representing time. Thus
the minimal set of variables is (x1 , x2 ) or (x1 , t), otherwise we have ordinary
di¬erential equations. We will assume that x ∈ „¦, where „¦ is a bounded
domain, e.g., a metal workpiece, or a groundwater aquifer, and t ∈ (0, T ] for
some (time horizon) T > 0. Nevertheless also processes acting in the whole
Rd — R, or in unbounded subsets of it, are of interest. One may consult the
Appendix for notations from analysis etc. used here. Often the function u
2 0. Modelling Processes in Porous Media with Di¬erential Equations

represents, or is related to, the volume density of an extensive quantity like
mass, energy, or momentum, which is conserved. In their original form all
quantities have dimensions that we denote in accordance with the Inter-
national System of Units (SI) and write in square brackets [ ]. Let a be
a symbol for the unit of the extensive quantity, then its volume density
is assumed to have the form S = S(u), i.e., the unit of S(u) is a/m3 . For
example, for mass conservation a = kg, and S(u) is a concentration. For
describing the conservation we consider an arbitrary “not too bad” sub-
set „¦ ‚ „¦, the control volume. The time variation of the total extensive
quantity in „¦ is then

‚t S(u(x, t))dx . (0.1)

If this function does not vanish, only two reasons are possible due to con-
” There is an internally distributed source density Q = Q(x, t, u) [a/m3 /s],
being positive if S(u) is produced, and negative if it is destroyed, i.e., one
term to balance (0.1) is „¦ Q(x, t, u(x, t))dx.
” There is a net ¬‚ux of the extensive quantity over the boundary ‚ „¦ of
„¦. Let J = J (x, t) [a/m2 /s] denote the ¬‚ux density, i.e., J i is the amount,
that passes a unit square perpendicular to the ith axis in one second in
the direction of the ith axis (if positive), and in the opposite direction
otherwise. Then another term to balance (0.1) is given by

’ J (x, t) · ν(x)dσ ,

where ν denotes the outer unit normal on ‚„¦. Summarizing the conserva-
tion reads

S(u(x, t))dx = ’ J (x, t) · ν(x)dσ +
‚t Q(x, t, u(x, t))dx . (0.2)
˜ ˜ ˜
„¦ ‚„¦ „¦

The integral theorem of Gauss (see (2.3)) and an exchange of time
derivative and integral leads to

[‚t S(u(x, t)) + ∇ · J (x, t) ’ Q(x, t, u(x, t))]dx = 0 ,

and, as „¦ is arbitrary, also to
‚t S(u(x, t)) + ∇ · J(x, t) = Q(x, t, u(x, t)) for x ∈ „¦, t ∈ (0, T ] . (0.3)
All manipulations here are formal assuming that the functions involved
have the necessary properties. The partial di¬erential equation (0.3) is the
basic pointwise conservation equation, (0.2) its corresponding integral form.
Equation (0.3) is one requirement for the two unknowns u and J , thus it
0.1. The Basic Partial Di¬erential Equation Models 3

has to be closed by a (phenomenological) constitutive law, postulating a
relation between J and u.
Assume „¦ is a container ¬lled with a ¬‚uid in which a substance is dis-
solved. If u is the concentration of this substance, then S(u) = u and a
= kg. The description of J depends on the processes involved. If the ¬‚uid
is at rest, then ¬‚ux is only possible due to molecular di¬usion, i.e., a ¬‚ux
from high to low concentrations due to random motion of the dissolved
particles. Experimental evidence leads to
J (1) = ’K∇u (0.4)
with a parameter K > 0 [m2 /s], the molecular di¬usivity. Equation (0.4)
is called Fick™s law.
In other situations, like heat conduction in a solid, a similar model occurs.
Here, u represents the temperature, and the underlying principle is energy
conservation. The constitutive law is Fourier™s law, which also has the form
(0.4), but as K is a material parameter, it may vary with space or, for
anisotropic materials, be a matrix instead of a scalar.
Thus we obtain the di¬usion equation
‚t u ’ ∇ · (K∇u) = Q . (0.5)
If K is scalar and constant ” let K = 1 by scaling ”, and f := Q is
independent of u, the equation simpli¬es further to
‚t u ’ ∆u = f ,
where ∆u := ∇·(∇u) . We mentioned already that this equation also occurs
in the modelling of heat conduction, therefore this equation or (0.5) is also
called the heat equation.
If the ¬‚uid is in motion with a (given) velocity c then (forced) convection
of the particles takes place, being described by
J (2) = uc , (0.6)
i.e., taking both processes into account, the model takes the form of the
convection-di¬usion equation
‚t u ’ ∇ · (K∇u ’ cu) = Q . (0.7)
The relative strength of the two processes is measured by the P´clete
number (de¬ned in Section 0.4). If convection is dominating one may ignore
di¬usion and only consider the transport equation
‚t u + ∇ · (cu) = Q . (0.8)
The di¬erent nature of the two processes has to be re¬‚ected in the models,
therefore, adapted discretization techniques will be necessary. In this book
we will consider models like (0.7), usually with a signi¬cant contribution
of di¬usion, and the case of dominating convection is studied in Chapter
9. The pure convective case like (0.8) will not be treated.
4 0. Modelling Processes in Porous Media with Di¬erential Equations

In more general versions of (0.7) ‚t u is replaced by ‚t S(u), where S
depends linearly or nonlinearly on u. In the case of heat conduction S is
the internal energy density, which is related to the temperature u via the
factors mass density and speci¬c heat. For some materials the speci¬c heat
depends on the temperature, then S is a nonlinear function of u.
Further aspects come into play by the source term Q if it depends linearly
or nonlinearly on u, in particular due to (chemical) reactions. Examples for
these cases will be developed in the following sections. Since equation (0.3)
and its examples describe conservation in general, it still has to be adapted
to a concrete situation to ensure a unique solution u. This is done by the
speci¬cation of an initial condition
for x ∈ „¦ ,
S(u(x, 0)) = S0 (x)
and by boundary conditions. In the example of the water ¬lled container
no mass ¬‚ux will occur across its walls, therefore, the following boundary
J · ν(x, t) = 0 for x ∈ ‚„¦, t ∈ (0, T ) (0.9)
is appropriate, which ” depending on the de¬nition of J ” prescribes the
normal derivative of u, or a linear combination of it and u. In Section 0.5
additional situations are depicted.
If a process is stationary, i.e. time-independent, then equation (0.3)
reduces to
∇ · J (x) = Q(x, u(x)) for x ∈ „¦ ,
which in the case of di¬usion and convection is speci¬ed to
’∇ · (K∇u ’ cu) = Q .
For constant K ” let K = 1 by scaling ”, c = 0, and f := Q, being
independent of u, this equation reduces to
’∆u = f in „¦ ,
the Poisson equation.
Instead of the boundary condition (0.9), one can prescribe the values of
the function u at the boundary:
for x ∈ ‚„¦ .
u(x) = g(x)
For models , where u is a concentration or temperature, the physical reali-
sation of such a boundary condition may raise questions, but in mechanical
models, where u is to interpreted as a displacement, such a boundary con-
dition seems reasonable. The last boundary value problem will be the ¬rst
model, whose discretization will be discussed in Chapters 1 and 2.
Finally it should be noted that it is advisable to non-dimensionalise the
¬nal model before numerical methods are applied. This means that both
the independent variables xi (and t), and the dependent one u, are replaced
0.2. Reactions and Transport in Porous Media 5

by xi /xi,ref , t/tref , and u/uref , where xi,ref , tref , and uref are ¬xed reference
values of the same dimension as xi , t, and u, respectively. These reference
values are considered to be of typical size for the problems under investiga-
tion. This procedure has two advantages: On the one hand, the typical size
is now 1, such that there is an absolute scale for (an error in) a quantity
to be small or large. On the other hand, if the reference values are chosen
appropriately a reduction in the number of equation parameters like K
and c in (0.7) might be possible, having only fewer algebraic expressions of
the original material parameters in the equation. This facilitates numerical
parameter studies.

0.2 Reactions and Transport in Porous Media
A porous medium is a heterogeneous material consisting of a solid matrix
and a pore space contained therein. We consider the pore space (of the
porous medium) as connected; otherwise, the transport of ¬‚uids in the
pore space would not be possible. Porous media occur in nature and man-
ufactured materials. Soils and aquifers are examples in geosciences; porous
catalysts, chromatographic columns, and ceramic foams play important
roles in chemical engineering. Even the human skin can be considered a
porous medium. In the following we focus on applications in the geosciences.
Thus we use a terminology referring to the natural soil as a porous medium.
On the micro or pore scale of a single grain or pore, i.e., in a range of µm
to mm, the ¬‚uids constitute di¬erent phases in the thermodynamic sense.
Thus we name this system in the case of k ¬‚uids including the solid matrix
as (k + 1)-phase system or we speak of k-phase ¬‚ow.
We distinguish three classes of ¬‚uids with di¬erent a¬nities to the solid
matrix. These are an aqueous phase, marked with the index “w” for water,
a nonaqueous phase liquid (like oil or gasoline as natural resources or con-
taminants), marked with the index “o,” and a gaseous phase, marked with
the index “g” (e.g., soil air). Locally, at least one of these phases has al-
ways to be present; during a transient process phases can locally disappear
or be generated. These ¬‚uid phases are in turn mixtures of several com-
ponents. In applications of the earth sciences, for example, we do not deal
with pure water but encounter di¬erent species in true or colloidal solu-
tion in the solvent water. The wide range of chemical components includes
plant nutrients, mineral nutrients from salt domes, organic decomposition
products, and various organic and inorganic chemicals. These substances
are normally not inert, but are subject to reactions and transformation
processes. Along with di¬usion, forced convection induced by the motion
of the ¬‚uid is the essential driving mechanism for the transport of solutes.
But we also encounter natural convection by the coupling of the dynamics
of the substance to the ¬‚uid ¬‚ow. The description level at the microscale
6 0. Modelling Processes in Porous Media with Di¬erential Equations

that we have used so far is not suitable for processes at the laboratory or
technical scale, which take place in ranges of cm to m, or even for processes
in a catchment area with units of km. For those macroscales new models
have to be developed, which emerge from averaging procedures of the mod-
els on the microscale. There may also exist principal di¬erences among the
various macroscales that let us expect di¬erent models, which arise from
each other by upscaling. But this aspect will not be investigated here fur-
ther. For the transition of micro to macro scales the engineering sciences
provide the heuristic method of volume averaging, and mathematics the
rigorous (but of only limited use) approach of homogenization (see [36] or
[19]). None of the two possibilities can be depicted here completely. Where
necessary we will refer to volume averaging for (heuristic) motivation.
Let „¦ ‚ Rd be the domain of interest. All subsequent considerations are
formal in the sense that the admissibility of the analytic manipulations is
supposed. This can be achieved by the assumption of su¬cient smoothness
for the corresponding functions and domains.
Let V ‚ „¦ be an admissible representative elementary volume in the
sense of volume averaging around a point x ∈ „¦. Typically the shape and
the size of a representative elementary volume are selected in such a manner
that the averaged values of all geometric characteristics of the microstruc-
ture of the pore space are independent of the size of V but depend on
the location of the point x. Then we obtain for a given variable ω± in the
phase ± (after continuation of ω± with 0 outside of ±) the corresponding
macroscopic quantities, assigned to the location x, as the extrinsic phase

ω± := ω±
|V | V

or as the intrinsic phase average

ω± := ω± .
|V± | V±

Here V± denotes the subset of V corresponding to ±. Let t ∈ (0, T ) be
the time at which the process is observed. The notation x ∈ „¦ means the
vector in Cartesian coordinates, whose coordinates are referred to by x,
y, and z ∈ R. Despite this ambiguity the meaning can always be clearly
derived from the context.
Let the index “s” (for solid) stand for the solid phase; then

φ(x) := |V \ Vs | |V | > 0

denotes the porosity, and for every liquid phase ±,

S± (x, t) := |V± | |V \ Vs | ≥ 0
0.3. Fluid Flow in Porous Media 7

is the saturation of the phase ±. Here we suppose that the solid phase is
stable and immobile. Thus
ω± = φS± ω±
for a ¬‚uid phase ± and

S± = 1 . (0.10)

So if the ¬‚uid phases are immiscible on the micro scale, they may be miscible
on the macro scale, and the immiscibility on the macro scale is an additional
assumption for the model.
As in other disciplines the di¬erential equation models are derived here
from conservation laws for the extensive quantities mass, impulse, and en-
ergy, supplemented by constitutive relationships, where we want to focus
on the mass.

0.3 Fluid Flow in Porous Media
Consider a liquid phase ± on the micro scale. In this chapter, for clarity, we
write “short” vectors in Rd also in bold with the exception of the coordinate
vector x. Let ˜± [kg/m3 ] be the (microscopic) density, q ± :=
˜ ˜
· ˜· v · ˜±
[m/s] the mass average mixture velocity based on the particle velocity v · of
a component · and its concentration in solution ˜· [kg/m ]. The transport
theorem of Reynolds (see, for example, [10]) leads to the mass conservation
‚t ˜± + ∇ · (˜± q ± ) = f±
˜ (0.11)
with a distributed mass source density f± . By averaging we obtain from
here the mass conservation law
‚t (φS± ±) ±q±) = f± (0.12)
with ± , the density of phase ±, as the intrinsic phase average of ˜± and
q ± , the volumetric ¬‚uid velocity or Darcy velocity of the phase ±, as the
extrinsic phase average of q ± . Correspondingly, f± is an average mass source
Before we proceed in the general discussion, we want to consider some
speci¬c situations: The area between the groundwater table and the imper-
meable body of an aquifer is characterized by the fact that the whole pore
space is occupied by a ¬‚uid phase, the soil water. The corresponding satu-
ration thus equals 1 everywhere, and with omission of the index equation
(0.12) takes the form
‚t (φ ) + ∇ · ( q) = f . (0.13)
8 0. Modelling Processes in Porous Media with Di¬erential Equations

If the density of water is assumed to be constant, due to neglecting
the mass of solutes and compressibility of water, equation (0.13) simpli¬es
further to the stationary equation
∇·q =f , (0.14)
where f has been replaced by the volume source density f / , keeping the
same notation. This equation will be completed by a relationship that
can be interpreted as the macroscopic analogue of the conservation of mo-
mentum, but should be accounted here only as an experimentally derived
constitutive relationship. This relationship is called Darcy™s law, which
reads as
q = ’K (∇p + gez ) (0.15)
and can be applied in the range of laminar ¬‚ow. Here p [N/m2 ] is the intrinsic
average of the water pressure, g [m/s2 ] the gravitational acceleration, ez the
unit vector in the z-direction oriented against the gravitation,
K = k/µ , (0.16)
a quantity, which is given by the permeability k determined by the solid
phase, and the viscosity µ determined by the ¬‚uid phase. For an anisotropic
solid, the matrix k = k(x) is a symmetric positive de¬nite matrix.
Inserting (0.15) in (0.14) and replacing K by K g, known as hydraulic
conductivity in the literature, and keeping the same notation gives the
following linear equation for
h(x, t) := p(x, t) + z ,
the piezometric head h [m]:
’∇ · (K∇h) = f . (0.17)
The resulting equation is stationary and linear. We call a di¬erential equa-
tion model stationary if it depends only on the location x and not on the
time t, and instationary otherwise. A di¬erential equation and correspond-
ing boundary conditions (cf. Section 0.5) are called linear if the sum or a
scalar multiple of a solution again forms a solution for the sum, respectively
the scalar multiple, of the sources.
If we deal with an isotropic solid matrix, we have K = KI with the d— d
unit matrix I and a scalar function K. Equation (0.17) in this case reads
’∇ · (K∇h) = f . (0.18)
Finally if the solid matrix is homogeneous, i.e., K is constant, we get from
division by K and maintaining the notation f the Poisson equation
’∆h = f , (0.19)
0.3. Fluid Flow in Porous Media 9

which is termed the Laplace equation for f = 0. This model and its more
general formulations occur in various contexts. If, contrary to the above as-
sumption, the solid matrix is compressible under the pressure of the water,
and if we suppose (0.13) to be valid, then we can establish a relationship
φ = φ(x, t) = φ0 (x)φf (p)
with φ0 (x) > 0 and a monotone increasing φf such that with S(p) := φf (p)
we get the equation
φ0 S(p) ‚t p + ∇ · q = f
and the instationary equations corresponding to (0.17)“(0.19), respectively.
For constant S(p) > 0 this yields the following linear equation:
φ0 S ‚t h ’ ∇ · (K∇h) = f , (0.20)
which also represents a common model in many contexts and is known from
corresponding ¬elds of application as the heat conduction equation.
We consider single phase ¬‚ow further, but now we will consider gas as
¬‚uid phase. Because of the compressibility, the density is a function of the
pressure, which is invertible due to its strict monotonicity to
p = P( ) .
Together with (0.13) and (0.15) we get a nonlinear variant of the heat
conduction equation in the unknown :
‚t (φ ) ’ ∇ · K( ∇P ( ) + 2
gez ) = f , (0.21)
which also contains derivatives of ¬rst order in space. If P ( ) = ln(± ) holds
for a constant ± > 0, then ∇P ( ) simpli¬es to ±∇ . Thus for horizontal
¬‚ow we again encounter the heat conduction equation. For the relationship
P ( ) = ± suggested by the universal gas law, ± ∇ = 1 ±∇ 2 remains
nonlinear. The choice of the variable u := 2 would result in u1/2 in the
time derivative as the only nonlinearity. Thus in the formulation in the
coe¬cient of ∇ disappears in the divergence of = 0. Correspondingly,
the coe¬cient S(u) = 1 φu’1/2 of ‚t u in the formulation in u becomes
unbounded for u = 0. In both versions the equations are degenerate, whose
treatment is beyond the scope of this book. A variant of this equation has
gained much attention as the porous medium equation (with convection) in
the ¬eld of analysis (see, for example, [42]).
Returning to the general framework, the following generalization of
Darcy™s law can be justi¬ed experimentally for several liquid phases:
q± = ’ k (∇p± + ± gez ) .
Here the relative permeability kr± of the phase ± depends upon the
saturations of the present phases and takes values in [0, 1].
10 0. Modelling Processes in Porous Media with Di¬erential Equations

At the interface of two liquid phases ±1 and ±2 we observe a di¬erence of
the pressures, the so-called capillary pressure, that turns out experimentally
to be a function of the saturations:

pc ±1 ±2 := p±1 ’ p±2 = F±1 ±2 (Sw , So , Sg ) . (0.22)

A general model for multiphase ¬‚ow, formulated for the moment in terms
of the variables p± , S± , is thus given by the equations

‚t (φS± ±) ± »± k(∇p± + ± gez )) = f± (0.23)

<< . .

. 2
( : 59)

. . >>