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-

servation:

” 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

condition

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

average

1

ω± := ω±

|V | V

or as the intrinsic phase average

1

±

ω± := ω± .

|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)

±:¬‚uid

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

3

a component · and its concentration in solution ˜· [kg/m ]. The transport

theorem of Reynolds (see, for example, [10]) leads to the mass conservation

law

˜

‚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

density.

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

1

h(x, t) := p(x, t) + z ,

g

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

2

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

2

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:

kr±

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)