<< Ļšåä. ńņš. ńņš. 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
on

ā Ā· 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

:= 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Ī· ,
Ė
Ė
where
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
ā‚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Ī· ,
(2)
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-
e
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
linear.
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
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
(1)
ļ¬‚uid phase, an analogous statement is true for the source term QĪ· .

Exercises
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
variables
q = qw + qg , ātotal ļ¬‚ow,ā

Ī»g ā’ Ī»w dpc
S
1 1
p= (pw + pg ) + dĪ¾ , āglobal pressure,ā
2 2 Ī»g + Ī»w dĪ¾
Sc

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
conditions.
Thus the pointwise formulation of a nonstationary equation (where S
does not vanish) requires the validity of the equation in the space-time
cylinder
QT := ā„¦ Ć— (0, T )
and the boundary conditions on the lateral surface of the space-time
cylinder
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
follows:
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
1
a(x, y) 2 b(x, y)
, (0.46)
1
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.

Exercises
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
inhomogeneous:
(a) Lu := uxx + xuy ,
(b) Lu := ux + uuy ,
ā
(c) Lu := 1 + x2 (cos y)ux + uyxy ā’ arctan x u = ln(x2 + y 2 ) ,
y
ā
(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)ĪĆĖĄĀĖÅĶČÅ Ńėåä. ńņš. >>