Fowler A. Mathematical Geoscience

Подождите немного. Документ загружается.

430 7 Groundwater Flow

Fig. 7.9 Data from borehole 102 at Rexco, May 2003, courtesy of David Lerner and Arnë

Huttmann, GPRG, University of Shefﬁeld. Units are metres for depth, mg l

−1

for concentrations.

Since the molecular weights of ammonium and nitrate are 18 and 62 (g mole

−1

), respectively, then

1mgNH

−1

mmol NH

−1

,1mgNO

−

−1

mmol NO

−

−1

. Nitrate shows a sharp

spike (of about 10 mmol l

−1

) at a depth of 19 m. There appears to be a second front at 23 m: the

nitrate produced at 19 m diffuses there and takes out (according to (7.201)) either the inorganic

carbon CH

O or the acid H

where the parameters ν

are deﬁned analogously to (7.194), i.e.,

ˆc

, (7.205)

ˆc

being the scale for c

. From the values given above, we can suppose that all the

 1.

We now mirror the Four Ashes discussion, supposing that Λ  1. The reaction

front at the plume fringe should then be thin as before, with oxygen diffusing to

the front from outside the plume, and ammonium diffusing to it from inside, and

both concentrations being zero at the front. The fringe position is unknown, and we

seek a conservation law to provide an extra condition for it analogous to (7.197).

Denote the right hand sides of (7.202)bye

–e

. There are ﬁve reactants but only

three reactions, therefore there are two (linear) relationships between the e

: these

provide suitable jump conditions across the fringe. The equations



=0are

satisﬁed for any β

and β

provided β

=−β

, β

=β

, and β

=−β

Selecting (β

,β

) =(1, 0) and (0, 1), we thus have

−e

=0,

−e

=0,

(7.206)

and it follows from this that in consideration of (7.204) there are two conserved

quantities across the fringe (i.e., ﬂux in equals ﬂux out); explicitly (with the same

7.8 Three Speciﬁc Remediation Problems 431

caveat concerning the normal n as expressed in (7.199)), we have



∂c

∂n

−ν

−

∂c

−

∂n

+ν

∂c

∂n



−

=0,



∂c

∂n

−

∂c

−

∂n

∂c

∂n

−ν

∂c

∂n



−

=0,

(7.207)

where [j ]

−

denotes the jump in j across the fringe.

Consulting (7.201), we can write the dimensionless reaction rates r

in terms of

Monod rates

(7.208)

in the following way:

−

(7.209)

where k

and k

are dimensionless constants (the ratios of the saturated maxima

max

of reactions 2 and 3 in (7.201) to that of reaction 1). The requirement that re-

action rates vanish on either side of the fringe is satisﬁed if c

= 0 outside the

plume (then r

=0), c

=0 within the plume (so r

=0) and if c

−

=0on

both sides (then r

=0). This simpliﬁes the situation to one where pre-existing

soil nitrate concentrations are low, and assumes all the nitrate produced in the re-

action front by reaction 1 is consumed there by reaction 2. An alternative would be

that the produced nitrate diffuse away on either side, but this is not consistent with

a fast reaction 2 unless there is a natural source for nitrate production (e.g., from

surface composting).

Since from (7.203), ν

∼ 10

−5

, this suggests that

∂c

∂n

is approximately continuous across the fringe, which is thus diffusive in character

(providing Pe 1). In that case the fringe is located simply by advection of ammo-

nium.

The second jump condition then provides a ﬂux boundary condition for c

which is not constrained to be zero at the front. In addition, existence of a spike of ni-

trate within the front must be determined by solution of the reaction equations within

the front. The reaction front is of thickness ∼(1/ΛPe)

1/2

(taking Pe



=Pe

⊥

for sim-

plicity), and within this front, the nitrate concentration is given approximately by

∂

−

∂N

+ν

−

−r

−0.6r

) =0, (7.210)

Alternatively, there may be more than one reaction front. This is suggested by Fig. 7.9.

432 7 Groundwater Flow

where N is a suitably rescaled normal variable. We can expect this to be solv-

able for c

−

subject to c

−

→ 0asN →±∞, since for example if we lin-

earise the nitrate Monod coefﬁcient, M

−

, then the solution can be written

using a Green’s function which depends also on all the other reactant concentra-

tions.

There remains the issue of determining a boundary condition for H

at the front.

Counting conditions, the ﬁrst ﬂux condition in (7.207) determines the location of the

fringe. We already have boundary conditions for c

and c

(both equal zero),

and the nitrate spike is determined from (7.210). Thus, in principle we know the

jump in ﬂux of ammonium, oxygen and nitrate, and by integrating (7.210) and its

equivalents for ammonium and oxygen through the front, we know that



∞

−∞

)dN =



∂c

−

∂N



−



∞

−∞

dN =



∂c

∂N



−

, (7.211)

−



∞

−∞

−(r

−r

−0.6r

)dN =0.

Thus we know the values of



∞

−∞

dN for i =1, 2, 3 and this tells us (by integrating

through the reaction front) the values of



∞

−∞

dN for i =4, 5, and thus the jump in

ﬂux of

∂c

∂N

; this provides the extra boundary condition we seek. The jump in ﬂux

of c

is also given by



∞

−∞

1.25ν

dN, but this is equivalent to (7.207)

.In

this way, the approximate model provides all the conditions necessary to determine

the solution.

St. Alban’s At St. Alban’s, there has been a petroleum spillage (at a ﬁlling station)

into an underlying chalk aquifer. The ﬂuids are LNAPLs: hydrocarbons, BTEX (a

cancer-forming aromatic hydrocarbon

and MTBE.

BTEX is retarded compared

to MTBE and thus forms a secondary plume within the MTBE plume. The LNAPLs

have seeped through the unsaturated zone and sit on top of the chalk aquifer, act-

ing as a source (via dissolution) of contaminant to the underlying groundwater

ﬂow.

The sequence of reactions appears to be similar to those of the other examples,

with oxidation by oxygen and nitrate at the plume fringe, and by Mn (manganese),

Fe (iron) and SO

2−

(sulphate) in the plume core. The sequence of reactions which

More speciﬁcally, BTEX refers to a suite of volatile hydrocarbons, the acronym referring to ben-

zene, toluene, ethylbenzene and xylene, with chemical formulae C

(benzene), C

(toluene),

(ethylbenzene and xylene); we use toluene in the chemical reaction model.

Methy tert-butyl ether, C

7.8 Three Speciﬁc Remediation Problems 433

are modelled is the following:

+9O

+3H

→ 7CO

2−

+7.5H

→ 2.5CO

2−

+4.5CH

+5H

+36Fe

→ 7CO

2−

+36Fe

+50H

+2O

→ CO

2−

+2H

+8Fe

+3H

→ CO

2−

+8Fe

+10H

→ Fe

+0.5H

+36FeOOH

(s)

+58H

→ 36Fe

+7CO

2−

+51H

+8FeOOH

(s)

+7H

→ 8Fe

+CO

2−

+5H

(7.212)

Maximum reaction rates vary over about four orders of magnitude, from about

−12

mol l

−1

(reactions 2, 3, 7) to 10

−8

mol l

−1

(reaction 8). The princi-

pal oxidising reaction rate r

∼10

−10

mol l

−1

. Contaminant levels are of order

−7

mol l

−1

, with oxygen level of order 10

−6

mol l

−1

. Plume depth and length are

d ∼ 20 m, l ∼ 300 m, dispersivity parameters are taken as α



∼ 1m,α

⊥

∼ 0.2m.

The hydraulic conductivity of the upper chalk aquifer is 25 m d

−1

(that of the lower

part is only about 0.1 m d

−1

), and the hydraulic gradient is 1.75 ×10

−3

, thus the

longitudinal ﬂux scale is U ∼16 m y

−1

. With these values, the parameters deﬁned

by (7.186) have values



∼300, Pe

⊥

∼10,Λ∼10

. (7.213)

As now seems monotonously to be the case, reactions are fast, longitudinal disper-

sion is small, but transverse dispersion may be effective because of the high aspect

ratio l/d.

The distinguishing feature of this particular site is that, although the chalk is very

porous (φ  0.3), it has a very low effective permeability, presumably due to chem-

ical adsorption by the chalk. On the other hand, the chalk matrix is dissected into

blocks by numerous fractures, and thus acts as a dual porosity system. Contaminant

can diffuse into the pore space of the matrix, and the issue of concern is whether and

how fast this happens, since storage in matrix blocks will act as a residual source of

contamination after the fracture system has been ﬂushed.

We have considered this problem before, in Sect. 7.6.1. Let us suppose for sim-

plicity that dispersivities of matrix and fractures are constant and isotropic but not

necessarily equal. (Additionally, we suppose the local fracture dispersivity D

equal to the macroscopic fracture dispersivity D

.) The point forms of the equations

describing a single reactant are given by (7.165), and the block scale and averaged

434 7 Groundwater Flow

equations for the matrix concentration are given by (7.167) and (7.170):



∂ ¯c

∂t

.∇

¯c



=∇

¯c

+Pe

n.∇

|

∂M

εPe



∂c

∂T

.∇



=∇

+ε

ΛS

(7.214)

where

,Λ=

. (7.215)

Analogous dimensionless forms for the fracture average and block scale equations

are, from (7.162) and (7.165)



∂ ¯c

∂t

.∇

¯c



=∇

¯c

+Pe

−

n.∇

|

∂M

εPe



∂c

∂T

.∇



=∇

+ε

ΛS

−

(n.∇

∂M

(7.216)

and Pe

=Ul/D

The question is how the interfacial transport term should be modelled. In princi-

ple we solve the block equation for c

with c

on ∂M, yielding the interfacial

term as an integral convolution of c

. Putting this into the block equation for c

and

solving this then determines c

and thus gives the interfacial term, which closes the

description of the averaged equations.

The upshot of our earlier discussion was that the details of the homogenisation

process depend on the relation between the parameters ε = d

/l, Pe and Λ.Clas-

sical homogenisation theory as in Sect. 7.6.1 assumes all the parameters are O(1)

apart from ε, but this is unlikely ever to be appropriate. Estimates for St. Alban’s

may be Pe

∼10

, Λ ∼10

, ε ∼10

−4

; in addition, u

u

and D

D

. Thus

the reaction terms in the macroscopic equations are always large, and this suggests

that, as before, reactions will be restricted to thin fronts. A complication in (7.214)

is that the interfacial transport term may be large also, so this now needs to be de-

termined.

To be speciﬁc, let us consider a simple two component reaction similar to (7.191),

i.e.,

+σ O

→products, (7.217)

and let p denote dimensionless BTEX concentration and c denote dimensionless

oxygen concentration. Just as in (7.192), we can put S

=νS in (7.221), and there

is a corresponding reaction term S in the equivalent equation for p: S can be taken

as the product of Monod rates in (7.193), and for simplicity we take the linear rates,

thus S =−pc.

7.8 Three Speciﬁc Remediation Problems 435

The critical parameter in (7.214)

is ε

Λ; we write

=ε

Λ, B

=ε

Λ; (7.218)

if Pe

∼10

then B

∼O(1), and B

1: reaction in the block is fast.

Just as in (7.196), c

−νp

satisﬁes an advective-diffusion conservation equa-

tion in the block, diffusion acting on a time scale t ∼ε

, where t is the macro-

scopic time variable. The parameter ε

is crucial. If we suppose that molec-

ular diffusion applies in the blocks, then D

∼ 10

−9

−1

, Pe

∼ 10

, and

∼ 10

−3

is small. Therefore on the long macroscopic time scale, diffusion

smoothes c

−νp

, and we can take it to be locally constant in the blocks (and thus

also the fractures). Since reaction is fast in the blocks, this implies

≈νp

, (7.219)

and therefore the block reaction rate for c

can be written as

ΛS

≈−B

. (7.220)

When B

is large there is a reaction boundary layer like a rind at the block

surface, within which c

satisﬁes

∂

∂N

−B

≈0. (7.221)

In the interior of a block inside the plume the oxygen is depleted. and we have

→0asN →−∞. The ﬁrst integral of (7.221) determines the ﬂux to the blocks

∂c

∂N



3/2

, (7.222)

using the fact that c

on ∂M. The matrix concentration is negligible, and the

average fracture concentration satisﬁes the Eq. (7.216)



∂ ¯c

∂t

.∇

¯c



=∇

¯c

−

¯c

−



3Pe

¯c

3/2

, (7.223)

where we take S

=−c

as in the matrix, and we have supposed that

=¯c

, c

3/2

=¯c

3/2

. (7.224)

With B

∼O(1), Pe

∼10

, Pe

∼10

, the ﬂux term to the blocks,

¯c

1 and

the fracture reaction term



3Pe

¯c

3/2

1; the implication is that the fracture

concentration of oxygen within the plume rapidly decreases within both blocks and

436 7 Groundwater Flow

fractures to very small levels. The dual porosity appears to have little effect on the

characteristics of the reactant distributions.

7.9 Precipitation and Dissolution

A sedimentary basin, as illustrated in Fig. 7.10, refers to an accumulation of sedi-

ments (of typical depth 10 km) derived from river outwash sands and silts, or marine

deposited muds and microfossils. Sedimentary basins are everywhere in continen-

tal rocks: we have the Paris basin, the London basin, the North Sea, and so on.

As sediments accumulate and are buried, they are subjected to increasing heat and

pressure, and these two factors enable the formation of intergranular cements, which

thus convert the sediments to rock (sand to sandstone; clay to shale, marine organ-

isms to limestone). As they are buried, the sediments also compact, expelling pore

water. Depending on the permeability, this can lead to pore pressures above hydro-

static, a situation which is of concern in oil-drilling operations (and is discussed

further in Sect. 7.11).

Diagenesis refers generally to the process of chemical alteration to rock, and in

this section we study the effects of diagenesis on the groundwater ﬂow of accumu-

lating sediments within a sedimentary basin. The particular type of diagenesis which

we discuss is the conversion of smectite (a form of hydrated clay) to illite (a dehy-

drated clay) via a dewatering reaction. The resultant release of water is also a poten-

tial cause of excess pore pressures, but the main purpose of the present discussion

is to show how the use of an approximation which we may call the weak solubility

limit allows enormous simpliﬁcation of quite complicated reaction schemes.

One view of the smectite–illite reaction is to treat it using ﬁrst order kinetics, thus

→I

+nH

O, (7.225)

where each mole conversion yields n moles of water: S denotes smectite, I denotes

illite, and the superscript

denotes the solid phase (likewise,

will denote an aque-

ous phase). Such a scheme is not inconsistent with at least some experimental data,

and the rate factor involved depends on temperature, with an activation energy in

the range 60–80 kJ mol

−1

However, it is likely that the transformation of smectite to illite occurs through

a compound sequence of precipitation and dissolution reactions; one possible de-

Fig. 7.10 Schematic of a

sedimentary basin. Sediment

accumulates from outﬂow

from rivers, and also through

the settlement of marine

organisms

7.9 Precipitation and Dissolution 437

scription is the following:

−→ X

+nH

KFs

−→ K

+AlO

−L

+s SiO

+AlO

−L

+f X

−→ f I

+SiO

SiO



−

Qz.

(7.226)

The smectite S

dissolves to form a hydrous silica combination X

, such as

(OH)

. Additionally, illite precipitation requires potassium ions, and these

may be obtained from the dissolution of potassium feldspar (in the second reaction);

the aluminium hydroxyl ions AlO

−L

act in the same way. The hydrous silica com-

bination now combines with the potassium and aluminium to form illite precipitate,

together with aqueous silica SiO

, which itself precipitates as quartz Qz. Taking

suitable multiples of the reactions and adding to eliminate the aqueous phases, the

overall reaction is found to be

S +f

−1

KFs

−→ I +nH

O +f

−1

(s +1)Qz. (7.227)

Ideally we would like to be able to write kinetics for (7.227) analogously to (7.225),

with a recipe for the effective reaction rate R. Note that all the reactions in (7.226)

are precipitation or dissolution reactions, and therefore the reaction rates are pro-

portional to grain surface area.

It turns out, at least for this reaction scheme (but we might suspect more gener-

ally), that the weak solubility limit allows such a recipe to be found. To illustrate the

method, note that conservation equations for the concentrations S, X, F , K, A, L,

Q, I of the substances smectite, hydrous silica, feldspar, potassium ions, aluminium

hydroxyl ions, aqueous silica, quartz and illite satisfy equations of the type

∂

∂t



(1 −φ)S



+∇.



(1 −φ)Su



=−r

∂

∂t

(φX) +∇.



φXu



−∇.(φD.∇X) =r

−fr

(7.228)

and so on, where the reaction rates r

are scaled with respect to the surface rates

by the speciﬁc interfacial surface areas Σ

(thus r

= Σ

). u

and u

are the

liquid and solid Darcy ﬂuxes (i.e., the volume ﬂuxes per unit area). We allow a

non-zero solid velocity in order to cater for the effects of compaction of the porous

matrix. There are six other equations of the type in (7.228), together with a water

conservation equation (this is the equation for φ). Of the total of nine equations,

four are for aqueous concentrations. That for A is identical to that for K, and we

ignore it henceforth.

438 7 Groundwater Flow

The weak solubility limit is associated with the observations that the aqueous

dissolved species X (hydrous silica), K (potassium) and L (silica) are typically

present in trace quantities of the order of 10–100 ppm (1 ppm = 10

−3

kg m

−3

and thus the concentrations of the aqueous phases are much less than those of the

solid phases. When the model is suitably non-dimensionalised, the result is that the

transport terms for the aqueous phases are very small, so that the corresponding

reaction terms can be taken to be in equilibrium. The reaction terms in the equations

for X, K and L are, respectively, r

−fr

(as above), r

−r

and sr

−r

−

From this, we obtain the three relationships

≈ fr

≈ r

, (7.229)

−r

−

≈ (s +1)r

Since the rate of dissolution of S is r

and the rate of precipitation of I is fr

this immediately shows that ﬁrst order kinetics of the form (7.225) does apply, with

the reaction rate being fr

. Since r

is a precipitation rate, and r

−

the dissolution

rate of the same mineral pair SiO

↔ Qz, and since either r

=0orr

−

=0, it is

apparent that (7.229)

determines r

and r

−

together. Thus all the reaction rates

can be written in terms of r

, and this then determines the aqueous phase pseudo-

equilibrium concentrations of X, L and K in terms of r

and various temperature

dependent rate factors, since the kinetic rates r

are prescribed in terms of these.

If the smectite equation (7.228)

is written in terms of smectite volume fraction

, then it becomes

∂φ

∂t

+∇.





=−

, (7.230)

where M

is the molecular weight of smectite, and ρ

its density. An equivalent

equation for illite is

∂φ

∂t

+∇.





, (7.231)

and there are two similar equations for φ

and φ

(with the right hand sides being

proportional to r

). In addition, the porosity φ satisﬁes

∂φ

∂t

+∇.



φu



. (7.232)

It remains to determine the rate constant r

in terms of the reactant concentrations.

Actually, 1 ppm =1mgkg

−1

,but1m

of water weighs 10

kg, so for aqueous solutions, this is

equivalent to 10

−3

kg m

−3

7.9 Precipitation and Dissolution 439

We assume that the reaction rates take the general form (in which D denotes

dissolution, and P precipitation)

=Σ



1 −



=Σ



−1



(7.233)

where c

is the relevant aqueous concentration of phase i and c

is the associated

solubility limit, i.e., the saturation concentration of aqueous phase i in the presence

of solid phase s. The rate factor R

generally depends on temperature. (7.233) states

that precipitation occurs when a solution is oversaturated, and dissolution occurs if

it is undersaturated.

For the speciﬁc case of the precipitation/dissolution scheme in (7.226), we sup-

pose that

=φ

[1 −ψ

]

=φ

[1 −ψ

]

[θ

−ψ

]

=φ

[ψ

−θ

]

[ψ

−θ

]

=φ

[ψ

−1]

=φ

[1 −ψ

]

(7.234)

In these equations we have assumed that speciﬁc surface area is proportional to

volume fraction, thus Σ

= φ

, and we have absorbed the grain size d

into

the reaction rates. The reaction rate subscripts describe what they represent: SD,

smectite dissolution; FD, feldspar dissolution; IP, illite precipitation; QP, quartz

precipitation; QD, quartz dissolution.

The quantities ψ

represent dimensionless aqueous concentrations (of phase D),

scaled with a suitable solubility limit: c

is scaled with c

, the solubility limit for

hydrous silica in the presence of smectite; c

is scaled with c

, the solubility

limit for silica in the presence of quartz. Potassium is more complicated, because

in (7.226), we see that dissolution of feldspar produces two aqueous phases, and

precipitation of illite requires two. We can in fact deﬁne four further solubility limits:

, that of silica in the presence of feldspar; c

, that of potassium in the presence

of feldspar; c

, that of potassium in the presence of illite; and c

, that of hydrous

silica in the presence of illite. The assumption is then that, for example, feldspar will

dissolve only if both potassium and silica are undersaturated, i.e., c

and

. In writing (7.234), we have scaled c

with c

, and the three solubility

ratios in (7.234) are therefore deﬁned by

,θ

. (7.235)