ńņš. 3 |

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

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Ī· ,

(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

reads

ā‚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

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

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