## 1. Introduction

### 1.1. Classification of solutes, pollutants and subsurface pollution

Solute transport is of importance in view of the movement of nutrient elements, e.g. towards the plant root system, and because of a broad range of pollutants. Pollution of the subsurface is often considered to be either point source pollution or diffuse source pollution. Point source pollution covers a limited area, and is often caused by accidental (or illegal) spills (e.g. leaking pipes, tanks, mine tailings, etc.). Diffuse source pollution covers a large area and is in general caused by large-scale application of both beneficial and hazardous compounds at the soil surface (manure and fertilizer, pesticides, atmospheric deposition of acids and radio nuclides, etc.). Pollution is not necessarily man induced, but may be due to geological or geohydrological causes, e.g. in the cases of pollution with arsenic, and salt.

For the polluting species, a distinction can be made between dissolved and immiscible, and between conservative and reactive. Dissolved pollutants (aqueous phase pollutants) will spread with the groundwater due to groundwater flow, diffusion and dispersion. Immiscible pollutants will spread as a separate phase (non-aqueous phase liquids, NAPL). They will contain components with very low solubility in the water phase. They constitute a long-term source for pollution.

Conservative pollutants are those that do not react with the solid soil material, do not react with other pollutants and will not be degraded by biological activity. Reactive solutes may enter or leave the water phase through adsorption/desorption, chemical reactions, dissolution/precipitation and/or biodegradation.

### 1.2. Some basic definitions

*Advection* : the spreading of a pollutant by groundwater flow.

*Diffusion* : the spreading of a species dissolved in the water phase by the Brownian motion of the ions (molecules).

*Dispersion* : the spreading of a species dissolved in the water phase by local variations in the water velocity.

*Adsorption/desorption* : interaction of species dissolved in the water phase with the solid matrix. This process can be physically based or chemically based, reversible or irreversible.

*Chemical reactions* : reactions of species dissolved in the water phase with other species, resulting in the occurrence of different species altogether.

*Biodegradation* : the degradation of species dissolved in the water phase by bacteria.

*Radioactive decay* : the degradation of species by radioactivity.

Concentrations of species in the water phase C_{i} (including pure water itself) are defined as the mass of the species per unit volume: *kg/m*^{3}*, g/l, mg/l,* etc.

The density of a multi-component fluid, consisting of N components, is then given as:

Mass fractions *ω* of the components (mass per unit of mass: *kg/kg, g/g, etc.)* are defined as:

For dilute solutions (tracer concentrations) all mass fractions *ω*_{i} *<<1*, except for the pure water. This means that the density of the fluid is close to the density of pure water, and can be assumed to be constant.

Water density is a function of pressure, temperature and composition. This last dependence is only important at high concentrations. E.g. in case of seawater intrusion, or in deep saline aquifers which are sometimes used to store waste or to produce energy. In these deep aquifers salt concentrations can be as high as 300 *g/l*, resulting in a water density of 1200 *g/l* (giving a salt mass fraction of 0.25). Water density fluctuations will also play a role in the subsurface storage of heat.

Water viscosity is a function of pressure, temperature and composition. This influences the hydraulic conductivity (see next section). The dependence on the temperature is by far the most important. Hence, this dependence must be taken into account in the analysis of subsurface storage of heat.

### 1.3. Groundwater flow

Groundwater flow is described by Darcy’s law. Darcy’s law is in principle the form of the momentum balance (Navier-Stokes equation), averaged over a large number of pores. It also follows from a balance of forces on water flowing through a porous medium.

Consider flow in the z-direction (see Figure 1). The net forces (positive upward) working on a body of water with dimensions Δx, Δy and Δz are:

Pressure forces: (p(z)-p(z+Δz))ΔxΔy

Gravity forces: -ρgΔxΔyΔz

Friction forces: -R’q_{z}ΔxΔyΔz

where *R’* is the resistance factor and *q*_{z} is the specific discharge (Darcy velocity) in the z-direction.

For the pressure forces, we can make the following approximation by using a first order Taylor series expansion:

Then setting up the force balance, we find:

where we have assumed that the resistance factor *R’* is proportional to the liquid viscosity *μ*. *κ*_{z} is the intrinsic permeability in the z-direction (L^{2}), which is assumed to be a property of the porous medium. The intrinsic permeability of a porous medium is largely determined by the pore sizes and shapes. A strong correlation between permeability and porosity exists. Similar expressions can be obtained for the flow in x and y direction:

Basic assumptions in this derivation are that the acceleration of the water can be neglected, and that the friction forces are linear dependent on the velocity.

The latter is not always true (especially at high water velocities, e.g. close to an abstraction or infiltration well), in which case Darcy’s law is not valid, but should be replaced by Forchheimer’s equation:

where *β* is again a property of the porous medium.

Define a piezometric head *h* as:

Basically, the piezometric head consists of a pressure head *p/ρg* and the vertical position *z* with respect to the reference level. It is the position of the top of the water column in an observation well with respect to the reference level (usually mean sea level). This is different from unsaturated flow, which is formulated in terms of the pressure head.

Substitution of equation (7) in equation (4), assuming that the density *ρ* is constant then gives Darcy’s law in terms of the groundwater head *h*:

and similar expressions can be obtained for *q*_{x} and *q*_{y}.

Consequently, if the density *ρ* and the viscosity *μ* are constant, we can define hydraulic conductivities as:

and the same for the y and z direction. This shows that the hydraulic conductivity *k* (L/T) is dependent on the fluid properties.

The groundwater flow equation follows from a mass balance for the complete water phase (including all dissolved species). Consider the element as depicted in Figure 2 with dimensions Δx,Δy and Δz. The net mass influx over a period *Δt* in the x-direction is given by:

A similar expression can be obtained for the net mass influx in the y and z directions. The change in the total mass in the element is given by:

where *n* is the porosity.

Equating the net mass influx and the change in mass gives the mass balance equation for the liquid phase:

For situations with varying fluid properties (salt water intrusion, storage of heat, etc.) this equation together with the pressure formulation of Darcy’s law (equations (4) and (5)) should be used. Note, that the flow equation in that case is non-linear. Note also, that in these cases, even though a piezometric head can be defined, it will not be the driving force for groundwater flow.

If the density of the liquid *ρ* and the porosity *n* is assumed to be dependent on the pressure *p* only, the time derivative in the mass balance equation can be written as:

where *S*_{s} is the specific storage. Combining this equation with the piezometric head formulation of Darcy’s law (equations (8) and (9)) and division by the (constant) density gives the well-known groundwater flow equation:

Note, that the average pore water velocity *v* is different from the specific discharge *q* :

*v=q/n.*

## 2. Simplified description of processes in reactive transport

### 2.1. General

Similar to the water balance, we can derive a general form for the mass balance of a dissolved component in groundwater. Assume that the mass fluxes in x, y and z-directions are given by *F*_{x}*, F*_{y} and *F*_{z} (M/L^{2}T) respectively (see Figure 3).

The net mass influx in the x-direction over a period *Δt* is then given by:

and similar expressions can be obtained for the net mass influx in the y and z-directions.

The change in mass of the component in the element over a period *Δt* is given by:

Due to the different processes occurring, mass of a component can be produced or lost in a period, e.g. because of adsorption/desorption, chemical reactions, decay, etc. The loss of mass due to these processes per unit volume and unit time will be indicated by *I* (M/L^{3}T). Note, that *I* can be either positive (loss of mass) or negative (gain of mass). Combining the different terms then gives the following general mass balance equation:

In the following, the mass fluxes and/or the mass production associated with the different processes playing a role will be given. For the time being, simplified (linear) expressions will be given, which will result in a mass balance equation in the form of the classical Advection-Dispersion (or Convection Dispersion) equation, CDE. Later, more complicated expressions will be covered.

### 2.2. Advection

Advection (or convection) is the transport of dissolved components by flowing groundwater. The mass transport per unit area of porous medium of a dissolved component by flowing groundwater is given by:

where *F*_{x} is the mass flux of the component in the x-direction (M/L^{2}T), *q*_{x} is the specific discharge of water (Darcy velocity) in the x-direction (L^{3}/L^{2}T) and *C* is the concentration of the component in the water phase (M/L^{3}). No mass is produced or lost, hence, *I=0*.

The underlying assumption is that the average velocity of the ions or molecules of the dissolved substance is the same as the average water velocity: if we move one liter of water over a certain distance, also all chemicals in that liter will have moved that distance. In most cases, this will be true, but there are exceptions. These exceptions occur e.g. when the molecules of the dissolved substance are very large (colloids, virus). If we consider the flow of water in a capillary, the water velocity *v* at a distance *r* from the centre is given by:

where *v*_{avg} is the average water velocity and *r*_{0} is the radius of the capillary. For large molecules, only part of the capillary is available for transport. That can be caused by either the size of the molecules or by the electrical charges on the surface. As a result, the average velocity of such particles in a capillary will exceed the average velocity of the water itself. If the radius of such particles is given by *r*_{c}, it can easily be inferred that the average velocity of the particles compared to the average water velocity of the water is given by:

(find the average water velocity by integrating the water velocity from 0 to *r*_{0}, and find the average particle velocity by integrating the water velocity from 0 to *r*_{0}*-r*_{c}). For instance, for particles with a size of 20% of the capillary diameter, the average velocity is some 40% larger than the water velocity. These effects have been observed in virus and colloid transport.

If only advective transport takes place, the mass balance for a component follows from (17) and (18):

Now consider the 1-dimensional mass balance equation with constant porosity *n* and constant specific discharge *q*. Division by *n* then gives the following equation:

where *DC/Dt* is the material derivative, i.e. the change in concentration when moving along with a water particle. Since this derivative is zero, advective transport only results in a displacement of the initial concentration distribution by *vt*, where *t* is the elapsed time. This is also true in 3 dimensions.

Equation (22) can be written in dimensionless form by defining the following dimensionless variables:

where *C*_{r}*, t*_{r} and *L*_{r} are reference or characteristic values for the system considered. Substitution of these dimensionless variables in the mass balance equation (22) gives:

We can now choose any one of the characteristic values *t*_{r} or *L*_{r} such that the coefficient in front of the spatial derivative in (24) is 1. This means that for a given characteristic time *t*_{r}, the characteristic length is given by *vt*_{r}, while for a given characteristic length *L*_{r}, the characteristic time is given by *L*_{r}*/v*. These characteristic values obviously are related to respectively the travel distance and the travel time of water particles.

The specific discharge that is required to quantify the advective fluxes follows from the mass balance equation for the water phase in combination with Darcy’s law. Basically, this means that (local) information about the value of the permeability (or hydraulic conductivity) is required.

### 2.3. Diffusion

Diffusion is the spreading of a component dissolved in the water phase by the Brownian motion of the molecules/ions. In open water, the mass flux due to diffusion is given by Fick’s first law:

where *D*_{m} is the molecular diffusion coefficient (L^{2}/T), which is typical for the component considered. In a porous medium, the mass flux due to diffusion is given by a similar expression:

The porosity *n* enters the equation to account for the area that is effectively available for mass transport. *τ* is the tortuosity of the porous medium (-), which accounts for the fact that the length of the path molecules or ions have to take in a porous medium to travel from one position to another is larger than the distance between these positions. For normal porous media, *τ* has a value in the order of 1.6 to 1.7. No mass is produced or lost due to diffusion, hence *I=0*.

In case only diffusion occurs, the mass balance equation reads:

Now consider a 1-dimensional mass balance equation with constant porosity *n* and diffusion coefficient *D*_{eff}:

For this equation, numerous analytical solutions dependent on the boundary conditions are known, both in Cartesian and radial coordinate systems. The steady state solution in a Cartesian coordinate system is:

where *A* and *B* are determined by the boundary conditions. Note, that this solution is not dependent on the effective diffusion coefficient.

Equation (28) can be written in dimensionless for by defining appropriate characteristic values for the concentration and the time and length scales:

Setting the coefficient in front of the spatial derivative to 1, characteristic values for the time and the length are found. For a given characteristic time *t*_{r}, the characteristic length is given by *L*_{r}=*√*(*D*_{eff}*t*_{r}), and for a given characteristic length *L*_{r}, the characteristic time is given by *t*_{r}=*L*_{r}^{2}*/D*_{eff}.

In general, molecular diffusion will not play an important role in porous media transport, unless the groundwater velocities are very small (which is e.g. the case for transport through very low permeable clays).

Measurement of effective diffusion coefficients in a porous medium is usually done in the laboratory by performing time dependent experiments

### 2.4. Dispersion

Dispersion is the spreading of a dissolved component due to local variations in the groundwater velocity. In general, we distinguish mechanical and hydrodynamic dispersion.

Mechanical dispersion takes place on the pore scale, and is caused by velocity variations across the cross section of the capillaries (or pores). Usually the groundwater velocities are so small, as are the pore diameters, that molecular diffusion is fast enough to balance concentration differences in the direction perpendicular to the flow (i.e. across the pores).

Hydrodynamic dispersion is the sum of molecular diffusion and mechanical dispersion. It usually occurs on a larger scale than a single pore, and is caused by all variations in the average groundwater velocity (i.e., averaged over a large number of pores) that we did not account for explicitly, including diffusion (Figure 4). Thus, if we consider layers with different values of the hydraulic conductivity (or permeability), this variation does not necessarily give rise to hydrodynamic dispersion. However, it is clear that such variation certainly may lead to variation of the rate of displacement of chemicals and to true mixing, if it combines with diffusion. This is commonly called macro or mega dispersion (see e.g. Dagan, 1987), and considered later in this chapter (section 7).

For both mechanic and hydrodynamic dispersion, the mass fluxes are assumed to be given by the following form of Fick’s first law:

where *D*_{xx}, *D*_{xy} and *D*_{xz} are elements of the dispersion tensor (L^{2}/T). Similar expressions are valid for the mass fluxes in y and z-direction. The dispersion tensor is symmetric, and consists of 6 different numbers, *D*_{xx}, *D*_{xy}*=D*_{yx}, *D*_{xz}*=D*_{zx}, *D*_{yy}, *D*_{yz}*=D*_{zy} and *D*_{zz}. The elements of the dispersion tensor are dependent on the groundwater velocity *v*, such that the dispersion coefficients in the direction of the flow and perpendicular to the flow are given by:

where *α*_{l} and *α*_{t} are the longitudinal and transversal dispersivities (L) respectively. These are assumed to be properties of the porous medium, and indicate the size of heterogeneities in the system that is not accounted for by variations in the (average) groundwater velocity. Because mass transport by hydrodynamic dispersion and by molecular diffusion is described by the same law, they are usually combined.

In a fully 3-dimensional system, with velocity components *v*_{x}, *v*_{y} and *v*_{z} respectively, the elements of the hydrodynamic dispersion (including molecular diffusion) are given by:

Because hydrodynamic dispersion occurs only in combination with groundwater flow, a mass balance for a component follows from the combination of mass fluxes as defined by equations (18) and (31), with the elements of the dispersion tensor given by equation (33). Using a short hand notation, this mass balance is given by:

Now, consider a one-dimensional system with constant porosity, velocity and dispersion coefficient:

Making this equation dimensionless by choosing appropriate characteristic values for the concentration, time and length gives:

Choosing the characteristic time *t*_{r}*=L*_{r}*/v* then gives the following dimensionless mass balance equation:

where *Pe* is the Peclet number. This number is characteristic for the ratio of advective transport and dispersive transport. The solution of (57) depends only on *Pe*, and large Peclet numbers indicate that advection dominates; small Peclet numbers indicate that dispersion dominates.

For mechanical dispersion, the longitudinal dispersivity *α*_{l} is in the order of the pore sizes (mm scale). For hydrodynamic dispersion, the longitudinal dispersivity is dependent on the scale of the problem. For laboratory experiments in columns, values of less than 1 mm to values larger than 1 cm have been reported. For field scale experiments values larger than 10 m have been reported, dependent on the size of the experiment and the heterogeneity of the aquifer in which the experiment was performed.

A cautioning remark is needed with regard to such large dispersivities, as such large values cannot possibly be due to complete true mixing of water of different composition in porous media. Such large values are commonly obtained based on several methodological complications: (i) the equipment used to measure ‘local’ concentrations (e.g. observation wells, geophysical methods) are themselves responsible for mixing, and (ii) modelling with large spatiotemporal discretization in view of computational efficiency may lead to numerical mixing, and (iii) dispersivities may be ‘fitted’ using a relatively simple transport equation (e.g. a one dimensional version of the transport equation), which leads to artifacts.

Based on field experiments, an empirical relation for the longitudinal dispersivity has proposed as:

where *L* is a characteristic length for the domain of interest. For large scales (large values of *L*), an upper bound for *α*_{l} is reached. Note, that relation (38) is based on the evaluation of a number of field experiments, and gives an estimate of the dispersivities only. Also note, that these estimates are based on the assumption that the aquifer permeability is homogeneous over the domain of interest.

As a guideline, the transversal dispersivity is commonly assumed to be 5 to 10% of the longitudinal dispersivity.

### 2.5. Adsorption/desorption

Adsorption/desorption creates a sink/source term for a component in the water phase. Two processes take place at the same time: molecules/ions will attach to the solid material, and attached molecules/ions can be released from the solid to the water phase. For the time being, we will adopt a linear description of the process, corresponding with case A of Figure 5:

where *k*_{a} and *k*_{d} are the attachment (1/T) and detachment (M/L^{3}T) coefficients respectively, and *C*_{s} is the concentration of the component adsorbed (M/M). *k*_{a} and *k*_{d} have different units because of the different units for *C* and *C*_{s} respectively. In this formulation, no equilibrium has been assumed. *K*_{d} is the distribution coefficient. In case of equilibrium, the expression in (39) between brackets is 0, hence:

which defines a linear adsorption isotherm. In principle, such a relation is only valid if the concentrations are very low, and if equilibrium between the water phase and the solid material exists.

If advection, dispersion and linear adsorption/desorption occur, the mass balance equation for a component in the water phase can be given as:

Note, that this equation is not based on the assumption of local equilibrium. Also, two unknowns are present in this equation: *C* and *C*_{s}. The other equation required to solve for the concentrations is given by a mass balance equation for the adsorbed component:

where *ρ*_{s} is the density of the solid material. Basically, this mass balance equation is comparable to equation (41), because the adsorbed component is not transported (no advection and dispersion), and that the source term due to adsorption/desorption has the opposite sign.

Adding equations (41) and (42) gives the total mass balance:

If we now assume equilibrium, equation (40) can be used to eliminate *C*_{s} from equation (43):

where *R* is the retardation factor. Note, that for non-reactive solutes (no adsorption/desorption) *R*=1. It is clear from this equation, that the retardation factor is only found in the time-derivative term. For this reason, it got its name, as this factor *R* implies that both ad/convection and dispersion are *R* times slower: they are retarded by a factor *R*.

Now consider a one-dimensional form of this mass balance, with constant porosity, velocity, dispersion coefficient and retardation factor:

Note that this equation is identical to equation (35), the mass balance equation for a non-reactive component, with both the velocity *v* and the dispersion *D* scaled by a factor *R*. If the molecular diffusion can be neglected, the dispersion coefficient is proportional to the velocity *v*, and equation (45) will give the same results as the mass balance equation for a non-reactive component with a velocity that is decreased by a factor *R*.

Now consider the one-dimensional form of equation (41) (non-equilibrium) with constant porosity, velocity, dispersion coefficient and attachment and detachment constants. This equation can be made dimensionless by choosing appropriate values for the characteristic concentration, time and length. If (as done before) we define the characteristic time *t*_{r}*=L*_{r}*/v*, the non-dimensional equation is given by:

The last two coefficients in this equation are two forms of the dimensionless Damkohler number. This number gives the ratio of the groundwater travel time and the time required to reach equilibrium. Large Damkohler numbers indicate that the assumption of local equilibrium is appropriate, while small Damkohler numbers indicate that adsorption/desorption should be described as a non-equilibrium process.

Measurement of the adsorption distribution coefficient *K*_{d} is commonly done in a laboratory batch experiment. A soil sample is mixed with water that contains a dissolved component at a certain concentration. This mixture is stirred gently for a long time in order to assure that equilibrium between the water phase and the soil is establihed. From a measurement of the resulting concentration in the water phase at equilibrium, the amount adsorbed can be determined, and the distribution coefficient calculated:

Assume we have a mass of soil *M*_{s}. We add a volume of water *V* which has dissolved in it a component at concentration *C*_{i}. At equilibrium, the concentration in the water phase is measured as *C*_{eq}. The amount of mass of the component added to the system is *VC*_{i}. At equilibrium, the total mass of the component in the water phase is *VC*_{eq}. Consequently, the total mass adsorbed is *V(C*_{i}*-C*_{eq}*)*, and the concentration of adsorbed component is *C*_{s}*=V(C*_{i}*-C*_{eq}*)/M*_{s} while the concentration in the water phase is *C*_{eq}. The value of the distribution coefficient follows directly from equation (40).

Measurement of the attachment and detachment coefficients can be done in batch experiments by measuring the concentration in the water phase as a function of time.

Interaction of solutes do not only take place with the solid material, but can also exist with colloïdal particles (e.g. natural organic material), which in itself are mobile. Consequently, a competition between adsorption on the solid matrix and adsorption on colloïds may occur. This may lead to an enhanced transport of species (e.g. heavy metals) that may otherwise be considered to be highly retarded.

### 2.6. Decay

For the simplified description given in this chapter, we will assume that the decay due to chemical reactions, biological activity and/or radioactivity is given by a first order expression:

where *λ* is the decay/degradation constant. *λ* is related to the half life *t*_{1/2} of the component by:

The halflife *t*_{1/2} is commonly measured in batch experiments by mixing a sample of soil material with water which has the component dissolved in it. Measuring the concentration in the water phase as a function of time will give an estimate of the decay. In these experiments, adsorption/desorption should also be taken into account. The process of first order decay/degradation is of great importance for much of the transport theory. As may be already apparent, radionuclides decay proportional to the total decaying mass present. For instance the recent tsunami accident with the Fukushima nuclear plant in Japan may have resulted in soil contamination with radionuclides, where the decay rate determines the period for which radiation problems may be acute. Likewise, the Chernobyl melt down resulted in continental scale contamination with radionuclides by different elements, that move towards groundwater with different rates, and different degradation rates. To appreciate the hazard for life, the rate of downward movement of chemicals in relation with the decay rate of hazardous radiation is a typical transport problem.

The first order degradation rate law is also the most commonly used rate law for describing the degradation of contaminants such as pesticides, nutrient chemicals such as nitrate, contaminants such as PAHs, BTEX, chlorinated hydrocarbons (the last under anaerobic conditions), and other contaminants (Keijzer et al., 1999, Jaesche et al., 2006, French et al., 2009), despite that it ignores that transformation products may be hazardous too.

The importance of degradation can be appreciated from an example of groundwater (rather than soil) contamination. About, say, a decade ago, the concept of natural attenuation has been developed. This concept proposes that the subsoil environment is able to cause natural degradation of contaminants, e.g. due to the intrinsic activity of microbial populations. Although dispersional mixing and dilution, as well as volatilization of chemicals may contribute to natural attenuation, degradation is a major process in this concept. The concept as such is important as it diminishes the environmental hazards of soil contamination, and therefore has become a major issue in soil and groundwater contamination strategies, management, and decision making.

### 2.7. Full simplified mass balance equation

A full mass balance equation assuming linear equilibrium adsorption and first order decay can now be written as:

where it has also been assumed that the component dissolved in the water phase as well as the component adsorbed onto the solid phase can both decay with the same decay constant. In reality, this is not necessarily true (Beltman et al., 2008).

This equation is the foundation of most software aimed at modelling soil and groundwater contamination, such as MODFLOW/MT3D and related models. As such, this equation is the core of much scientific as well as management supporting investigations done at international, national and local levels. Few, if any, predictions and prognoses are made on the fate of contaminants, that are not based on equation (49).

For some one-dimensional and simple two- or een three-dimensional problems, analytical solutions exist for equation (49). A number of these solutions have been programmed and are available at the web site:

www.cee.uiuc.edu/transport

With those solutions, it is easy to obtain an impression of the effects of the different processes and their parameters on the transport behaviour of dissolved chemicals.

## 3. Some effects in the numerical solution of the transport equation

### 3.1. General

In many cases, analytical solutions of the simplified transport equation are not available. That is e.g. the case for heterogeneous systems or complicated boundary conditions. In those cases one has to resort to the numerical solution of the groundwater flow equation and the solute transport equation.

In order to analyse the behaviour of the numerical solution of the partial differential equation describing the transport of a solute, we consider a simplified 1-dimensional system with constant porosity *n*, constant specific discharge *q* and constant dispersion coefficient *D*:

For the evaluation of the spatial derivatives in a finite difference approach, we use the following Taylor series expansion:

where *Δx* is the blocksize in the x-direction.

The first order derivative of *C* with *x* can now be approximated in two ways. A backward finite difference approximation follows from the second equation given by (51):

Another approximation for the first order derivative can be obtained by taking the difference of the two equations given in (51):

An approximation for the second order spatial derivative is obtained by adding the two equations (51). After some manipulation, the following approximation is obtained:

where a higher order term not given in equation (51) has been taken into account. The truncation error for the approximation of the first order derivative is for the backward difference (52) of the order *Δx*, and for the central difference of the order *(Δx)*^{2}. In other words, the central difference approximation is more accurate. The truncation error for the approximation of the second order derivative is of the order *(Δx)*^{2}.

For the time derivative the following approximation can be obtained:

where *Δt* is the time step and the derivative in the higher order term is evaluated at time *t+Δt* (i.e. at the end of the time step). Note, that the truncation errors in the approximation for the spatial and temporal derivatives indicate, that small grid blocks (small *Δx*) should be used where the second order derivative of the concentration with respect to *x* is large, and small time steps should be used when the second order derivative with respect to *t* is large.

In the following, we will adopt the notation: *C*_{i}=*C(x)*, *C*_{i-1}*=C(x-Δx)*, *C*_{i+1}*=C(x+Δx)*, while values evaluated at the beginning of a time step will have a superscript *o*, and values evaluated at the end of a time step will have a superscript *n*.

If we evaluate the spatial derivatives at a time level between the beginning of the time step and the end of the time step, the discretised mass balance equation can then be written as:

where *θ* is a factor between 0 and 1. For *θ*=1, all spatial derivatives are evaluated at the new time level (the end of the time step). This is a fully implicit scheme. For *θ*=0, all spatial derivatives are evaluated at the old time level (the beginning of a timestep). This is a fully explicit scheme. A mixed scheme (known as the Crank-Nicholsen scheme) is obtained by setting *θ*=0.5.

Note that in equation (56) a backward difference for the advective term is used. A similar expression can be obtained for a central difference approximation for the advective term.

### 3.2. Numerical dispersion

Numerical dispersion is an extra dispersion in the numerical solution of the transport equation which is caused by the discretisation of the advective term.

Consider the discretised equation (56), and evaluate all spatial derivatives at the end of the time step (fully implicit, *θ=*1):

One could now ask the question which partial differential equation is approximately solved by these equations if higher order terms are taken into account. Using equations (52), (54) and (55), keeping all terms with derivatives of order 2 then gives the following partial differential equation:

In order to obtain an expression for the second order derivative with respect to time, we will differentiate equation (50) with respect to time:

An expression for the cross derivative term is obtained by differentiating equation (50) with respect to *x*:

Substitution of equation (60) in (59) then gives:

Substitution of (61) in (58), collecting the second order derivative terms and neglecting higher order terms then gives:

Equation (62) is up to second order terms identical to the discretised equation (57). In other words, the discretised equation is an approximation to a solute transport equation with enhanced dispersion (cf. the term between brackets). Note, that in the analysis higher order terms are neglected. As a consequence, the discretised equation (57) is not completely identical to equation (62).

The extra dispersion *vΔx/2+v*^{2}*Δt/2* is called the numerical dispersion. If we assume that molecular diffusion can be neglected, the physical dipersion is given by *D=αv*, where *α* is the dispersivity of the medium. The numerical dispersion can now be neglected if the block size *Δx* and the time step *Δt* are chosen such that:

Carrying out the same analysis with the advective term in equation (57) approximated by a central difference, will show that the numerical dispersion in that case is *v*^{2}*Δt/2*, which is smaller than the one for the backward difference of the advective term. There can be, however, reasons for adopting the backward difference approximation (see next section).

There are a number of ways in which we can neutralize the effect of numerical dispersion. The most obvious way is to correct the dispersivity for the numerical dispersion. Suppose the physical dispersivity is given by *α*_{f}, while the dispersivity defined for the numerical calculations is given by *α*_{m}. For a finite difference approximation with a backward difference for the advective term, the following choice:

will result in a total dispersion in the numerical calculations equal to the physical dispersion. This approach can, however, only be adopted if the model dispersivity *α*_{m} remains positive because negative values of the model dispersivity will generate instabilities in the numerical solution.

More complicated ways to minimise numerical dispersion in the numerical solution of the solute transport equation can be thought of. In all cases one should consider the fact that numerical dispersion is solely caused by the first order spatial derivative (the advective term).

One way of avoiding numerical dispersion is by what is called operator splitting. In that approach the change in concentration is split in two parts: one part due to advective transport, and one part due to dispersion:

and adding the two contribution gives the full transport. The first part in equation (65), the advective part is now solved by a characteristic method. In this method, water particles are followed as they are transported with velocity *v*. Each water particle represents a certain mass of solute, which in a time step *Δt* is transported over a distance *vΔt*. Once the advective transport has been solved, the dispersion has to be added. That again can be done in a number of ways, the most simple being a finite difference approximation. Another way is the random walk method, where the effect of dispersion is simulated by random displacements of the water particles around the mean displacement given by the velocity *v*. These random displacements are related to the dispersion coefficient.

Using a characteristic method to simulate the advective transport has the disadvantage that only a discrete number of particles can be followed, and that interpolation is required to transform mass per particle to concentration distribution. A smooth concentration distribution from the distribution of the particles can only be obtained if a very large number of particles is used. This is especially true for the relative low concentration contours.

A method that strongly resembles the characteristic method is the Eulerian-Lagrangian method. In this method, the time derivative is approximated by taking the difference in concentration not at the same place, but at different places. Using a Taylor series expansion, it can be shown that the expression:

is second order correct, i.e. neglected derivatives are of third order or higher. Consequently, no numerical dispersion is generated by this method. In a standard finite difference method the second term in the left side of (66) is evaluated by evaluating the concentration at position *x-vΔt* at the beginning of a time step. That can usually be one by interpolation, although special precautions have to be taken close to boundaries.

Each of the methods mentioned here can easily be extended to two or three dimensions.

### 3.3. Oscillations in the solution

Some of the finite difference schemes can, under certain conditions, generate oscillations in the solution, i.e. negative concentrations or concentrations larger than the maximum value defined by the initial and boundary conditions may occur. It should be pointed out that these are not instabilities, where errors in the solution can grow unbounded.

Consider the discretised, implicit equation, where the advective term is approximated by a backward difference (equation (57)). With some manipulation, this equation can be written as:

Equation (67) shows that the concentration at the new time level is a weighted average of the concentration at the old time level and the concentrations in the adjacent grid blocks. Inspection of the weighing coefficients *β* shows that:

As a consequence *C* obeys a maximum principle, i.e. it can never become smaller that the smallest value given in the initial and boundary conditions, or larger than the largest value given in the initial and boundary conditions.

If we would have used a central difference for the advective term an expression similar to equation (67) can be written, however, with weighting coefficients:

Inspection now reveals that the sum of the weighting factors *β* again equals 1, but that all weighting factors are positive only under the condition:

Condition (70) is often given in a slightly different form:

where *Pe*_{cell} is called the cell Peclet number (cf. equation (37)).

Basically, a backward difference for the advective term generates more numerical dispersion than a central difference. However, a central difference approximation might result in oscillations in the solution, which will not occur for a backward difference.

It should be pointed out that the global mass balance will in all cases by perfect, irrespective of numerical dispersion or oscillations. This also indicates that damping oscillations in a solution by simply not allowing the concentration to become larger than a predefined value or smaller than another predefined value will ultimately result in mass balance errors.

### 3.4. Stability of the explicit solution

The explicit formulation of the discretised equations has the advantage that an explicit expression for the concentrations in each grid block is obtained, which therefore does not require matrix manipulation to obtain the solution. However, under certain conditions, such an explicit formulation may become unstable, i.e. small errors in the solution may grow uncontrolled and unbounded in time, resulting in very large positive and negative concentrations.

Appendix A gives the derivation of the criteria for the explicit formulation to be stable. In general, conditions (139) and (140) are given in a slightly different form. Condition (139) can be written as:

Substitution of this relation in condition (140) then gives:

which can then be given as the well known Courant condition:

and the Neumann condition:

Both conditions have a physical interpretation. The Courant condition states that in one time step, a water particle cannot travel further than the length of a grid block. The Neumann condition relates the characteristic length associated with the dispersion over a time step to the block size.

For a stable solution it is required that both conditions are satisfied, which means that the most restrictive condition determines the time step size that will still result in a stable solution.

Even an instable solution will give a perfect global mass balance (provided we are able to calculate the mass balance with enough significant digits). Consequently, a perfect mass balance (very small errors) is no guarantee for a good solution. However, large errors in the global mass balance is a guarantee for errors in the concentration distribution.

## 4. Initial and boundary conditions

In order to be able to model the transport of (reactive) solutes in groundwater it is necessary to define both the initial and boundary conditions. Initial conditions are (mathematically) only required for transient or time-dependent problems.

For local pollution problems, the initial conditions (a clean soil) can usually be given for the time before the pollution or spill occurred. For diffuse sources of pollution, the initial condition (present day situation in cases predictions have to be made) are derived by simulating long periods before the present day, using known (or estimated) mass inflow.

Boundary conditions for solute transport can be defined similar to the boundary conditions for groundwater flow. We distinguish three types of boundary conditions:

**Dirichlet:** the concentration on the boundary is fixed. Although mathematically this is a valid boundary condition, it is physically almost always impossible to create such a boundary condition. Nevertheless it is often applied at inflow boundaries, where the concentration of the solute in the water phase is known.

**Neumann:** the total mass flow of a solute across a boundary is defined. This type of boundary condition (also called mass loading) is often applied to inflow boundaries in e.g. experiments, but also in field situations dealing with point sources of pollution where the total mass of solute entering the groundwater is known (or can be estimated). For outflow boundaries, this type of boundary condition is physically not possible, because at outflow boundaries the total mass leaving the system is unknown (dependent on the concentration in the aquifer).

**Cauchy:** the mass flux of solute across the boundary is dependent on the concentration in the water phase itself. This type of boundary condition is often applied at outflow boundaries, defining the total mass flux as *qC*, where q is the water flux and *C* is the concentration of the solute in the water. Implicit in this definition of the boundary condition is the assumption that the dispersive flux across the boundary is zero.

For local point source pollution problems, we choose boundaries far enough away from the point source to make sure that the boundary conditions do not influence the concentration distribution.

## 5. Non-linear, non-equilibrium processes

### 5.1. General

Diffusion and radio-active decay have already been described in section 2, and do not need further elaboration.

### 5.2. Advection and dispersion

In heterogeneous systems, advection and dispersion are closely related. Dispersion is used to describe the transport of a contaminant due to variations in the groundwater flow which are not described by our “model” (where model can mean anything from complex numerical systems to the assumption of uniform flow). These variations are caused by local heterogeneities which have not been taken into account.

Basically, this means that the more information on convective or advective transport (in fact the variation in hydraulic conductivity) are taken into account, the smaller the dispersion will be.

From large numbers of field measurements, it is known that within one geological formation, the hydraulic conductivity may show a log-normal distribution with for instance an exponential covariance function defining the spatial correlation:

where *Y=ln(k)*, *σ*_{Y} is the standard deviation in *Y*, *r* is the distance between the points where *Y*_{1} and *Y*_{2} were measured, and *l*_{Y} is the correlation length. The basic assumption is that the co-variance is dependent on the distance between the measurement points only. This relation can be extended to account for direction dependence. For instance, the fact that sedimentation in geological formation have taken place in a certain direction will generate this direction dependence (correlation length in the *z*-direction will be smaller than the correlation length in the *x* and *y*-directions). Further, more practical aspects are discussed later in this chapter.

If we consider the formation to be homogeneous, all heterogeneities in the formation will have to be described by dispersion. In such a system, the dispersivity will be related to the correlation length (relation dependent on the flow pattern). That is, however, only true if the plume “has seen” all heterogeneities, i.e. if the size of the plume is (much) larger than the correlation length. In that case we will call the plume “ergodic”. For small plumes (non-ergodic), that is not the case, and the behaviour of such a plume cannot be described on a large scale by global dispersion. For the behaviour of such plumes it is necessary to incorporate information on local heterogeneities.

Note, that hydrodynamic dispersion generates a spreading in the average values of the concentrations, where the averaging volume is determined by the scale on which we assume the system to be homogeneous. If that scale is large, we are in principle not allowed to make a comparison of calculated concentrations with local measurements. Another reason to use spatial moments.

### 5.3. Adsorption/desorption

Equilibrium adsorption/desorption, especially for larger concentrations, can often be described by either a Freundlich equation (Figure 5, cases B and C)

where *K* is a constant, and *p*<1, or by a Langmuir isotherm:

where *C*_{smax} the maximum amount is that can be adsorbed (occurs when *kC*>>1). For low concentrations (where *kC*<<1) relation (78) becomes linear. The Langmuir isotherm is typical for soils and solid surfaces that have a limited number of sites available for adsorption.

For the non-linear adsorption, we can still define a retardation factor:

which will, however, be dependent on the concentration *C*. The mass balance equation is then given by:

where *F* is the mass flux (by advection and dispersion).

Inspection of the Freundlich and Langmuir isotherms shows that the retardation factor becomes smaller for larger concentrations. For the displacement of a front of solute, this means that the higher concentrations are less retarded than the lower concentrations. In other words, the higher concentrations try to “overtake” the lower concentrations. This effect is counteracted by dispersion. After some time, equilibrium will occur between these competing processes and a travelling wave will develop (Bosma and Van der Zee, 1993). An example of such displacement is shown in Figure 6.

In case of non-equilibrium, the adsorption/desorption may also be non-linear due to the limited number of sites available. For non-equilibrium adsorption/desorption we also need to solve for a mass balance of the adsorbed solute. For instance, a Langmuir type non-equilibrium interaction takes the form:

where *θ* is the fraction of the adsorption sites that is occupied (Van der Zee et al., 1987). Non-equilibrium often results in tailing in breakthrough curves.

If more solutes are present, competition for the adsorption sites may occur. One of the effects that may occur, is, that some solute adsorbs fast, but will later be replaced by a solute that adsorbs slower, but has a higher affinity for the adsorption sites.

### 5.4. Chemical reactions

Chemical reactions are usually described by equilibrium reactions. These are in general highly non-linear, while a large number of species play a role. However, they are local in nature (there is no spatial partial derivative in these equations). Solving for chemical reactions in combination with transport can be done in two ways:

Operator splitting: solve for the transport of all species required. The result is a redistribution of the concentrations. Starting with these concentrations, calculate the new chemical equilibrium, etc. Advantage is that the number of equations is limited. Disadvantage is that the time step size is limited.

Combine the equations for chemical equilibrium with the transport equations, and solve simultaneously. Advantage is that the time step size is less limited than in the operator splitting method. Disadvantage is the large system of non-linear equations that need to be solved.

Note, that negative concentrations (however small they are) cannot be allowed when dealing with chemical reactions.

Dissolution/precipitation reactions play a slightly different role. These reactions are dependent on threshold values. Consider a mineral the dissolves in the water phase. As long as the mineral is present, the concentration in the water phase is constant, and the mass transfer from the solid to the liquid phase is unknown. When the mineral is not present, it will act in the water phase as any other species. Until the concentration becomes large, and precipitation starts to occur.

One of the biggest problems in chemical reactions is to limit the number of species (and hence the number of mass balances that need to be considered) that have to be taken into account (Schroeder, 2005).

### 5.5. Biodegradation

Biodegradation needs at least, beside the present of bacteria, also the presence of a carbon source and the presence of an electron acceptor (if degradation is aerobic). The degradation can also be anaerobic (e.g. of inflammable NAPLs such as chlorinated hydrocarbons), in which a chemical is needed that supplies electrons for the degradation process. In general, aerobic degradation is faster than anaerobic degradation.

Biodegradation can be assumed to take place in the water phase only. Dependent on the number of bacteria present, we can distinguish between:

single bacteria in the water phase; these can also be transported with the flowing groundwater;

colonies on the solid matrix; these bacteria are not transported, but they still have access to all species in the water phase;

biofilms on the solid matrix; for bacteria in biofilms, mass transport from the free water phase to the biofilm is usually diffusion controlled, and this (slow) mass transfer has to be taken into account when describing biodegradation.

If a carbon source is available, the biomass (number of bacteria) will grow, and in time we can have a transition from free bacteria to colonies to biofilm. Two other effects play a role in biodegradation: 1) bacteria die, which can usually be described by a first order decay, and 2) species may be present that inhibit the biodegradation.

If a first order description of biodegradation is not sufficient, a typical way to describe biodegradation then is e.g. by Monod type relations:

where *I*_{max} is the maximum amount that can be degraded, *C* is the concentration of the carbon source, *O* the concentration of the electron acceptor, *I* the concentration of the inhibitor, and *k* are constants. Equation (82) can be extended to include more species that play a role in the biodegradation. Note, that for each of these species a mass balance needs to be solved. Equation (82) does not include the biomass, but that can be done in the same way, provided a mass balance for the biomass is solved as well.

If biodegradation takes place, the products of the degradation may in itself be degraded again by the same (type of) bacteria. In these cases it is required to analyse (or predict) the transport of multiple solutes, where the degradation of one solute results in a source term for the daughter product. These multiple solutes will in general have different adsorption/desorption characteristics and different decay characteristics.

### 5.6. Non Aqueous Phase Liquids (NAPL)

Non-aqueous phase liquids or oil does not mix with water. Infiltration of these liquids in the subsurface (e.g. through leaking tanks or pipes) will create multiphase systems.

If a spill occurs, oil will be transported through the unsaturated zone, leaving behind a residual oil saturation (around 20-30% of the pore volume), which is not mobile due to capillary forces.

When the oil reaches the water table the oil will float on the water table if it is lighter than water (LNAPL). Such floating lenses have been found on many places worldwide, often due to spills of gasoline at gasoline stations. The analysis of such lenses has been studied experimentally (Wipfler et al., 2004) but also numerically and analytically (Van Dijke and Van der Zee, 1998). If the NAPL is denser than water (DNAPL), it will move downwards through the groundwater until it reaches an impermeable layer. The transport of DNAPL through the saturated zone is controlled by instabilities (local heterogeneities), and it is therefore very difficult to predict where exactly the DNAPL will be present. However, residual oil will be present at those locations where the DNAPL has passed through.

Components of the oil will dissolve in the water phase, although generally in very small quantities. If the oil phase is in equilibrium with the water phase, the concentration of the oil components in the water phase can be described by Raoult’s law:

where *C*_{max} the maximum concentration is if the water phase is in contact with the pure component, and *X* is the mol fraction of the component in the oil phase.

For a layer of the pure oil, there is in general no equilibrium between the water phase and the oil phase. In these cases, the mass transfer is diffusion controlled, and described by a first order relation:

where *k* is the mass transfer coefficient, and *C*_{eq} is the equilibrium concentration, which is given by equation (83). The mass transfer coefficient is dependent on a.o. the diffusion coefficient, the oil saturation, average pore diameter and the groundwater velocity. A number of empirical relations exist that define these relations.

If a free oil phase is present, the transport of oil components requires both the solution of the oil and water flow equations and the transport equation of each species in both phases. Note, that the oil and water flow are coupled due to the capillary forces, and the dependence of the water hydraulic conductivity on the oil saturation.

## 6. Density dependent flow and transport

### 6.1. General

In many cases, the properties of groundwater, and in particular the density, are influenced by the concentration of dissolved species. That is e.g. the case for sea water intrusion in coastal aquifers, or for the infiltration of leachate from a landfill into a fresh water aquifer. In many of these cases, we may assume that the density is the only property of the water phase that is influenced by the concentrations of the dissolved species. In other cases, we have to take into account the fact that the viscosity of the water phase is also influenced by the concentration of the dissolved species. That is in particular true for very high solute concentrations (or for large temperature differences). For instance, if we consider the disposal of waste in deep saline aquifers, the effect of a changing viscosity cannot be neglected. The viscosity of water varies by a factor of 2 from fresh water to concentrated brine. These situations differ from those of NAPLs, as the different types of groundwater are miscible, whereas NAPL and water are immiscible.

In case of changing fluid properties, the governing flow equation cannot be formulated in terms of groundwater potential but have to be given in terms of pressures and a gravity term. In the following, the governing equations for density dependent flow and transport will be given, followed by a simplified set of equations.

### 6.2. Basic equations

For density dependent flow, the general form of Darcy’s law needs to be used:

The mass balance equation for the water phase is given by:

and the mass balance for the solute by:

where it has been assumed that there are no external sources and sinks, and that the solute creating the density differences is a conservative one.

We can now define a “fresh water potential” as:

where *ρ*_{f} is the density of fresh water (*ω=0*) and *z* is the vertical position with respect to a reference level. Substitution of (88) in (85) gives the following form of Darcy’s law:

which shows that the hydraulic conductivities are dependent on the fluid properties. It is also obvious that the fresh groundwater potential is not the only driving force for the groundwater flow. The system of equations (85) or (89) with (86) and (87) have to be supplemented by equations of state for both the density and the viscosity of the water.

The mass balance equations with Darcy’s law form a set of coupled, non-linear equations. The coupling is a two-way coupling: the flow is dependent on the concentration distribution (through the dependence of the density and the viscosity on the concentration), while the concentration distribution is dependent on the flow through the advective term in equation (87).

### 6.3. Simplified equations

In many cases, the equations governing the density dependent flow and transport can be simplified by making a number of assumptions. First of all, for relatively low concentrations (e.g. the salt concentration in seawater) we may assume that the changes in the viscosity with the concentration are negligible. Furthermore, the density can be assumed to be linear dependent on the salt mass fraction. That is certainly the case for salt water:

where for salt, *γ* has a value of 0.7.

Secondly, we will assume the (temporal) changes in the porosity can be neglected.

Finally, we can adopt Boussinesq’s approximation, which states that the variations in the liquid density can be neglected everywhere, with the exception in the gravity term of Darcy’s law. With all these assumptions, the governing equations can be written as:

Mass balance of the water phase:

Mass balance of the salt:

Darcy’s law:

Equation of state:

where *e*_{z} is a unit vector in the z-direction (positive upward).

Substitution of the equation of state (94) in Darcy’s law (93), and assuming that the porosity and the dispersion coefficients are constants, even further simplifies the set of equations:

These equations can be made dimensionless by choosing appropriate reference values for the different variables in the equations:

where *L* is a characteristic dimension of the system considered.

Now, if we choose the following reference values:

the system of equations reduces to:

where it is understood that all variables are dimensionless (subscript *d* has been omitted), and:

is the Rayleigh number. This number defines the ratio of gravity and dispersive forces. Together with the initial and boundary conditions it fully controls the solution of equations (100) through (102).

For the stable situation, where fresh water (derived from rainfall) is situated above saline groundwater, Eeman et al. (2011, 2012) provided analyses for the behaviour of such rainfall lenses in coastal delta areas. Experimental evidence of such lenses was given by Lebbe et al. (2008), Vandenbohede et al. (2008) and Eeman et al. (2011, 2012). Their work was preceded by over a century of investigations of the dynamics of fresh water pockets in coastal dunes (Herzberg, 1901, Van der Veer, 1977, Maas, 2007), barrier dunes, and atols (see Eeman et al. 2012).

For systems that are possibly instable, i.e. systems where fresh water is overlain by water with a higher density, the Rayleigh number controls whether such instabilities will occur. These situations arise e.g. under landfills, where the (heavier) leachate infiltrates in a fresh water aquifer, or at transgression of the coast line by rise of the seawater level.

For reasonable simple geometries, a perturbation analysis can show us when instabilities may occur. For the Rayleigh-Bernard problem, i.e. a rectangular vertical slab, with no-flow boundary conditions at the sides, and *ω*=1 at the top boundary and *ω*=0 at the bottom boundary, it can be shown that instabilities will occur for *A>4π*^{2}.

## 7. Transport in heterogeneous media

### 7.1. General considerations

If we wish to convey our understanding of transport of solutes to real soils, we have to account for spatiotemporal variability. After all, it is well known that soils vary spatially (layers, horizons, pores size distribution), and as a function of time, e.g. due to time varying weather, groundwater flow and many other causes. Such variability leads to rather complex behaviour in space and time. Instead of the rather simple, smooth concentration profiles and breakthrough curves, that we obtain for homogeneous soil and simple initial and boundary conditions, quite involved transport trajectories and concentration distributions result.

These are difficult to convey to the stakeholders, that have to base decisions on such results. For such reasons, spatial and temporal moment theory can be used to capture the transport phenomena in relatively robust terms. The moment theory is explained and for the case of spatial moments, it is presented mathematically. This is not done for temporal moments, because the mathematical details are completely in analogy to the spatial moments, and therefore obsolete here. Before going into details, first a qualitative impression is given of the moments.

The first spatial moment is related to the solute velocity, or, in case of a conservative solute, the groundwater velocity. For an instantaneous release of a non-reactive tracer in a steady state groundwater flow field, the first spatial moments as a function of time tell us exactly where the tracer is located and therefore also what the groundwater velocity is. Comparison of the first spatial moment of a reactive solute with the first spatial moment of a non-reactive tracer gives an estimate of the retardation factor. For a solute plume, these first moments characterize therefore the velocity of the entire plume, not of individual solute particles.

Following the zero’th spatial moment of a solute in time gives information about the degradation of the reactive solute. It is even better to compare these moments with the zero’th moment of a non-reactive (inert) tracer. If we indicate the properties of the reactive tracer with superscript *r* and the properties of the non-reactive tracer with *n*, we can define:

The decay constant *λ* of the reactive solute is then given by:

We could have used the zero’th spatial moment of the reactive solute only (as a function of time). However, comparing it with the zero’th spatial moment of a non-reactive solute may give a better answer because possible errors due to limited available information may cancel if the concentrations required for the determination of the zero’th spatial moment are measured at the same locations.

The second central spatial moment is related to the dispersion on the scale of the plume. For an instantaneous release of a solute in a steady state groundwater flow field, one of the elements of the dispersion coefficient is given by:

and similar expressions can be given for the other 5 elements of the dispersion tensor (cf equation (33)).

The accuracy with which the spatial moments can be determined is of course dependent on the amount of data available. Higher order spatial moments are often less accurate than the lower order spatial moments.

Characterizing the behaviour of a plume of contaminant on the basis of spatial moments is usually more robust than trying to characterize it by individual measurements of concentrations. The latter are very dependent on local heterogeneities. Apart from that, decision makers are usually interested in global measures like the transport of the plume (first spatial moment) and the spread around the mean travel distance (second central spatial moment) rather than in local concentrations.

### 7.2. Moment theory

As flow and transport conditions are spatiotemporally variable, complex patterns can develop, as illustrated in Figure 7 for the distribution of a chemical that moves through soil. This pattern is difficult to describe to someone who cannot see it. However, often that is not really important, because we are only interested in simpler information. Simpler information, that is also statistically robust, are the so-called moments. Moments can be described for any property that is distributed (in space, in time). For instance, if we consider the concentration of a solute as in Figure 7, that is distributed spatially, we can regard this distribution as a probability density function pdf (or: frequency function), in particular: a pdf of travelled distances, in this case. Its discrete analogue is a histogram, provided the lengths of all bars are normalized in such a way that their sum is equal to 1. A pdf *f* is related with a cumulative distribution (*P*) according to:

and for integration to *x*=+∞, *P* becomes one: it is the probability that the value of *x* lies between the lower and the upper boundary, and this probability is one for this *x*=+∞ upper boundary. Therefore, if we deal with a concentration distribution, we first have to make a pdf out of it. We do so with the zero’th moment.

Instead of doing so, we will illustrate moments on mass instead of concentration. A concentration is a mass (of solute) divided by a mass or volume of water. You can of course calculate what happens if you add one glass of water with a particular concentration to another glass with another concentration. However, physically it is much better to go to more basic properties: how much mass of solute is in each glass, how much water, and what is the final mass of solute and of water if both are combined? Hence: never calculate moments of concentrations, or water fractions. Instead, calculate moments of quantities (kg, volume,...). For the present purpose, we calculate moments of solute mass in solution (

The zero’th moment of the solute mass distribution in solution is given by:

if water fraction and concentration are only distributed in the x-direction (otherwise, we have to integrate for x, y, and z). The zero’th moment is also known as the mass of the distribution, and is a ‘normalizing’ property: division of higher order moments by the mass, renders the distribution of water fraction times concentration a pdf!

The first moment is given by:

and written this way, illustrates that division by the 0^{th} moment is needed to obtain a pdf. The first spatial moment divided by the 0^{th} spatial moment is known as the average position of the plume of solute (*x*_{av}), whereas the first temporal moment divided by the 0^{th} temporal moment is the mean breakthrough time. The second moment is obtained by multiplying with *x*^{2} instead of *x*. More common to use is the second central moment, given by:

from which the variance can be obtained as:

This variance is a good measure for the spreading in space (here the x-direction), and is related with the concept of dispersion, as mentioned before. To describe the transport behaviour of solutes, the moments up to the variance are not always sufficient. The main reason is that the transport velocities, i.e., including retardation effects, of parcels of solute may not be symmetrically (e.g. normally) distributed. In that case, higher order moments may be needed, or the assumption that the pdf has another than normal shape is required. An important alternative to the normal (symmetric) pdf is the lognormal pdf.

### 7.3. The lognormal pdf

The normal pdf of *X* is given by:

where *s*_{X} is the standard deviation and *m*_{X} is the mean. The lognormal pdf of *Y* is very similar, namely:

where the mean *m*_{Y} and standard deviation *s*_{Y} are those of the lognormally transformed parameter *X*=ln(*Y*). That means that if *Y* is lognormally distributed, then *X* is normally distributed, and the pdf of lognormally distributed *Y* is characterized with the statistics of *X*. The following equation shows that this is correct:

So, if we wish to characterize the lognormal pdf, we first have to log-transform, to calculate the statistics, and then we know the pdf. There is another way, as the statistics of both pfd’s are related. The following relationships can be derived if *Y* is lognormally distributed and if *X*=ln(*Y*):

So, to calculate the mean, you first have to calculate the variance (of *X*), from the statistics (mean and variance) of *Y*.

Median and mode for the normally distributed *X* are equal to the mean.

If you know the statistics of *X* (for instance log(hydraulic conductivity)), but not those of hydraulic conductivity itself (*Y*), then you can calculate those easily by inverting the above equations. The result is given by:

and for completeness, also the median and modus:

The lognormal pdf is very attractive for several reasons. An important reason is that many properties appear to be well lognormally distributed, whereas they are not normally distributed. For instance, the saturated hydraulic conductivity *K*_{s} has been found experimentally to be often lognormally distributed, and Miller and Miller (1955) also give a theoretical basis with the similar media or similitude theory.

### 7.4. Approaches to model transport in heterogeneous media

Modeling of heterogeneous media can be done in several ways, that vary in complexity. In all cases, an impression of the heterogeneity is gained experimentally. This gives us information regarding the statistics of the important soil properties, such as hydraulic conductivity. This information can be directly used as input for any model. However, it is also possible to generate a field of hydraulic conductivities and use that for input. This is attractive, because often we have data only for a very limited number of positions, leaving large volumes for which we do not have data and for which we have to find a good strategy to parameterize.

If we generate random fields, we can do a calculation for each generated field. If we do so repeatedly, for a designated set of statistics, then the calculations will be similar but different: just as when you take a picture of a meadow every five minutes: the meadow will be the same, but each picture is different. Repeated calculations for the same statistics is what we call Monte Carlo simulation. With this technique, we can determine the uncertainty that is intrinsic to our model results.

Monte Carlo simulation is a numerical approach, with which we can check whether analytical solutions are sound. But Monte Carlo simulation may also be done using analytical solutions, if they are available. Besides Monte Carlo, a common way to deal with heterogeneity is in the context of a GIS system, where each cell or pixel gets its own value, depending on overlay maps and transfer functions. Then, for each cell or pixel the same set of calculations is done, only with different parameters. An example of the output obtained this way is shown in Figure 8 for pesticide leaching. In these calculations, geospatial data of soil type, weather conditions, geohydrology, land use (which pesticides, which crops, which growing season and pesticide application date), and other information can be linked. Clearly, then it is possible to determine which regions have a major hazard of leaching a particular pesticide (or more generally, a contaminant of interest).

This type of modelling is very important for risk analysis and screening of pesticides in the admission of these chemicals. Whether or not a pesticide is admitted to the market, or should be taken out of the market (as happened about 1990 with quite a number of pesticides in the EU), may depend on the outcome of fate calculations with models as described here. For instance, the pesticide screening in EU uses the so-called FOCUS protocols, as mentioned by Beltman et al. (2008).

In a relatively simple approach, the spatial configuration is neglected. This is suitable, if the interest is not where something happens but what its effect is on the entire system. An example is leaching of solute from a field: instead of wishing to know where the leaching occurs, it may be sufficient to know what the average leaching is. This is the essence of the parallel stream tube PST model.

### 7.5. Parallel stream tube model

This approach has been developed by Bresler and Dagan in two highly innovative papers in 1979, and their analyses for a conservative tracer was later extended to reacting chemicals (Van der Zee and Van Riemsdijk, 1987, Destouni and Cvetkovic, 1992). It is very effective to capture the effect that fronts do not move with the same velocity for all places. It is used more for leaching in unsaturated soil, than for groundwater transport, because stream tubes are more unidirectional in soil (vertical) than in groundwater.

For the vertical transport of solute, the simplest model is the purely convective model. For the case of linear sorption, it is given by:

for each concentration that at time 0 was located at depth (*z*) equal to zero. For gravity leaching, the velocity is determined by the unit gradient and the hydraulic conductivity, where the latter has been demonstrated to be strongly spatially variable. The retardation factor is likewise known to vary spatially. This implies, that even if all concentrations ‘start’ at time 0 at the same position (say *z*=0), then after a certain time, they will have moved more in some stream tubes than in others. This is illustrated in Figure 9. Due to spatial variability of flow velocity and retardation factor, also the front moves with a distributed velocity!

Now, it is important to recognize with which complexity we wish to describe such spatiotemporal variable transport. In case of the PST model, we are satisfied if we know how the transport is for the entire area (for instance field) on average, and we are less interested in which part of the field the transport is fast or slow. In that case, the configuration of the solute front in *x* and *y* is of less importance, and we consider the field average front. For that front, we can be interested in the mean front position, its distribution, and its variance. In this chapter, we focus on the mean position and the distribution around this mean.

The field average front can be easily determined using Monte Carlo simulation, where random numbers are drawn for the distributed parameters, combined in the convective model, to determine the distributed front position at designated time (Van der Zee and Van Riemsdijk, 1987). However, it is also possible to determine the field average front analytically, in case that the parameters *v* and *R* are lognormally distributed. This analytical procedure is as follows, for the simple case of a Heaviside solute input at the soil surface.

First, we experimentally or theoretically determine the mean (*m*) and standard deviation (*s*) of the lognormally transformed *v* and *R*, i.e. of ln(*v*) and ln(*R*). The so-called reproductive properties of the lognormal distribution are for the purely convective transport model of use to determine the statistics of front depth: if *z*=*vt*/*R*, then:

As for the Heaviside input, the concentration at the soil surface changes abruptly from one value to another, we obtain different results depending on the initial and boundary conditions. As a first step, we determine, for a distributed ln(*z*) (for front depth), what is the probability that the front has not passed a certain depth. This probability is given by:

with the pdf given in section 8.3. Inserting the lognormal pdf, with statistics as just determined, we obtain:

This probability is the probability that the front has not passed the depth z^{*}. To interpret that as a concentration, we have to reason what that means. If the front has not passed a certain depth, at that depth we still find the original (initial) concentration. Hence, *Pr* represents the fraction of the area, where we still find the initial concentration at depth *z*=*z*^{*}. We now consider two simplified systems. In the first case, the initial concentration is 1 and the incoming concentration is 0. In that case, the areally averaged concentration at depth *z*^{*} is:

If, on the other hand, the initial concentration is equal to 0 and the incoming concentration is 1, then:

More generally, if the initial and final (or input) concentrations are not equal to 0 or 1, the mean concentration can be transformed in the real value, using transformations such as:

Due to the logarithmic transforms, the mean concentration profile commonly looks as in Figure 10.

In the presentation so far, we assumed that the important properties are lognormally distributed, i.e., taking the logarithmic transform, we obtain numbers that are normally distributed. Often, this assumption is appropriate (Van der Zee et al., 1988), but also convenient, if the convective PST model is adopted: the reproductive properties hold for multiplication and division in case of a lognormal pdf, whereas they are additive/subtractive for the normal pdf: see Eq. (119). Also, lognormally distributed parameters cannot be negative, whereas normally distributed parameters can. This is important, because for instance the hydraulic conductivity and the retardation factor cannot be negative.

Although transport in the water unsaturated soil is predominantly vertically downward or upward, some horizontal displacement always happens. This transversal transport can also be important, for instance, if a slug of contaminant moves downwards, but dispersional mixing and movement in the transversal directions occurs for a contamination event of limited areal extent (in the horizontal x and y directions). In that case, the explained theory is valid but spatial moments have to be calculated in the other directions than only depths, as in the previous sections.

### 7.6. Multidimensional stochastic modelling

In the period 1980-1995, stochastic groundwater hydrology made major advances. The basis of those advances was that the saturated hydraulic conductivity was a lognormally distributed random space function RSF, that could be characterized with e.g. an exponential autocorrelation function. For this assumption, analytical expressions were derived for the hydraulic head and the specific discharge RSF. Water flow understanding was compiled in the milestone book by Gedeon Dagan (1987), and an insightful paper on water flow in unsaturated soil was written by Kurt Roth (1995), with attractive illustrations of water channeling due to hydraulic conductivity variability, under steady state flow.

Whereas numerical and analytical work has been done with regard to solute transport in 2D heterogeneous conductivity fields, only in 1993 two papers came out that addressed transport of solutes that adsorb linearly, hence with RSF for the hydraulic conductivity as well as for the (linear) adsorption coefficient, and retardation factor (Bellin et al., 1993, Bosma et al., 1993). Results were produced in terms of second central moments of the displacement distance.

At the moment, the techniques to generate RSF for hydraulic parameters has advanced enormously compared to the early days of stochastic groundwater hydrology and contaminant hydrology. Therefore, it has become possible to look at real cases of soil and groundwater contamination, taking random variability of properties into account. This versatility is made a bit more modest, as it remains a critical, often subjective and much debated issue of which values to give to core properties, such as the dispersivities.

## 8. Conclusion

Solute transport has originally been considered for the local scale, where small scale heterogeneity resulted in a diffusion-dispersion type of mixing that at the scale of interest (often columns of porous media, soil between soil surface and groundwater level, or aquifers) resulted in a normally distributed concentration in space and time. Because of this normal distribution, this mixing could be described by a normally (or Gaussian, Fickian) distributed process, which is in full agreement with a linear parabolic partial differential equation.

In the past few decades, particularly between 1979, with the pioneering work of Bresler and Dagan (1979), and today, it has become apparent, that also at larger scales than the sand grain, heterogeneity leads to variability in flow and transport. This variability has been approached predominantly from a stochastic point of view, by considering it as a spatially distributed stochastic process. This approach appeared to be a rewarding one to address both real mixing and uncertainty at scales of interest: the decision and management scale. At the same time, it is still debatable which values of the dispersivities are appropriate, and in view of the theory provided in this chapter, section 3, this also holds for the discretization used in numerical modelling. For instance, the stochastic theory for water flow and solute transport resulted in equations for the macro dispersivities. However, these dispersivities do not necessarily represent real mixing. Briefly, this issue is discussed by both Janssen et al. (2006) and Eeman et al. (2012).

Whereas conventional solute transport often considered only one scale of heterogeneity (the grain or soil sample scale), the stochastic approach addressed so far mainly two scales, the microscopic and one larger scale, the last characterized by statistics of macroscopic properties such as the hydraulic conductivity or the retardation factor. In reality, we have to deal with a whole hierarchy of scales, from grain, to sample, horizon/layer, geological strata, to watershed. Comprehensive theory cannot address all these scales in a simple theory, so it considers those that are deemed most important.

Scientifically, further advances are being made. Whereas the mentioned work in section 7 particularly addresses spatial variability, in the research on transport of tracers and reactive solutes in soil and groundwater, also the variability in time is getting more interest. A prominent example in this respect has been that since 2000, eco-hydrological theory has been added to the spectrum, that emphasizes variability as a function of time, particularly through atmospheric forcing: rainfall varies erratically as a function of time (Rodriguez-Iturbe and Porporato, 2004, Vervoort and Van der Zee, 2008). In near future, these developments of, on one side, spatial, and on the other side temporal variability, will be no doubt combined.

In soil and groundwater management, modelling as discussed in this chapter is essential to prioritize decision making. Recently, a paper considered policy making in the Netherlands (Witte et al., 2012), and was quite influential nationally and internationally. This paper emphasized local variations of several conditions that affect the fate of factors that influence ecological impacts. It is a clear example of the potential and limitations of transport modelling.

Meanwhile, throughout the world, this type of modelling is central in risk assessments of soil and groundwater contamination, in dimensioning soil and groundwater remediation, and in evaluations of how policy decisions work out in practice for the cases of changed and unchanged policies. For these reasons, awareness and in some cases experience, with models as described in this chapter, is important.

## Appendix A: Stability analysis of explicit finite difference transport equation

Consider the explicit finite difference approximation of the 1-dimensional solute transport equation, where the advective term has been approximated by a central difference:

Collecting terms, this can be rewritten as:

where:

Now assume that at a certain time a small error *ε* (perturbation) is introduced in the solution to the equations. That can e.g. be caused by roundoff errors in a computer calculation. If we indicate the “correct” solution with *C* in the finite differencd equation (126) leads to:

Since the correct solution *ε*:

The fact that the perturbation *ε* is given by the same equation as the correct value of *C* is caused by the fact that the equation in *C* is a linear one.

Now consider one Fourrier component of the perturbation given by:

where *ω* is the wavenumber, and *λ* the time dependent amplification factor. In order to have a stable solution, we will require that any perturbation *ε* in the solution will decrease with time for any wave number. Basically this means that *λ* will decrease with time.

Substitution of (130) in (129) gives:

Division by λ^{n}e^{iωx} then gives:

Now we will use the following relations:

Substitution in equation (132) then leads to:

If we now take the abolute value, and substitute for shorthand θ=ωΔx, the following relation is obtained:

From well known goniometric relations:

Substitution in equation (135) gives, after rearranging terms:

In order to have a stable solution, the absolute value of the ratio of λ^{n}/λ^{o} should be <1. Bearing in mind that both sin^{2} and sin^{4} are always positive, it is sufficient to require that:

The second relation in (138) can be written as:

and the first relation as:

Note, that in order to obtain a stable solution, both requirements (A-15) and (A-16) need to be fulfilled.