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
Concentrations of species in the water phase Ci (including pure water itself) are defined as the mass of the species per unit volume:
The density of a multi-component fluid, consisting of N components, is then given as:
For dilute solutions (tracer concentrations) all mass fractions
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
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’qzΔxΔyΔz
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
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:
Define a piezometric head
Basically, the piezometric head consists of a pressure head
and similar expressions can be obtained for
Consequently, if the density
and the same for the y and z direction. This shows that the hydraulic conductivity
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
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:
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
Note, that the average pore water velocity
2. Simplified description of processes in reactive transport
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
The net mass influx in the x-direction over a period
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
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
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.
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:
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
(find the average water velocity by integrating the water velocity from 0 to
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
Equation (22) can be written in dimensionless form by defining the following dimensionless variables:
We can now choose any one of the characteristic values
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.
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:
In case only diffusion occurs, the mass balance equation reads:
Now consider a 1-dimensional mass balance equation with constant porosity
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:
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
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
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:
In a fully 3-dimensional system, with velocity components
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
For mechanical dispersion, the longitudinal dispersivity
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:
As a guideline, the transversal dispersivity is commonly assumed to be 5 to 10% of the longitudinal dispersivity.
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:
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:
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
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
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
Assume we have a mass of soil
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.
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:
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:
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
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
For the evaluation of the spatial derivatives in a finite difference approach, we use the following Taylor series expansion:
The first order derivative of
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
For the time derivative the following approximation can be obtained:
In the following, we will adopt the notation:
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:
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,
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
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
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
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
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
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
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
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
As a consequence
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
Condition (70) is often given in a slightly different form:
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:
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
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:
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.
Equilibrium adsorption/desorption, especially for larger concentrations, can often be described by either a Freundlich equation (Figure 5, cases B and C)
For the non-linear adsorption, we can still define a retardation factor:
which will, however, be dependent on the concentration
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:
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).
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:
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:
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:
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
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:
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,
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:
Equation of state:
These equations can be made dimensionless by choosing appropriate reference values for the different variables in the equations:
Now, if we choose the following reference values:
the system of equations reduces to:
where it is understood that all variables are dimensionless (subscript
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
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
The decay constant
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
and for integration to
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 () instead of concentration.
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 0th moment is needed to obtain a pdf. The first spatial moment divided by the 0th spatial moment is known as the average position of the plume of solute (
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
where the mean
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
So, to calculate the mean, you first have to calculate the variance (of
Median and mode for the normally distributed
If you know the statistics of
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
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 (
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
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
First, we experimentally or theoretically determine the mean (
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(
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,
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.
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:
Now assume that at a certain time a small error
Since the correct solution obeys equation (126), we obtain the following equation for the perturbation
The fact that the perturbation
Now consider one Fourrier component of the perturbation given by:
Division by λneiω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 sin2 and sin4 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.