This is an updated version of Joachim Gruber, "Contaminant accumulation during transport through porous media", Water Resourc. Res., 26, 99 - 107, 1990.

The analysis suggests that the potential hazard associated with a contaminated soil or sediment cannot be derived from the contaminant adsorption isotherms alone. The additional data needed are the velocities of the possible waves.

The spectrum of examples is as diverse as the nature of possible local processes:

- In the metallurgical process of zone refining, a solid is purified by passing a liquid zone through it which collects the unwanted components [Pfann, 1958].
- Hunt and coworkers [Hunt, Sitar, and Udell, 1988] have shown that a similar process can be used to clean up hydrocarbon spills in soils: The local process in this case is evaporation of the trapped immobile contaminant phase by a steam front.
- In geology, precipitation/dissolution at migrating fronts leads to ore body formation [Lichtner, 1985; Walsh, Bryant, Schechter and Lake, 1984], skarns [Guy, 1981] and other geochemical self-organization processes [Ortoleva et al., 1986, 1987].
- Frontal analysis and displacement development [Glueckauf, 1955] are examples from the field of multicomponent ion chromatography. Here, the local process is ion exchange. The latter processes accumulate traces.

These examples have the following in common: A local process changes the ratio C

This paper is concerned with adsorption as a local process and aqueous solute transport in a porous medium containing a deposit of adsorbed M. M accumulates at the fronts of waves of solution composition changes propagating through the deposit and triggering spatially correlated remobilization [Gruber, 1994].

The concentration maxima observed in soils by Killey et al. [1984], Serne, Peterson and Gee [1983], Peterson, Martin and Serne [1986], Dubrovsky et al. [1985] and White et al. [1984] can be interpreted as a result of such remobilizing waves.

The former group [Killey et al., 1984] investigated a desposit of Co-60 near a low-level waste infiltration pit at Chalk-River Nuclear Laboratories, Ontario, Canada. About two decades before the authors started their investigations of the contaminated porous medium (a sandy aquifer), nitric, tartaric, and oxalic acids had been discharged at the waste site. The then five years old Co-60 deposit had thus been remobilized by complexation. The authors report fronts of elevated concentrations of Co-60 at some depth beneath the surface.

Of the latter groups Serne, Peterson, and Gee [1983] and Peterson, Martin and Serne [1986] worked with laboratory soil columns and the other groups [Dubrovsky et al., 1985; White et al., 1984] reported on mill tailings. They observed the effect of a migrating low-pH front on the spatial element distribution in the tailing piles. The front acts as "chemical dam" [Guy, 1986], behind which elements pile up, among others cobalt, lead, cadmium, zinc, nickel and the uranyl ion. Similar accumulations arise when exhausted chromatographic columns are regenerated by infiltration of a remobilizing solute [Helfferich and Klein, 1970]. This demonstrates that precipitation as well as sorption can lead to element reconcentration.

If in addition to the physical laws of solute transport the equilibrium adsorption isotherms C

Rather than presenting site specific transport simulations of accumulation scenarios, this paper will discuss which chemical factors lead to the development of accumulations.

The first part of this paper reviews one possible mechanism of accumulation of a trace M ("contaminant"), known in chromatography [Helfferich and Klein, 1970]. Here the terminology will be defined, and the description will use basic concepts and the language of multicomponent ion chromatography with interference [Hwang, Helfferich, and Leu, 1988]. The accumulation conditions will be stated in terms of distribution coefficients kd of M and a suitable secant of the isotherm of the remobilizing chemical f. These system properties can either be taken from empirical compilations or can be computed based on a suitable adsorption model.

In the central part of this paper ("Accumulation Potential") it is shown that once these simple system properties are known, one can decide whether or not a constant kd model (or a model assuming a stochastic kd distribution) may be used to compute contaminant transport, thus ignoring the possibility of accumulation. Contaminant kd's alone do not provide sufficient information to make that decision. Cation adsorption on oxide mineral surfaces will be used to illustrate this point.

(1)

where
C_{ti} = C_{i} + c_{i} is the total concentration of i, v_{0} is the velocity of the liquid and D is the coefficient of hydrodynamic dispersion. Influx of non soluble phases of i is not considered in (1). In this paper it will be assumed that the dispersive influx term is small compared to the advective one represented by the first term on the right hand side of (1). Disregarding the dispersive term, (1) reduces to

(2), a hyperbolic set of conservation laws

In the case of adsorption with an isotherm C_{i} = f_{i}(c_{1}, c_{2}, ..., c_{n})

where
the summation Σ_{j} is done over all chemical components j. This can be written in matrix form as

Here, **C** and **c** are vectors with components C_{i} and c_{i}, respectively. **Q**(**c**) is the Jacobi matrix of the multicomponent isotherm **C**(**c**), with elements

In this notation the hyperbolic set of conservation laws (2) becomes

(3)

with **I** being the unit matrix. (In multicomponent chromatography and contaminant transport **I** + **Q** = **R** is called "retardation matrix".)

The set of quasilinear conservation laws (3) is hyperbolic, meaning that the inverse of the matrix (**I** + **Q**(**c**))/v_{0} has n real, distinct eigenvalues λ_{j}(**c**), (j = 1, n). Customarily, they are arranged in increasing order

A line in {x, t} space (Fig. 1) with slope λ_{k}(**c**) is called "k-characteristic at concentration **c**". It is the path in this space of a concentration vector **c** being a solution of (3) corresponding to the eigenvalue λ_{k}(**c**). This solution is called k-rarefaction wave, the name referring to the properties

- "wave", i.e. the concentration
**c**varies continuously across the wave (e.g. along a path perpendicular to the characteristics) and propagates with time. Usually the variation is not periodic. - "rarefaction", i.e. the fan of characteristics spreads with increasing time t.

- λ
_{1}(**c**) with**c'**≤**c**≤**c**,_{1} - λ
_{2}(**c**) with**c**≤_{1}**c**≤**c**,_{2} - λ
_{4}(**c**) with**c**≤_{3}**c**≤**c**=_{4}**c"**.

**FIGURE 1**: Schematic of the solution of hyperbolic system (2) with 4 components (vectors **c**_{i} have 4 components: **c**_{i} = {c_{i1}, c_{i2}, c_{i3}, c_{i4}}), represented in {x, t} space. System solution consists of 4 constant states, **c**_{1}, **c**_{2}, **c**_{3}, **c**_{4} = **c"**, the characteristics of which are the heavy lines. Constant states 1 through 3 are bounded by the 1-, 2-, and 4-rarefaction waves (heavy solid lines) and the 3-shock ("line of discontinuity, L", heavy, dashed). Rarefaction waves are pictured by a fan of some arbitrarily chosen characteristics (light lines), each representing the path of a concentration **c** within the wave.

Insert: Characteristics 1', 2', 3', 3", 4" are defined as

- 1' is the line with slope λ
_{1}(**c'**), - 2' is the line with slope λ
_{2}(**c'**), - 3' is the line with slope λ
_{3}(**c'**), - 3" is the line with slope λ
_{3}(**c"**), - 4" is the line with slope λ
_{4}(**c"**) = λ_{4}(**c**)._{4}

There exist also **piecewise** smooth solutions of the system of conservation laws (2). Here the concentration changes by a finite amount rather than by differential amounts as it does across a rarefaction wave. Let P be a point on the line of discontinuity L and let **c'** and **c"** be the values of the solution on the left and right side of the discontinuity, respectively. Accordingly, the mass conservation law involves the differences

Δ**c** = **c'** - **c"**

Δ**C** = **C'** - **C"**

rather than differentials ∂**C** and ∂**c**. For every component i we have a conservation law similar to (2)

(5)

After multiplying with Δx

we can solve for Δx/Δt = σ_{k}, the speed of the k-discontinuity:

(6)

_{}

In {x, t} space, σ_{k} is the slope of the line of the k-discontinuity. Equation (6), the so-called "Rankine-Hugoniot relation", has n physical (piecewise smooth) solutions. Any such solution is called a k-shock (also called "front", "self sharpening wave" or "migrating boundary") if it meets the following "entropy" conditions (from Greek εντρεπειν (entrepein) = to turn -or run- into) [Lax, 1957]:

- Draw, issuing from P in positive t-direction, the k-characteristics of constant states
**c'**and**c"**(in the example of Fig. 1 (insert) these are the dashed thin lines 3' and 3"). The characteristics need to cross line L.

(7)λ

_{k}(**c'**) > σ_{k}> λ_{k}(**c"**), (k = fixed, in our example: k = 3).

- Draw, issuing from P in positive t-direction, those k-characteristics with respect to the state
**c'**which stay to the left of L, and those with respect to the state**c"**which stay to the right of L (thin lines in box in Fig. 1). The total number of characteristics drawn in this way needs to be n-1.

(8)λ

_{k-1}(**c'**) ≤ σ_{k}≤ λ_{k+1}(**c"**), (k = fixed, in our example: k = 3).

Because the λs are arranged in increasing order this implies

λ

_{1}(**c'**) ≤ λ_{2}(**c'**) ≤ ... λ_{k-1}(**c'**) ≤ σ_{k}≤ λ_{k+1}(**c"**) ≤ ... ≤ λ_{n}(**c"**)or in our example:

λ

_{1}(**c'**) ≤ λ_{2}(**c'**) ≤ σ_{3}≤ λ_{4}(**c"**).

Condition (7) makes sure that the k-wave "breaks", meaning that it resembles a water wave breaking on its path onto the shore, its top running faster than its bottom part. The physics of concentration waves (c ≥ 0) allows this breaking process to proceed only until all parts of the wave have arrived at the same line. At this line the concentration jumps from **c _{k-1}** to

In chromatography [Helfferich and Klein, 1970, Charbeneau, 1988] it is customary to use the matrices

The eigenvalues

- of
**I**+**Q**(**c**) are the retardations, R_{j}(**c**) with j = 1, 2, ... n, and - of
**Q**(**c**) are the generalized distribution coefficients, kd_{j}(**c**) with j = 1, 2, ... n.

(9)

λ_{j}(**c**) = v_{0}/(1 + kd_{j}(**c**)) = v_{0}/R_{j}(**c**),

thus

(10)

kd_{j}(**c**) = v_{0}/λ_{j}(**c**) - 1.

The retardation ρ_{k} of a k-shock is defined analogously to (9)

(11)

σ_{k}(**c**) = v_{0}/ρ_{k}(**c**).

The Rankine Hugoniot condition (6) tells us that the retardation of the k-shock is

(12')

ρ_{k} = 1 + ΔC/Δc.

where C is any component C_{i} of the vector **C** and c is any component c_{i} of the vector **c**. (In a k-shock ΔC_{1}/Δc_{1} = ΔC_{2}/Δc_{2} = ... ΔC_{n}/Δc_{n})

In the rest of this chapter we will apply the theory to show how a contaminant will accumulate during migration.

The retardation of this shock is

ρ_{f} = 1 + ΔC_{f}/Δc_{f}.

Because the concentration c_{M} of the contaminant M is small, trace M moves on the background of the natural components and f, without being able to influence the chemistry of the background. Its speed is approximated by λ_{M} given in (9) with the vector c in kd(c) not containing the concentration c_{M}:

(9')

λ_{M}(c) = v_{0}/(1 + kd_{M}(c)),

c = {c_{1}, ..., c_{M-1}, c_{M+1}, ..., c_{n}}

Condition (7) for the contaminant is

(14)

kd_{M}(c') ≤ ΔC_{f}/Δc_{f} = kd_{f} ≤ kd_{M}(c").

This means: An observer located at P on the remobilizing shock will see species M coming towards her both from upstream and downstream locations. The front carries the inventory of the remobilized contaminant with it. The degree of accumulation can be calculated [Gruber, 1984].

The following limiting case of this situation is interesting as well: The remobilizing front moves without retardation (σ_{f} = v_{0} = λ_{M}(c') or kd_{f}(c') = ΔC_{f}/Δc_{f} = kd_{f} = kd_{M}(c') = 0):

(15)

0 = kd_{M}(c') = kd_{f}(c') ≤ kd_{M}(c").

The contaminant is remobilized by the moving front but cannot travel across the front. Thus, the remobilized contaminant inventory will stay in solution.

Whenever (14) or (15) is met, the decrease of the distribution coefficient (14) must not be disregarded, and transport models based on the assumption of a constant kd or a spectrum of stochastic kd 's are inappropriate. In that case the transport model has to contain a chemical adsorption model. Any of the chemical models meeting that requirement can be chosen, if it is used within the composition space in which it has been calibrated.

In surface complexation models adsorption is described as complex formation on hydrolyzed surface sites, S-OH, for example

(16)

S-OH + M^{++} + H_{2}O = S-OMOH + 2H^{+}.

Consistent with the assumption of M being a trace element, only mononuclear surface complexes of M need to be considered. The concentration of a surface complex s containing M, C_{M,s}, is then expressed as the product of a side reaction coefficient, σ_{s}, and the free concentration of the adsorbing ion, [M^{++} ], (definition of σ_{s} ):

(17)

C_{M,s} = σ_{s}[M^{++} ].

Similarly, the concentration of a soluble complex mononuclear in M, [Q_{l}], is the product of a corresponding side reaction coefficient, σ_{l}, and [M^{++} ]. This coefficient is a function of the concentration of a background element L that forms a soluble complex with M:

(18)

[Q_{l}] = c_{M,l} = σ_{l}(L) [M^{++} ].

The major chemical components of the system are also subject to adsorption and complexation but without restriction to mononuclear reactions. The appendix gives the sources of the data considered in the following examples.

Above pH = 7 adsorption of the divalent transition elements Pb, Cd, Cu, Zn, Ni and of UO_{2}^{++} on a metal oxide can be represented by one single surface complexation reaction of type (16) [Hohl and Stumm, 1976; Vuceta, 1976, James and MacNaughton, 1977; Benjamin, 1978; Davis and Leckie, 1978; Richter, 1978; Balistrieri and Murray, 1979, 1981, 1982; Altmann, 1984; Tripathi, 1984]. The corresponding adsorption side reaction coefficient is

(19)

where {H^{+}} is the activity of the proton in solution. Because of its surface charge, the oxide is surrounded by an electrostatic potential field, Ψ. The exponential term in (19) takes account of the work necessary to assemble the surface complex. There are two major adsorbing locations on the surface, the 0-plane (the surface itself), which loses one proton upon adsorption of one M, and the β-plane (located at some distance from the surface), which receives the (single) positive charge of the MOH complex. e, Ψ_{0} and Ψ_{β} are the charge of the electron and the values of the electrostatic potential in the 0- and β-plane, respectively. K is the equilibrium surface complex formation constant. kT is the product of the Boltzmann constant and the absolute temperature of the system. With

(20)

C_{M} = σ_{s} [M^{++} ]

and Σ_{l} meaning the summation over l, the soluble complexes of M and

(21)

c_{M} = (1 + Σ_{l} σ_{l} ) [M^{++} ]

the distribution coefficient of M can be rewritten as a ratio of the side reaction coefficients:

(22)

Combining (22) with (19) gives the final expression for kd of M:

(23, Triple Layer Model)

Major processes and environmental variables influencing kd are now represented by individual terms in (23):

- (a) soluble complex formation (Σ
_{l}in denominator),

- (b) pH ({H
^{+}}^{2}in denominator),

- (c) electrostatic repulsion or attraction due to the presence of major ions, such as H
^{+}, Ca^{2+}, or SO_{4}^{2-}, in the 0- or β-planes (exponential term),

- (d) competition for adsorption sites ([S-OH]).

When for example major ions in solution forming surface complexes [S-X

(24)

S_{r} = [S-OH] + Σ_{l} [S-X_{l}] + σ_{s} [M^{++}].

The right coordinate represents the factor S

S

(25)

S_{r} = n_{s} f v ρ (1 - Φ) = 0.1 moles/liter.

where n_{s} is the number of sites per unit area of the oxide. To simplify the discussion n_{s} = 10^{-5} moles/m^{2} is used for all oxides in Fig. 2, because the variation of S_{r} due to the variation of n_{s} [James and Parks, 1982] is small compared to the variability of the other terms in (23). f is the fraction of the soil particle covered with the oxide (0.2 is assumed). v is the surface to weight ratio of the soil particle (30 m^{2}/g is assumed). ρ is the density of the soil particle (2.6 g/cm^{3} assumed) and Φ is the average volumetric water content of the porous medium (to simplify notation, Φ = 1 was used in chapter "Mechanism of Accumulation"). S_{r} = 0.1 moles/liter is consistent with the estimate used by Anderson [1979].

**FIGURE 2.** Distribution coefficients (kd, right ordinate, (-) means dimensionless) calculated from adsorption constants (K, left ordinate) for 6 divalent cations and 5 oxide adsorbents ("contaminants M"). Data and assumptions (see (23)):

- pH = 7,
- [S-OH] = S
_{r}= 0.1 moles/l (major ions in solution do not compete with contaminant M for adsorption sites), - σ
_{l}= 0 (no soluble complex formation), and - exponential term exp{-e(Ψ
_{β}-Ψ_{0})/kT} = 1 (no electrostatic forces on M due to the presence of major ions).

Accumulation of M occurs when a remobilizing agent shifts the corresponding symbol in Fig. 2 from its pre-remobilization position kd" = kd_{M}(c") (e.g. the one used in the figure) down to a position kd' = kd_{M}(c'). According to (14) the kd of the shock kd_{f} = ΔC _{f}/Δc _{f} is intermediate between kd' and kd". So for the majority of cations M and adsorbents in Fig. 2 retardations of shocks range between 1 and 100:

1 ≤ ρ_{f} ≤ 100.

The remainder of this chapter will show that such shifts of the symbol locations in Fig. 2 are well possible (processes (a) - (d)). The element nickel is chosen as a representative example (because (16) describes adsorption in the pH range of interest, 4 <= pH <= 10 [Richter, 1978]). The computer code used is MINEQL-SGMA (Stanford Generalized Model of Adsorption), a computerized version of the surface complexation model [Davis, James and Leckie, 1978; Davis and Leckie, 1978].

*Processes (a), (b):*

Reducing the pH in Fig. 2 by one unit shifts the symbol of Ni downward by two orders of magnitude (Fig. 3). Increasing the pH by one unit moves the symbol upward and will therefore not remobilize nickel. Among the natural ligands other than OH^{-}capable of complexing contaminants M, carbonate is likely to be most effective in reducing kd: For the group of divalent cations in Fig. 2 the bracket in the denominator of (23) can vary between 1 and approximately 10^{3}under natural conditions [Turner, Whitfield and Dickson, 1981]. But also phosphate and sulfate concentrations above 0.01 moles/l and traces of organic complexing agents, for example EDTA, have the potential of reducing kd by orders of magnitude.

**FIGURE 3.**Distribution coefficient of nickel, calculated as a function of pH and major ion concentrations. Numbers on curves indicate composition of water, as explained in Tab. 1. The other data used in the computations are given in the appendix. See also Fig. 1 in [Gruber, 1988].

*Process (c):*(regard only the curves for EDTA-free pore waters, i.e curves labelled "no EDTA")

Electrostatic forces in the 0- or β-planes are less effective in remobilizing the considered contaminant cations. Therefore all (no EDTA) curves in Fig. 3 are close to each other.In detail:

- Curve 2 gives the computed nickel kd for a very dilute pore water, i.e. for a positive charge strongly repelling (0-plane in the) oxide surface (the β-plane being empty).
- All other curves represent nickel adsorption in the presence of elevated anion (chloride, sulfate) and cation (calcium) concentrations (Tab. 1).

These anions and cations reside in the 0- or β-plane and shield or buffer the 0-plane charges present in very dilute pore waters.- Curves 1 and 3 show that the net effect of negative ions adsorbed (in the β-plane) onto goethite is an increase of kd by up to one order of magnitude. Animation of accumulation of strontium based on this effect.
- Curve 4: The large calcium concentration in an equimolar solution of calcium and carbonate in equilibrium with calcite and the resulting large concentration of adsorbed calcium ions changes the exponential term in (23) by a similar factor as the presence of the mentioned negative ions, but in the opposite direction.

**FIGURE 4.**Distribution coefficient of nickel, calculated as a function of pH for the major ion concentrations of the undisturbed environment (curve IC) and of the remobilizing water (curve Remob) in Fig. 5 and 6.

*Process (d):*

Competition for adsorbing sites by major cations, e.g. calcium or magnesium on the one hand and M on the other hand is negligible when the total concentration of adsorbing sites S_{r}is considerably larger than the largest total concentration C_{i}of a competing cation. This is the case of Fig. 3. In Fig. 4, however, S_{r}is two orders of magnitude smaller than in Fig. 3. In addition, in calculating curve Remob in Fig. 4, precipitation of calcite at higher pH has been precluded. This increases the concentration of the competing calcium ions in solution above the one of curve IC. Thus, especially the bending of curve Remob away from curve IC shows us the extent to which competition influences adsorption.

EDTA-curves in Fig. 3:

These curves give a quantitative impression of another form of competition, the competition for complexing ions, in this case EDTA. In the presence of enough calcium, EDTA is bound largely to this cation. Consequently, less EDTA is available to complex the contaminant. Therefore, when a small amount of EDTA is added to a calcium rich solution, the kd of the contaminant is not changed nearly as much as in the absence of calcium (compare the shift of curve 4 with the ones of curves 2 and 3).

Finally, also the model parameter K and the exponent of {H^{+}} in (19) vary slightly with changing environmental conditions ("residual variation") [Benjamin, 1978; Honeyman, 1984; Altmann, 1984].

**TABLE 1.** (Negative logarithm of) Total concentrations of major ions (moles/l) used in computations of Figs. 3 - 6. Total concentration of adsorbing sites is 0.25 moles/l except for Figs. 4 and 6, where it is 0.0025 moles/l. Total concentration of Ni is 10^{-8} moles/l, except in Fig. 5 (intruding water), where it is 3.8 x 10^{-10} moles/l.

Figure |
Ca |
Na |
CO_{3} |
Cl |
SO_{4} |
pH |

3 and 4 (curves 1, Remob) | 2 | 6 | 6 | 2 | 6 | - |

3 and 4 (curves 2, IC) | 6 | 6 | 6 | 6 | 6 | - |

3 (curve 3) | 6 | 2 | 6 | 2 | 2 | - |

3 (curve 4) | * | 6 | * | 6 | 2 | - |

5 (intruding water) | 6 | 6 | 6 | 6 | 6 | 4 |

5 (undisturbed environment) | # | 2 | # | 2 | 3 | 7 |

6 (intruding water) | 2 | 6 | 6 | 2 | 6 | 8 |

6 (undisturbed environment) | 6 | 6 | 6 | 6 | 6 | 8 |

Legend:

- - : 3 <= pH <= 10
- * : water is saturated with calcite, total calcium and carbonate concentrations are equal.
- # : same as *, with calcite concentration in cells = 3 10
^{-4}moles/l.

Based on the preceding quantitative discussion of processes (a) - (d) one can expect critical kd shifts in Fig. 2 (i.e. the ones necessary to meet accumulation conditions (14)) even in the case of a poorly retarded front (ΔC _{f}/Δc _{f} = kd_{f} (appr.=) 0). The deposits on γ-Al_{2}O_{3} and amorphous Fe_{2}O_{3} will resist accumulation less than the ones on α-SiO_{2}, because the former are more easily remobilized.

In this paper a very crude difference scheme will be shown to be able to illustrate two simple cases (a much more sophisticated model based on the philosophy presented here has been published by Amir and Kern, 2010, in cache). In both cases, remobilization of a contaminant deposit leads to a migrating contaminant wave carrying the remobilized inventory. The difference scheme meets the requirements of von Neumann stability and small truncation error [Lax and Wendroff, 1960], while using a negligible fraction of the total computing time.

The set of difference equations used in the model (MINTRANS) is

(26)

(C_{tik}^{m+1} - C_{tik}^{m}) / Dt + v_{0} (c_{ik}^{m} - c_{ik-1}^{m})/Dx = 0.

k denotes a position and m a time step:

c_{i}(kDx, mDt) = c_{ik}^{m}.

The set (26) expresses the mass balance for transport in a chain of cells k, each cell having a length Dx and spatially uniform equilibrium concentrations C_{tik}^{m}, c_{ik}^{m} (with C_{tik}^{m} = C_{ik}^{m}+c_{ik}^{m}). Dt is the time necessary to exchange the soluble phases of the cells:

(27)

Dx/Dt = v_{0} ≥ v_{j} for j = 1, 2, ..., n.

v_{j} is any possible concentration velocity in the chain. Because of (27) the scheme (26) satisfies the Courant-Friedrichs-Lewy convergence condition [Lax and Wendroff, 1960].

Equation (26) is an appropriate difference scheme for (2) in the following two situations:

*Vertical transport in unsaturated homogeneous field plots under periodic rain or irrigation.*

Jury and coworkers [Jury, 1982; Jury and Sposito, 1985; Butters, Jury and Ernst, 1989] have shown for their field that up to a depth of several meters the solute can be viewed as moving in a system of non-interacting stream tubes or flow paths. The lengths of the stream tubes leading to a given depth are statistically distributed. The solute moves only when water enters a flow path at the soil surface. In the time between rain events the solute is stationary. If the duration Dt of single rain events is short compared to the interval between them [Eagleson, 1978] and compared to desorption times, the spatial concentration profile (along a flow path) is better represented by a series of steps than by a continuous function. The distance between adjacent steps is Dx = i_{o}Dt/Φ, where i_{o}Dt is the height per rain storm. A stream tube can then be approximated by a chain of cells with uniform concentrations, and (26) is applicable for every flow path. Transport in the entire field plot is the linear superposition of the transport through individual stream tubes.

*Transport over microscopic distances.*

While the set of difference equations (26) describes a macroscopic profile in the above case of a series of discrete rain events, it illustrates the microscopic profile when water input is continuous. In this case it is the discrete spatial structure of the solid matrix of the porous medium which allows the simple set of difference equations.

The set of hyperbolic differential equations (2) is valid for continuous space. Let d be the distance at which the solution of (2) breaks down because of the discrete structure of the pore network. In this paper it is assumed that solute transport within d is random. This includes that the starting and ending points of transport trajectories larger than d are randomly distributed within volumes V_{min}= F d, where F is the cross section of a flow channel. This causes mixing of the pore water within d. The mixing process is assumed not to allow the development of any structure within V_{min}. (The cells V_{min}resemble "continuously stirred reactors"). The microscopic concentration profile along a flow path is therefore rather a sequence of steps than a continuous curve, and consequently the set of mass balance equations (2) can be approximated by (26) with Dx = d.

- In the first example condition (14) is met,
- in the second example condition (15) is approximately valid (ΔC
_{f}/Δc_{f}= kd_{f}(appr.=) 0).

The cells contain an adsorbent (goethite) and an aqueous solution (no gas phase). The chemical components are Ca^{2+}, Na^{+}, CO_{3}^{2-}, Cl^{-}, SO_{4}^{2-}, and Ni^{2+}, the latter being the contaminant. The appendix gives the data used in the calculations. The chemical model is the triple layer model as implemented in MINEQL-SGMA. The initial states of the chain of cells is given in Tab. 1.

- 1 <= k <= 24 (initial deposit) : C
_{tik}^{o}= concentration as given for "undisturbed environment". Very little nickel is found in solution, the nickel inventory being almost entirely adsorbed onto the goethite surfaces, i.e.c

_{Ni,k}^{o}<< C_{Ni,k}^{o} - -20 <= k <= 0 (remobilizing water): C
_{tik}^{o}= concentration as given for "intruding water".

In the first scenario (Figs. 5a and b)

- remobilization is caused by water with pH = 4. The intruding water is assumed to be undersaturated with respect to calcite.
- The low-pH front is retarded by a precipitation/ dissolution buffer: 3 x 10-4 moles of calcite per liter of system volume are assumed to be initially present in each cell of the deposit.

The dissolution of calcite in the system (closed to the atmosphere) consumes protons and therefore raises the pH. When calcite has totally disappeared in a cell, the pH drops to the value in the incoming water and the contaminant (nickel) is mobilized. Next the remobilized inventory crosses the pH step and becomes re-adsorbed immediately thereafter. Thus nickel accumulates at the downstream flank of the pH step.

At pH = 4 the kd of nickel is well below 1, thus the accumulated nickel can easily follow the relatively slowly advancing pH front (condition (14) is met). In a more realistic model, accumulation would be limited by precipitation of Ni and other processes not considered here, such as the concentration dependency of the parameters of the chemical model, dispersion or kinetic effects.

**FIGURE 5.** Spatial concentration, kd and pH distribution in a nickel deposit adsorbed on goethite. Environment contains calcite as pH buffer. Upper part of figure: undisturbed environment. Lower part: environment after intrusion of remobilizing water having pH = 4 (see Tab. 1). The front of the intruding water has arrived in cell 21. Dotted curve gives C_{Ni}(c)/c_{Ni}(c). Curve representing soluble nickel coincides with abscissa. Shaded area marked (-) is remobilized Ni-inventory, shaded area marked (+) is accumulated Ni-inventory (secondary repository).

Contrary to scenario 1, in scenario 2 (Figs. 6a and b) the remobilizing front is poorly retarded. The intruding remobilizing water has the same composition as the water in the undisturbed deposit except for its higher calcium chloride concentration (10^{-2} M). Fig. 4 shows the variation of the kd of the contaminant in the deposit prior to and after remobilization (curves IC and Remob, respectively).

As mentioned above, calcium competes here with the contaminant for adsorption sites at high pH (process (d), above), because C_{tCa} (appr.=) S_{r}. At low pH, the adsorbed calcium ions affect contaminant adsorption via the electrostatic field they build up at the surface (process (c), above).

When the calcium rich water intrudes into the initially calcium poor deposit, calcium adsorbs and increases the (positive charge repelling) surface potential. This, in turn, leads to the desorption of protons into the pore water. The result lies between two extremes, depending on the efficiency of the proton buffer in the system:

- The system has an effective proton buffer: Each adsorbed calcium ion displaces one proton from the surface. The desorbed protons disappear in the buffer pool. The concentration of protons that this pool will have to accommodate is approximately equal to the total concentration of soluble calcium ions in the intruding water, C
_{tCa}, because at the pH of interest (pH > 7) a large fraction of the intruding calcium ions become adsorbed. - The system has no proton buffer: Because C
_{tCa}is large compared with the activity of the free proton, the pH drops rapidly. This stops desorption of protons even before the surface proton concentration has changed by more than a small fraction of C_{tCa}.

Solute transport is considered here for case 2. In Fig. 4 one can see that kd_{Ni} - 1 << 1, i.e. the remobilization meets conditions (15), indicating that a soluble concentration wave of the contaminant is to be expected (in contrast to an adsorbed one as in the previous scenario). This is indeed what the simulations show (Fig. 6).

A similar scenario would be one in which the goethite surfaces release protons upon desorption of sulfate or chloride ions. This happens when the intruding (remobilizing) water has a smaller concentration of these anions than the initial pore water in the deposit.

**FIGURE 6.** Spatial concentration, kd and pH distribution in a nickel deposit adsorbed on goethite, no pH buffer. Upper part of figure: undisturbed environment. Lower part: deposit after intrusion of remobilizing water. Intruding water front is in cell 21. pH in intruding water is the same as in undisturbed environment. Ni inventory has been remobilized (shaded area, -) and put into solution (shaded area, +).

Oxides of iron, aluminum, titanium, and silicium are of particular interest, because their adsorption properties vary strongly with solution composition. The waves of major species, such as protons, released from oxide surfaces upon disturbance by infiltration of remobilizing water can amplify contaminant remobilization if buffers are insufficient.

An experimental example of a contaminant wave has possibly been shown in silicate soil underlying contaminated soils (S. Eberle (local link), see also: Sickerwasserprognose in "Sickerwässern auf der Spur: Tagung der Gesellschaft für Hydrogeologie über Wasserwirtschaft und Umweltschutz in München", Deutschlandfunk, Forschung Aktuell, 21.4.2004)

This paper has shown that investigations of the groundwater pollution potential of contaminated adsorbents have to start by determining which remobilizing fronts can develop in the environments of concern and with what velocities v_{f} these fronts might migrate. Subsequently, the compositions upstream and downstream of these fronts have to be used in batch experiments to determine the contaminant adsorption isotherms. This sequence of investigations will enable us to single out critical remobilization processes, i.e. those leading to contaminant build-up at migrating remobilizing fronts.

C_{j} = adsorbed concentration of component j, unit: M, first time used

C_{tj} = total concentration of component j = c_{j} + C_{j}, unit: M, first time used

**c** = {c_{1}, c_{2}, ..., c_{n}} = vector of soluble concentrations, unit: M, first time used

c_{j} = soluble concentration of component j, unit: M, first time used

c' = soluble concentration upstream of the remobilizing front, unit: M, first time used

c" = soluble concentration downstream of the remobilizing front, unit: M, first time used

D = combined coefficient of diffusion and hydrodynamic dispersion, unit: m2/s, first time used

Δ = difference operator, Δc = c' - c", ΔC = C' - C", first time used

D = cell size operators in difference scheme (26), Dx = size of spatial cell, Dt = size of time step, first time used

∂ = partial derivative, first time used

f = remobilizing agent, first time used

f_{j} = multicomponent adsorption isotherm of component j, unit: M, first time used

I = unit matrix, an (n x n) diagonal matrix with 1's in the diagonal, 0's off diagonal, first time used

i = index meaning chemical component, first time used

K = equilibrium surface complexation constant, unit: M, first time used

k = index for wave, k = 1, 2, ..., n, first time used

k-characteristic = line in {x, t} space on which c is constant, i.e. distance x that a fixed concentration c_{j} travels on k-wave during time t, first time used

kd_{k}(c) = k-eigenvalue of **Q**(**c**), generalized distribution coefficient in k-wave at concentration **c** defined as normalized adsorbed concentration (normalized with respect to soluble concentration), (unit: L/kg), first time used

kd = distribution coefficient of contaminant M, dimensionless, first time used

k-rarefaction wave = continuous solution of (3), first time usedkd = (μg of contaminant M on solid phase per liter of system volume) / (μg of contaminant M in liquid phase per liter of system volume), system volume = volume of a cell of porous medium, containing a solid phase and a liquid phase, the volume fraction occupied by the liquid phase being φ, the volume fraction of the solid phase being 1 - φ. Note:

The partitioning coefficient Kd = (μg of contaminant M per kg of solid) / (μg of contaminant M per liter of liquid at equilibrium with the solid), dimensions: L/kgKd = kd φ/(ρ (1 - φ))

Using this Kd, the retardation R = 1 + kd can be written as

R = 1 + ρ (1 - φ)/φ Kd.

For a discussion of the limited usefulness of the constant partitioning coefficient Kd see also chapter 2.3.1 "Constant Partition Coefficient (Kd) Model", page 2.16 in" "Understanding Variation in Partitioneing Coefficient, Kd, Values", Volume I: The Kd Model, Methods of Measurement, and Application of Chemical Reaction Codes, United States Environmental Protection Agency, Office of Air and Radiation, EPA 402-R-99-004A, August 1999 (in Cache)

k-shock = piecewise smooth solution of (3), first time used

kT = product of Boltzmann constant and absolute temperature, unit: eV, first time used

λ_{k}(c) = speed of concentration (vector) c in k-rarefaction wave, unit: m/s, first time used

l = index meaning soluble complex, first time used

M = unit mole/L, first time used

M = contaminant, present at trace concentration, first time used

[M^{++}] = concentration of free ion M, unit: M, first time used

n = number of components in geochemical system (n = maximum of j), first time used

Q = Jacobian matrix of multicomponent isotherm, unit: dimensionless, first time used

[Q_{l}] = concentration of soluble complex l containing contaminant M, unit: M, first time used

q_{fj} = partial derivative of isotherm C_{f} with respect to concentration c_{j}, unit: dimensionless, first time used

R_{k}(c) = 1 + kd_{k}(c) = k-eigenvalue of (I + Q(c)), retardation of concentration (vector) c in k-rarefaction wave, unit: dimensionless, first time used

ρ = mass density of soil particle, unit g/cm^{3}, first time used

ρ_{k}(**c**) = retardation of concentration (vector) **c** in k-shock wave, unit: dimensionless, first time used

s = index meaning adsorbed complex, first time used

σ_{k}(**c**) = speed of concentration (vector) **c** in k-shock wave, unit: m/s, first time used

σ_{l}(L) = side reaction coefficient for formation of a soluble complex ML_{x} containing M and a background element L. It is the concentration of the soluble complex ML_{x} divided by the concentration [M^{++}] of the free (i.e. non complexed) ion M^{++}, unit: dimensionless, first time used

σ_{s} = side reaction coefficient for adsorption of M. It is the concentration of the adsorbed complex S-OMOH divided by the concentration [M^{++}] of the free ion M^{++}, unit: dimensionless, first time used

S-O = adsorbing surface site, S being the oxide lattice and O the surface oxygen layer, unit: M, first time used

S-OH = hydrolyzed adsorbing surface site, unit: M, first time used

S-OMOH = adsorbing surface site the proton of which has been replaced with the MOH complex, unit: M, first time used

S_{r} = total number of adsorption sites, unit: M, first time used

Σ_{j} = summation over j, first time used

t = time variable, unit: s, first time used

Φ = is the average volumetric water content of the porous medium (to simplify notation, Φ = 1 was used in chapter "Mechanism of Accumulation"), first time used

Ψ_{0} = electrostatic potential at surface, unit: V, first time used

Ψ_{β} = electrostatic potential in beta plane, unit: V, first time used

v_{0} = flux of water in direction of x, unit: m^{3}/(m^{2} s) = m/s, first time used

x = spatial variable, unit: m, first time used

{x, t} space = space with abscissa x and ordinate t, first time used

http://epubs.siam.org/doi/abs/10.1137/130907185

http://www.researchgate.net/publication/259044134_HYPERBOLIC_THEORY_FOR_FLOW_IN_PERMEABLE_MEDIA_WITH_PH-DEPENDENT_ADSORPTION/file/9c960529d47a873def.pdf

E-Mail: valentina.prigiobbe@austin.utexas.edu

Complexation & Precipitation | |

Element |
Reference |

Ca | 1, 2 |

Cd | 1, 3, 4, 5 |

Cu | 2, 4 |

Mg | 1, 2, 3 |

Ni | 1, 2, 3, 6, 7, 8 |

Pb | 1, 2, 3 |

Zn | 1, 2, 4 |

Cl, CO_{3}, PO_{4}, SO_{4} |
4 |

UO_{2} |
9 |

Adsorption |
||

Element |
Adsorbent |
Reference |

Ca, Mg, Na | G | 11 |

Cd, Cu, Pb, Zn | G | 12 |

Cd, Cu, Zn, Pb | A, F, Q | 13 |

Cd, Mg | T | 13 |

Ni | G | 14, 15 |

Cl, SO_{4} |
G | 11 |

UO_{2} |
G | 9 |

Legend

- G: adsorption on goethite
- F: adsorption on amorphous Fe
_{2}O_{3}

- A: adsorption on γ-Al
_{2}O_{3}

- T: adsorption on TiO
_{2}

- Q: adsorption on α-SiO
_{2}

- References
- 1: Ball, Nordstrom, and Jenne, 1980
- 2: Baes and Mesmer, 1976
- 3: Smith and Martell, 1976
- 4: Turner, Whitfield, and Dickson, 1981
- 5: Stumm and Morgan, 1981
- 6: Mattigod and Sposito, 1977
- 7: Martell and Smith, 1974
- 8: Sillen and Martell, 1964, 1971
- 9: Tripathi, 1984
- 10: James and Parks, 1982
- 11: Balistrieri and Murray, 1981
- 12: Balistrieri and Murray, 1982
- 13: Davis and Leckie, 1978
- 14: Richter, 1978
- 15: Gruber, 1988

Address of this Page

Physics Home

Last Changed: 25 May 2017

Joachim Gruber