<< . .

. 3
( : 59)

. . >>

with the mobilities »± := kr± /µ± , and the equations (0.22) and (0.10),
where one of the S± ™s can be eliminated. For two liquid phases w and g,
e.g., water and air, equations (0.22) and (0.10) for ± = w, g read pc =
pg ’ pw = F (Sw ) and Sg = 1 ’ Sw . Apparently, this is a time-dependent,
nonlinear model in the variables pw , pg , Sw , where one of the variables can
be eliminated. Assuming constant densities ± , further formulations based

∇ · q w + q g = fw / + fg / (0.24)
w g

can be given as consequences of (0.10). These equations consist of a sta-
tionary equation for a new quantity, the global pressure, based on (0.24),
and a time-dependent equation for one of the saturations (see Exercise 0.2).
In many situations it is justi¬ed to assume a gaseous phase with constant
pressure in the whole domain and to scale this pressure to pg = 0. Thus
for ψ := pw = ’pc we have

φ‚t S(ψ) ’ ∇ · (»(ψ)k(∇ψ + gez )) = fw / (0.25)

:= w , and S(ψ) := F ’1 (’ψ) as a strictly
with constant pressure
monotone increasing nonlinearity as well as ».
With the convention to set the value of the air pressure to 0, the pressure
in the aqueous phase is in the unsaturated state, where the gaseous phase is
also present, and represented by negative values. The water pressure ψ = 0
marks the transition from the unsaturated to the saturated zone. Thus
in the unsaturated zone, equation (0.25) represents a nonlinear variant
of the heat conduction equation for ψ < 0, the Richards equation. As
most functional relationships have the property S (0) = 0, the equation
degenerates in the absence of a gaseous phase, namely to a stationary
equation in a way that is di¬erent from above.
Equation (0.25) with S(ψ) := 1 and »(ψ) := »(0) can be continued in a
consistent way with (0.14) and (0.15) also for ψ ≥ 0, i.e., for the case of a
sole aqueous phase. The resulting equation is also called Richards equation
or a model of saturated-unsaturated ¬‚ow.
0.4. Reactive Solute Transport in Porous Media 11

0.4 Reactive Solute Transport in Porous Media
In this chapter we will discuss the transport of a single component in a
liquid phase and some selected reactions. We will always refer to water
as liquid phase explicitly. Although we treat inhomogeneous reactions in
terms of surface reactions with the solid phase, we want to ignore exchange
processes between the ¬‚uid phases. On the microscopic scale the mass con-
servation law for a single component · is, in the notation of (0.11) by
omitting the phase index w,
‚t ˜· + ∇ · (˜· q ) + ∇ · J · = Q· ,
J · := ˜· (˜ · ’ q ) [kg/m2 /s]
v (0.26)
represents the di¬usive mass ¬‚ux of the component · and Q· [kg/m3 /s] is
its volumetric production rate. For a description of reactions via the mass
action law it is appropriate to choose the mole as the unit of mass. The
di¬usive mass ¬‚ux requires a phenomenological description. The assump-
tion that solely binary molecular di¬usion, described by Fick™s law, acts
between the component · and the solvent, means that
J · = ’ ˜D· ∇ (˜· / ˜) (0.27)
with a molecular di¬usivity D· > 0 [m2 /s]. The averaging procedure applied
on (0.26), (0.27) leads to
‚t (˜c· ) + ∇ · (qc· ) + ∇ · J (1) + ∇ · J (2) = Q(1) + Q(2)
· ·

for the solute concentration of the component ·, c· [kg/m3 ], as intrinsic
phase average of ˜· . Here, we have J (1) as the average of J · and J (2) ,
the mass ¬‚ux due to mechanical dispersion, a newly emerging term at the
(1) ˜
macroscopic scale. Analogously, Q· is the intrinsic phase average of Q· ,
and Q· is a newly emerging term describing the exchange between the
liquid and solid phases.
The volumetric water content is given by ˜ := φSw with the water
saturation Sw . Experimentally, the following phenomenological descriptions
are suggested:
J (1) = ’˜„ D· ∇c·
with a tortuosity factor „ ∈ (0, 1],
J (2) = ’˜Dmech ∇c· , (0.28)
and a symmetric positive de¬nite matrix of mechanical dispersion D mech ,
which depends on q/˜. Consequently, the resulting di¬erential equation
‚t (˜c· ) + ∇ · (qc· ’ ˜D∇c· ) = Q· (0.29)
12 0. Modelling Processes in Porous Media with Di¬erential Equations

(1) (2)
with D := „ D · + D mech, Q· := Q· + Q· .
Because the mass ¬‚ux consists of qc· , a part due to forced convection, and
of J (1) + J (2) , a part that corresponds to a generalized Fick™s law, an equa-
tion like (0.29) is called a convection-di¬usion equation. Accordingly, for
the part with ¬rst spatial derivatives like ∇ · (qc· ) the term convective part
is used, and for the part with second spatial derivatives like ’∇ · (˜D∇c· )
the term di¬usive part is used. If the ¬rst term determines the character of
the solution, the equation is called convection-dominated. The occurrence
of such a situation is measured by the quantity Pe, the global P´clet num-
ber, that has the form Pe = q L/ ˜D [ - ]. Here L is a characteristic
length of the domain „¦. The extreme case of purely convective transport
results in a conservation equation of ¬rst order. Since the common mod-
els for the dispersion matrix lead to a bound for Pe, the reduction to the
purely convective transport is not reasonable. However, we have to take
convection-dominated problems into consideration.
Likewise, we speak of di¬usive parts in (0.17) and (0.20) and of (nonlin-
ear) di¬usive and convective parts in (0.21) and (0.25). Also, the multiphase
transport equation can be formulated as a nonlinear convection-di¬usion
equation by use of (0.24) (see Exercise 0.2), where convection often dom-
inates. If the production rate Q· is independent of c· , equation (0.29) is
In general, in case of a surface reaction of the component ·, the kinetics of
the reaction have to be described . If this component is not in competition
with the other components, one speaks of adsorption. The kinetic equation
thus takes the general form
‚t s· (x, t) = k· f· (x, c· (x, t), s· (x, t)) (0.30)
with a rate parameter k· for the sorbed concentration s· [kg/kg], which is
given in reference to the mass of the solid matrix. Here, the components
in sorbed form are considered spatially immobile. The conservation of the
total mass of the component undergoing sorption gives
Q(2) = ’ b ‚t s· (0.31)

with the bulk density b = s (1’φ), where s denotes the density of the solid
phase. With (0.30), (0.31) we have a system consisting of an instationary
partial and an ordinary di¬erential equation (with x ∈ „¦ as parameter). A
widespread model by Langmuir reads
f· = ka c· (s· ’ s· ) ’ kd s·
with constants ka , kd that depend upon the temperature (among other
factors), and a saturation concentration s· (cf. for example [24]). If we
assume f· = f· (x, c· ) for simplicity, we get a scalar nonlinear equation in
c· ,
‚t (˜c· ) + ∇ · (qc· ’ ˜D∇c· ) + = Q(1) ,
b k· f· (·, c· ) (0.32)
0.4. Reactive Solute Transport in Porous Media 13

and s· is decoupled and extracted from (0.30). If the time scales of transport
and reaction di¬er greatly, and the limit case k· ’ ∞ is reasonable, then
(0.30) is replaced by
f· (x, c· (x, t), s· (x, t)) = 0 .
If this equation is solvable for s· , i.e.,
s· (x, t) = •· (x, c· (x, t)) ,
the following scalar equation for c· with a nonlinearity in the time
derivative emerges:
+ ∇ · (qc· ’ ˜D∇c· ) = Q(1) .
‚t (˜c· + b •· (·, c· )) ·

If the component · is in competition with other components in the sur-
face reaction, as, e.g., in ion exchange, then f· has to be replaced by a
nonlinearity that depends on the concentrations of all involved components
c1 , . . . , cN , s1 , . . . , sN . Thus we obtain a coupled system in these variables.
Finally, if we encounter homogeneous reactions that take place solely in the
¬‚uid phase, an analogous statement is true for the source term Q· .

0.1 Give a geometric interpretation for the matrix condition of k in (0.16)
and D mech in (0.28).

± ∈ {w, g})
0.2 Consider the two-phase ¬‚ow (with constant ±,

‚t (φS± ) + ∇ · q ± = f± ,
’»± k (∇p± +
q± = ± gez ) ,
Sw + Sg = 1,
pg ’ pw = pc
with coe¬cient functions
± ∈ {w, g}.
pc = pc (Sw ) , »± = »± (Sw ) ,
Starting from equation (0.23), perform a transformation to the new
q = qw + qg , “total ¬‚ow,”

»g ’ »w dpc
1 1
p= (pw + pg ) + dξ , “global pressure,”
2 2 »g + »w dξ

and the water saturation Sw . Derive a representation of the phase ¬‚ows in
the new variables.
14 0. Modelling Processes in Porous Media with Di¬erential Equations

0.3 A frequently employed model for mechanical dispersion is
Dmech = »L |v|2 Pv + »T |v|2 (I ’ Pv )
with parameters »L > »T , where v = q/˜ and Pv = vv T /|v|2 . Here2
»L and »T are the longitudinal and transversal dispersion lengths. Give a
geometrical interpretation.

0.5 Boundary and Initial Value Problems
The di¬erential equations that we derived in Sections 0.3 and 0.4 have the
common form
‚t S(u) + ∇ · (C(u) ’ K(∇u)) = Q(u) (0.33)
with a source term S, a convective part C, a di¬usive part K, i.e., a total
¬‚ux C ’ K and a source term Q, which depend linearly or nonlinearly
on the unknown u. For simpli¬cation, we assume u to be a scalar. The
nonlinearities S, C, K, and Q may also depend on x and t, which shall be
suppressed in the notation in the following. Such an equation is said to be
in divergence form or in conservative form; a more general formulation is
obtained by di¬erentiating ∇ · C(u) = ‚u C(u) · ∇u + (∇ · C)(u) or by

introducing a generalized “source term” Q = Q(u, ∇u). Up to now we have
considered di¬erential equations pointwise in x ∈ „¦ (and t ∈ (0, T )) under
the assumption that all occurring functions are well-de¬ned. Due to the
applicability of the integral theorem of Gauss on „¦ ‚ „¦ (cf. (3.10)), the
integral form of the conservation equation follows straightforwardly from
the above:

(C(u) ’ K(∇u)) · ν dσ = Q(u, ∇u) dx
‚t S(u) dx + (0.34)
˜ ˜
„¦ „¦

with the outer unit normal ν (see Theorem 3.8) for a ¬xed time t or also
in t integrated over (0, T ). Indeed, this equation (on the microscopic scale)
is the primary description of the conservation of an extensive quantity:
Changes in time through storage and sources in „¦ are compensated by the
normal ¬‚ux over ‚ „¦. Moreover, for ‚t S, ∇ · (C ’ K), and Q continuous
on the closure of „¦, (0.33) follows from (0.34). If, on the other hand, F is
a hyperplane in „¦ where the material properties may rapidly change, the
jump condition
[(C(u) ’ K(∇u)) · ν] = 0 (0.35)
for a ¬xed unit normal ν on F follows from (0.34), where [ · ] denotes the
di¬erence of the one-sided limits (see Exercise 0.4).
Since the di¬erential equation describes conservation only in general,
it has to be supplemented by initial and boundary conditions in order to
0.5. Boundary and Initial Value Problems 15

specify a particular situation where a unique solution is expected. Boundary
conditions are speci¬cations on ‚„¦, where ν denotes the outer unit normal
• of the normal component of the ¬‚ux (inwards):
’ (C(u) ’ K(∇u)) · ν = g1 on “1 (0.36)
(¬‚ux boundary condition),
• of a linear combination of the normal ¬‚ux and the unknown itself:
’ (C(u) ’ K(∇u)) · ν + ±u = g2 on “2 (0.37)
(mixed boundary condition),
• of the unknown itself:
u = g3 on “3 (0.38)
(Dirichlet boundary condition).
Here “1 , “2 , “3 form a disjoint decomposition of ‚„¦:
‚„¦ = “1 ∪ “2 ∪ “3 , (0.39)
where “3 is supposed to be a closed subset of ‚„¦. The inhomogeneities
gi and the factor ± in general depend on x ∈ „¦, and for nonstationary
problems (where S(u) = 0 holds) on t ∈ (0, T ). The boundary conditions
are linear if the gi do not depend (nonlinearly) on u (see below). If the gi
are zero, we speak of homogeneous, otherwise of inhomogeneous, boundary
Thus the pointwise formulation of a nonstationary equation (where S
does not vanish) requires the validity of the equation in the space-time
QT := „¦ — (0, T )
and the boundary conditions on the lateral surface of the space-time
ST := ‚„¦ — (0, T ) .
Di¬erent types of boundary conditions are possible with decompositions
of the type (0.39). Additionally, an initial condition on the bottom of the
space-time cylinder is necessary:
for x ∈ „¦ .
S(u(x, 0)) = S0 (x) (0.40)
These are so-called initial-boundary value problems; for stationary prob-
lems we speak of boundary value problems. As shown in (0.34) and (0.35)
¬‚ux boundary conditions have a natural relationship with the di¬erential
equation (0.33). For a linear di¬usive part K(∇u) = K∇u alternatively
we may require
‚νK u := K∇u · ν = g1 on “1 , (0.41)
16 0. Modelling Processes in Porous Media with Di¬erential Equations

and an analogous mixed boundary condition. This boundary condition is
the so-called Neumann boundary condition. Since K is symmetric, ‚νK u =
∇u · Kν holds; i.e., ‚νK u is the derivative in direction of the conormal Kν.
For the special case K = I the normal derivative is given.
In contrast to ordinary di¬erential equations, there is hardly any general
theory of partial di¬erential equations. In fact, we have to distinguish dif-
ferent types of di¬erential equations according to the various described
physical phenomena. These determine, as discussed, di¬erent (initial-)
boundary value speci¬cations to render the problem well-posed. Well-
posedness means that the problem possesses a unique solution (with certain
properties yet to be de¬ned) that depends continuously (in appropriate
norms) on the data of the problem, in particular on the (initial and)
boundary values. There exist also ill-posed boundary value problems for
partial di¬erential equations, which correspond to physical and technical
applications. They require special techniques and shall not be treated here.
The classi¬cation into di¬erent types is simple if the problem is lin-
ear and the di¬erential equation is of second order as in (0.33). By order
we mean the highest order of the derivative with respect to the variables
(x1 , . . . , xd , t) that appears, where the time derivative is considered to be
like a spatial derivative. Almost all di¬erential equations treated in this
book will be of second order, although important models in elasticity the-
ory are of fourth order or certain transport phenomena are modelled by
systems of ¬rst order.
The di¬erential equation (0.33) is generally nonlinear due to the nonlin-
ear relationships S, C, K, and Q. Such an equation is called quasilinear if
all derivatives of the highest order are linear, i.e., we have
K(∇u) = K∇u (0.42)
with a matrix K, which may also depend (nonlinearly) on x, t, and u.
Furthermore, (0.33) is called semilinear if nonlinearities are present only
in u, but not in the derivatives, i.e., if in addition to (0.42) with K being
independent of u, we have
S(u) = Su , C(u) = uc (0.43)
with scalar and vectorial functions S and c, respectively, which may depend
on x and t. Such variable factors standing before u or di¬erential terms are
called coe¬cients in general.
Finally, the di¬erential equation is linear if we have, in addition to the
above requirements,
Q(u) = ’ru + f
with functions r and f of x and t.
In the case f = 0 the linear di¬erential equation is termed homoge-
neous, otherwise inhomogeneous. A linear di¬erential equation obeys the
superposition principle: Suppose u1 and u2 are solutions of (0.33) with the
0.5. Boundary and Initial Value Problems 17

source terms f1 and f2 and otherwise identical coe¬cient functions. Then
u1 + γu2 is a solution of the same di¬erential equation with the source
term f1 + γf2 for arbitrary γ ∈ R. The same holds for linear boundary
conditions. The term solution of an (initial-) boundary value problem is
used here in a classical sense, yet to be speci¬ed, where all the quantities
occurring should satisfy pointwise certain regularity conditions (see De¬ni-
tion 1.1 for the Poisson equation). However, for variational solutions (see
De¬nition 2.2), which are appropriate in the framework of ¬nite element
methods, the above statements are also valid.
Linear di¬erential equations of second order in two variables (x, y) (in-
cluding possibly the time variable) can be classi¬ed in di¬erent types as
To the homogeneous di¬erential equation
‚2 ‚2 ‚2
Lu = a(x, y) 2 u + b(x, y) u + c(x, y) 2 u
‚x ‚x‚y ‚y (0.44)
‚ ‚
+ d(x, y) u + e(x, y) u + f (x, y)u = 0
‚x ‚y
the following quadratic form is assigned:
(ξ, ·) ’ a(x, y)ξ 2 + b(x, y)ξ· + c(x, y)· 2 . (0.45)
According to its eigenvalues, i.e., the eigenvalues of the matrix
a(x, y) 2 b(x, y)
, (0.46)
2 b(x, y) c(x, y)
we classify the types. In analogy with the classi¬cation of conic sections,
which are described by (0.45) (for ¬xed (x, y)), the di¬erential equation
(0.44) is called at the point (x, y)
• elliptic if the eigenvalues of (0.46) are not 0 and have the same sign,
• hyperbolic if one eigenvalue is positive and the other is negative,
• parabolic if exactly one eigenvalue is equal to 0.
For the corresponding generalization of the terms for d + 1 variables and
arbitrary order, the stationary boundary value problems we treat in this
book will be elliptic, of second order, and ” except in Chapter 8 ” also
linear; the nonstationary initial-boundary value problems will be parabolic.
Systems of hyperbolic di¬erential equations of ¬rst order require partic-
ular approaches, which are beyond the scope of this book. Nevertheless,
we dedicate Chapter 9 to convection-dominated problems, i.e., elliptic or
parabolic problems close to the hyperbolic limit case.
The di¬erent discretization strategies are based on various formulations
of the (initial-) boundary value problems: The ¬nite di¬erence method,
which is presented in Section 1, and further outlined for nonstationary prob-
lems in Chapter 7, has the pointwise formulation of (0.33), (0.36)“(0.38)
18 0. Modelling Processes in Porous Media with Di¬erential Equations

(and (0.40)) as a starting point. The ¬nite element method, , which lies in
the focus of our book (Chapters 2, 3, and 7), is based on an integral formu-
lation of (0.33) (which we still have to depict) that incorporates (0.36) and
(0.37). The conditions (0.38) and (0.40) have to be enforced additionally.
Finally, the ¬nite volume method (Chapters 6 and 7) will be derived from
the integral formulation (0.34), where also initial and boundary conditions
come along as in the ¬nite element approach.

0.4 Derive (formally) (0.35) from (0.34).

0.5 Derive the orders of the given di¬erential operators and di¬er-
ential equations, and decide in every case whether the operator is
linear or nonlinear, and whether the linear equation is homogeneous or
(a) Lu := uxx + xuy ,
(b) Lu := ux + uuy ,

(c) Lu := 1 + x2 (cos y)ux + uyxy ’ arctan x u = ln(x2 + y 2 ) ,

(d) Lu := ut + uxxxx + 1 + u = 0 ,
(e) utt ’ uxx + x2 = 0 .

0.6 (a) Determine the type of the given di¬erential operator:
(i) Lu := uxx ’ uxy + 2uy + uyy ’ 3uyx + 4u ,
(ii) Lu = 9uxx + 6uxy + uyy + ux .
(b) Determine the parts of the plane where the di¬erential operator Lu :=
yuxx ’ 2uxy + xuyy is elliptic, hyperbolic, or parabolic.
(c) (i) Determine the type of Lu := 3uy + uxy .
(ii) Compute the general solution of Lu = 0.

0.7 Consider the equation Lu = f with a linear di¬erential operator of
second order, de¬ned for functions in d variables (d ∈ N) in x ∈ „¦ ‚ Rd .
The transformation ¦ : „¦ ’ „¦ ‚ Rd has a continuously di¬erentiable,

<< . .

. 3
( : 59)

. . >>