FP-Experiment: Water and solute transport in porous media

23 Pages·2008·0.719 MB·English

by  Klaus Schneider, Ute Wollschläger, and Kurt Roth

FP-Experiment
Water and solute transport in porous media
Klaus Schneider, Ute Wollschläger, and Kurt Roth
Institut für Umweltphysik, Universität Heidelberg
Version: April 2008
Location of the experiment: INF 229, room 204
Contents
1 Overview 2
2 Part I: Estimation of hydraulic properties by inverse modelling 3
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.3 Dynamics of water movement . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.3.1 Hydraulic properties . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.3.2 Experimental determination of hydraulic properties by inverse mod-
elling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 Multi-Step Outflow Experiment . . . . . . . . . . . . . . . . . . . . . . . . 11
2.4.1 Setup of the experiment . . . . . . . . . . . . . . . . . . . . . . . . 11
2.4.2 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4.3 Schedule of the practical course (part I) . . . . . . . . . . . . . . . 13
2.5 What to learn in this experiment . . . . . . . . . . . . . . . . . . . . . . . 15
2.6 Questions you should be able to answer before starting . . . . . . . . . . . 16
3 Part II: Measurement and parameterisation of solute transport 17
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2.1 Convection-Dispersion model (CD) . . . . . . . . . . . . . . . . . . 17
3.2.2 Extension to multi-domain models (MIM) . . . . . . . . . . . . . . 20
3.3 Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3.1 Setup of the experiment . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3.2 Schedule of the practical course (part II) . . . . . . . . . . . . . . . 21
3.4 What to learn in this experiment . . . . . . . . . . . . . . . . . . . . . . . 23
3.5 Questions you should be able to answer before starting . . . . . . . . . . . 23
1 Overview
Background: Porous media are part of our daily life and include technical devices like
filter materials and fuel cells as well as natural systems like geologic formations and soils.
An important characteristic of porous media are capillary forces that play a dominant role
for the dynamics of such materials. It is a challenge for physicists to model fluid flow and
transport of dissolved chemicals in porous media e.g. in order to predict the behaviour
of a system as response to various boundary conditions or to construct porous media
with the properties desired. The focus here is on soil, which is typically characterised by a
complex structure. It consists of minerals, organic material and the pore space in-between.
Moreover, this structure is typically not static but continuously changed by plants, animals
and various physical boundary conditions due to temperature and precipitation. Crucial
questions that require a quantitative understanding of the processes in porous media are
forexample: Howtostoreinfiltratingwater insoiltowarrant thegrowthofplants? Howto
prevent pollution in the environment of waste disposals? How to treat soils contaminated
by anthropogenic wastes? What is the role of soils for global mass and energy fluxes?
Subject of the course: In principle, the microscopic structure of the pore space deter-
mines the behaviour of porous media. However we are typically not in the position to
describe this 3-dimensional structure in all detail. We aim at modelling and predicting
fluid flow and solute transport in porous media at a scale of some 10−1m within a soil
up to some 103m within a catchment. This range of scales is significantly larger than
the typical size of pores in soil that are in the range of 10−6m to 10−3m. (Just think
about the amount of computer memory you would need to store the structure of 1dm3
soil at a resolution of 1µm). This leads to the introduction of “effective” models based on
“effective” material properties that translate the effects of the microscopic structure into
macroscopic properties.
In the first part of this course we only consider water flow in porous media and we intro-
duce two effective material properties: 1) The relation between the energy density of water
(capillary pressure) and the water saturation which is referred to as water characteristic
and 2) the relation between the water saturation of the media and its hydraulic conduc-
tivity which is termed the hydraulic conductivity function. An experiment is performed to
determine these basic material properties which are essential for modelling the dynamics
of fluids in porous media at larger scales.
In the second part we consider solute transport through porous media which also may
be related to macroscopic material properties as for instance a dispersion tensor. Another
experiment is performed to identify a convenient transport model and to determine the
related parameters.
Preparation: To prepare this course, this short script should be sufficient. Here, the
basic theory of water dynamics in porous media is described and some practical de-
tails on the set-up of the experiment are provided. Finally the evaluation of the mea-
sured data through inverse modelling is described. For a deeper study of the subject
you may consider Roth (2006) at www.iup.uni-heidelberg.de/institut/forschung/
groups/ts/soil_physics/students/lecture_notes05.
2
1cm
(a) (b)
(c)
Figure 1: 2D and 3D structure of soil at the pore scale with a resolution of
0.04mm/pixel obtained from serial sections (a and b). Structure of porosity
including few large pores (black spots) at a larger scale with a resolution of
0.4mm/pixel obtained by X-ray tomography (c).
2 Part I: Estimation of hydraulic properties by inverse
modelling
2.1 Introduction
To model the dynamics of water in porous media at a scale much larger thanthe individual
pores, we introduce so called “effective material properties” which are based on variables
that are visible at the larger scale. These properties and the related variables are obtained
by averaging over the complicated porous structure at the microscopic scale. Thereby, the
details of the porous structure are lost, but the aim is to preserve those properties that
are important to quantitatively understand the processes at the larger scale – this is why
we call them “effective”.
At the pore scale (Figure 1b), flow and transport processes may be described using clas-
sical concepts of fluid dynamics where the acceleration of a fluid element is related to the
forces acting on it, most importantly to gravity, pressure gradients and friction. This leads
to the well-known Navier-Stokes equation. This equation is not invoked here, however, be-
cause (i) the inertia term can be neglected for the low flow velocities typically encountered,
3
(ii) implementing the boundary conditions for the pore space is next to impossible and
the computational effort is excessive, and (iii) for most phenomena a detailed pore-space
description is not required. Hence, we go for a continuum description that is obtained from
averaging. In formal analogy to the averaging used in thermodynamics, the dynamics of
water flow is described over a representative region of the pore space. The remaining task
then is to estimate pertinent effective material properties.
2.2 Theory
The complicated structure of the soil cube in Figure 1b is replaced by “effective” measures
and state variables with respect to the fluids. This can be interpreted as averaging over
microscopic properties or upscaling from microscopic to macroscopic properties.
Consider a porous medium which consists of the solid and undeformable “soil matrix”,
i.e. the material the soil is made of, and the intermediate pore space. Pores are filled with
a fluid, e.g. water or air. At every point of the pore space, there is precisely one phase:
soil matrix, water, or air.
State Variables
We define the state of a fluid element of a fluid i in a porous medium by the following
state variables: its height z , its pressure p , its temperature T and the concentrations C
i i i ij
of solute chemicals j. In the following, only pure fluids are considered: C = 0.
ij
Fluid Content
The distribution of a phase i (soil matrix, air or water) is averaged over a volume ∆V
around the spatial location r. One obtains the volume fraction
θ := ∆V /∆V . (1)
i i
It specifies the fraction of the volume ∆V around r which is taken up by the phase i.
θ (r) is a continuous function. Note that at the pore scale, every point belongs exactly to
i
one phase, while in contrary, in the macroscopic description, one point belongs to several
phases.
The volume fraction of the pore space is called porosity φ = θ + θ = 1 − θ . This
a w m
corresponds to the water content at complete water saturation (θ = 0, θ = φ). If we
a w
know the density of the solid material ρ (for quartz it is 2.65g/cm3) we can calculate the
s
porosity from bulk density: φ = 1−(ρ /ρ ).
b s
Representative Averaging Volume
The volume ∆V should be chosen such that θ does not depend on the actual size or shape
i
of ∆V, i.e. that it is large enough to contain all the microscopic heterogeneities. If ∆V is
about the size of a pore, θ will hardly differ from the microscopic description. When ∆V
i
is enlarged, the average value changes abruptly, depending crucially on the choice of ∆V,
and eventually becomes stable if ∆V is above some critical value ∆V . This is the point
c
where the averaging begins to make sense.
Such a volume which leads to a representative value of θ is called Representative Ele-
i
mentary Volume (REV). If no REV can be found, a reasonable averaging is not possible.
4
Potential
The density of potential energy ψ (r) of a fluid i is defined by the energy that is necessary
i
to move a unit volume of fluid from the reference state z = z , p = p , T = T to the state
0 0 0
at position r. p is normally chosen to be the ambient air pressure.
0
By convention, z points downwards into the ground. In an isothermal situation with
no dissolved chemicals, only pressure and gravitation contribute to the energy density.
Considering an incompressive medium we obtain
ψ (r) = p (r)−p −(z −z )ρ g . (2)
i i 0 0 i
=:ψim =:ψig
| {z } | {z }
ψ is called matric potential of the fluid i. In unsaturated porous media, this poten-
im
tial describes the energy densities associated with the boundary layers solid-liquid and
liquid-gas, i.e. caused by capillary forces. In mineral soils, the solid soil matrix is typi-
cally wettable so that the solid material is completely covered by a thin water film (i.e.
the energy density of the boundary layer solid-liquid is smaller than that of liquid-gas).
Hence, the matric potential of water ψ is only determined by the interface liquid-gas σ
m wa
(0.0725Jm−2 for pure water at 20◦C). In the following, only water is considered as fluid
and the fluid index i is omitted. Therewith the total energy density below is denoted by
ψ (w for water).
w
At equilibrium, when the free energy is minimal, the surface area of the liquid-gas
interface is also minimal and one can show that this surface is described by the Young-
Laplace Equation (Landau and Lifschitz, 1991)
1 1
ψ = p −p = σ + , (3)
m w a wa
r r
(cid:20) 1 2(cid:21)
where p and p are the pressures in water and air, respectively, and r , r are the principal
w a 1 2
radii of curvature (the two radii of a local elliptic approximation of the curved interface
area). We use the convention that the sign of the radius is negative if it is within the air
phase. In the simplest case of a cylindrical capillary (3) reduces to ψ = 2σ /r. In the
m wa
convention used here, ψ is negative if the fluid is bound (p < p ) and positive if it is
m i 0
free. The energy density ψ [Jm−3] equals the pressure jump across the interface p −p
m w a
[Nm−2].
In soil physics, often the energy density per unit weight is used:
ψ
m
h := (4)
m
ρ g
w
h is called “hydraulic head” which is usually expressed in “cm water column” [cmWC].
The idea is that h can be interpreted as equivalent height of a fluid column; negative values
correspond to a hanging column.
Additional remarks
1. We restrict the description of the matric potential to the most simple case. In
reality the situation is much more complicated. In particular, σ is not a constant
wa
but depends on temperature and on the composition of the water phase (dissolved
chemicals). Moreover, thewettability ofsoilconstituents mayvaryintimeandspace.
5
Finally, σ changes with the thickness of the water layer when going to very dry
wa
conditions.
2. Besides gravity and capillarity, there may be additional forces: osmotic forces due
to gradients in salt concentration, mechanical forces in case the solid material is
deformable (this is the rule especially in clay rich soils), and the air pressure may be
locally increased when air is entrapped by the water phase.
2.3 Dynamics of water movement
Flux law
In porous media, typical pore dimensions and typical velocities are so small compared
to the viscosity of water and the time scale of typical external forcing that there is no
turbulence. With this prerequisite, the volume flux of water j [ms−1] is linearly related
w
to the gradient of the potential energy density ψ of the water,
w
j = −K(θ)∇ψ , (5)
w w
where the tensorial quantity K(θ) is the hydraulic conductivity that depends on the vol-
umetric water content and ψ is given by (2) applied for water. With decreasing water
w
content, less area is available for water transportation. As large pores are drained first and
flow velocities are higher in larger pores, the conductivity rapidly decreases with water
content. K(θ) is a material property. In a homogeneous medium, the tensor K reduces to
a scalar K.
This flux law is called the Buckingham-Darcy law. Notice the formal analogy with flux
laws for heat (Fourier), dissolved substances (Fick), and electric charge (Ohm).
Richards’ Equation
We formulate the conservation of mass of water – actually of volume since water may be
considered as incompressible – as
∂ θ+∇·j = 0 , (6)
t w
where θ is the volumetric water content. Inserting the flux law (5) into the conservation
equation yields Richards’ equation
∂θ ∂ψ
m −∇·[K(θ)∇ψ ] = 0 (7)
w
∂ψ ∂t
m
for the dynamics of water movement.
2.3.1 Hydraulic properties
Soil water characteristic
The matric potential ψ = p − p is related to the curvature of the water-air interface
m w a
via the Young-Laplace Equation (3). This radius cannot be smaller than the radius of the
pore the interface is located in.
6
6.6 mm = (cid:0)1
= (cid:0)38 hPa = (cid:0)16 hPa
= (cid:0)11 hPa = (cid:0)8 hPa
Figure 2: Distribution of different phases (water: grey, air: black) within a 2D
section through a porous media at different matric potentials ψm. Note the
different radii of curvature of the interfaces according to the Young-Laplace
equation (3).
If ψ is made more negative (e.g. by changing ∆p by applying an external pressure),
m
the radius r of the interface must become smaller according to (3). If r gets smaller than
the radius r¯ of the pore i, r < r¯, there cannot be an interface any more, and the pore will
i i
be drained. This point is called the air entry point. This line of thoughts reveals that the
water content θ is a function of the matric potential ψ : θ = θ (ψ ). This relation is
w m w w m
illustrated in Figure 2: water in wettable porous media tends to fill the small pores while
air resides in the large pores.
Figure 3 shows typical examples for a loamy and a sandy soil. The loamy soil has a
higher water content at saturation, θ , because of a higher porosity. Going from complete
s
saturation to drier conditions by lowering the matric potential the water content of the
sandy material decreases very rapidly below a potential of h ≈ −30cmWC. This is
m
7
)
C℄
W )

m 4 m/h℄ 0
[
m [
h K
g((cid:0) 3 og( -5
o l
l y
ntial 2 tivit -10
te 1 u
o d
p n -15
tri
0 Co
a
m 0 0.20 0.40 0 1 2 3 4
water
ontent(cid:18) matri
potentiallog((cid:0)hm [
mWC℄)
(a) Water characteristic (b) Hydraulic conductivity function
Figure 3: Typical hydraulic properties of a loamy soil (solid line) and a sandy
soil (dashed line). The thin dashed line represents a sand with larger grains
and thus, larger pores.
because this sand is a rather homogeneous material where all pore sizes are within a
narrow range. According to the Young-Laplace equation (3) and using (4) we find an
equivalent pore radius for h = −30cmWC of r ≈ 0.05mm. When going from −30cmWC
m
to about −100cmWC the water within these pores is removed. In the loamy soil the water
content decreases more continuously with decreasing water potential, meaning that the
pore size distribution is much broader compared to the sand. Clearly, the shape of the soil
water characteristic is an effective description of the underlying pore structure. However,
it should be noted that not only the pore size distribution determines the shape of θ(ψ )
m
but also the connectivity of the pores.
Note that we assume that there is an instantaneous equilibrium between θ and ψ which
m
is only the case at high water contents and herewith small resistivity to water flow. In the
dry range however, water in adsorbed films is barely mobile and equilibration has to be
reached through vapour diffusion which is a very slow process.
Thecomplicatedconnectivity ofporesinsoilsalsoleadstothephenomenon ofhysteresis.
Tounderstandthisphenomenon,lookatasingleporeinporespace(Figure4a–c). Initially,
the poreissaturated withwater. Now, pressure atthelower end oftheporeiscontinuously
decreased. The water-air interface stays constant until the potential drops below ψ which
1
corresponds to the radius of the smaller capillary according to (3). Reducing the pressure
any further empties the whole pore: the possible radii in the cavity are too large to sustain
the interface as the pressure difference is too high. This leads to Figure 4 (b).
At this point the process is reverted and the pressure is increased again. This does not
lead back to (a); instead, at ψ , the configuration (c) will be found, which has the same
2
curvature radius as (a) but a different water content θ.
The water characteristic is often described by the parameterisation of van Genuchten
(van Genuchten, 1980) which can be fitted to a wide range of different materials:
θ(ψ )−θ
Θ(ψ ) = m r = [1+[αψ ]n]−1+1/n (8)
m m
θ −θ
s r
with the parameters
Θ [-] water saturation [0,1]
8
drainage
wetting
(a) (b) (c)
ψ
m
ψ
1
ψ
2
θ
(d)
Figure 4: Effect of pore connectivity on the hysteretic behaviour of the water
characteristic shown for an idealised pore: at ψ1 > ψm > ψ2 we can find
very different water contents during drainage ((d), blue curve) and wetting
(red curve) at the same matric potential (same interface radius). Subframes
(a)–(c) show the water content of that pore at the different potentials: (a)
shows the pore at ψm > ψ1, sub-frame (b) at ψm < ψ1 and sub-frame (c) at
ψ1 > ψm > ψ2 after re-irrigation, with the same curvature radius as (a), but
different water content. Without the effect of connectivity (in case of parallel
capillaries) we would find the non-hysteretic dotted line.
θ [-] water content at saturation
s
θ [-] residual water content
r
α [cm−1] empirical parameter
n [-] empirical parameter
The residual water content θ represents the small amount of water which is adsorbed
r
in thin films at low water potentials and which is barely mobile (see Fig. 3 for log(−ψ ) >
m
3 log(cmWC) in the sand). The parameter α scales h and thus determines the position
m
of the curve relative to the axis of matric potential and the parameter n determines the
shape of the curve. Thus, α is related to the mean pore size and n to the width of the
pore size distribution. However both parameters are fitting parameters without immediate
physical meaning. The parameters for the curves in Figure 3 are give in Tab. 1.
Hydraulic conductivity function
The conductivity K depends on the geometry of pores as well: large connections are better
conductors than narrow necks. Pores which are filled with air are not available for liquid
water transportation.
9
With decreasing water content, large pores are drained first and the remaining water
is located in smaller pores. Because the conductivity K depends on the radius of the
flow channel (j ∼ r4), it decreases rapidly with decreasing water content. As additional
w
effect, the fraction of pores which are used for water conduction decreases as well, and
the way between two points A and B which only consists of water filled pores gets longer.
Consequently, the conductivity of water is a function of the water content: K = K(θ ).
w
As θ is a function of ψ , K depends on ψ as well.
w m m
In Figure 3 the hydraulic conductivity functions K(ψ ) are shown for the sandy and the
m
loamy soil. Although the porosity of the sand is lower, the hydraulic conductivity at water
saturation is higher in the sand. This is because the pores are larger (see θ(ψ)). However,
for decreasing water potential the conductivity in the sand drops very rapidly. This is
intuitively clear, because at the critical water potential below h ≈ −30 cmWC also the
m
water content drops significantly. Hence, the shape of θ(ψ ) tells us something about the
m
shape of K(ψ ) which is also used for the parameterisation. Using arguments derived
m
from statistical pore models Mualem (1976) proposed to derive the hydraulic conductivity
curve through integration of the water characteristic
Θh (θ)−1dθ 2
K(Θ) = K Θτ 0 m . (9)
0 "R01hm(θ)−1dθ#
R
Inserting the van Genuchten parameterisation (8) this leads to
K(Θ) = K Θτ 1−[1−Θn/[n−1]]1−1/n 2 (10)
0
h i
where now the conductivity is given as a function of Θ. Note that the shape of the curve
is determined by the same shape parameter n as for the soil water characteristic. The
hydraulic conductivity at water saturation, K , is introduced as a parameter to fix the
0
absolute height of the hydraulic conductivity. Moreover an additional parameter, τ, is
used which is termed “tortuosity” and which is thought to account for the change in the
topology of the water phase with decreasing water content. However, τ is considered to
be a pure fitting parameter which frequently is fixed at a value of 0.5.
At this point we arrive at the complete set of hydraulic parameters, describing the hy-
draulic properties of porous media. These parameters are listed in Tab. 1 for the examples
shown in Figure 3. The aim of the first part of this course is to experimentally determine
such a parameter set for a given soil.
2.3.2 Experimental determination of hydraulic properties by inverse modelling
Direct measurements of hydraulic properties are very time consuming and, in case of the
hydraulic conductivity, also very difficult. In this practical course, we apply an indirect
method, the approach of inverse modelling, which is based on the theory given in the
previous chapter.
The basic idea is that we can describe the dynamics of the system through Richards’
equation (7). Thus, for a given set of initial and boundary conditions, we can predict the
behaviour of the system in case the hydraulic properties areknown — which is not the case
here. But we can observe the behaviour of the system and choose the hydraulic properties
10