This example demonstrates the capabilities of PHREEQC to calculate flow in a dual-porosity medium with exchange among the mobile and immobile pores via diffusion. The flexible input format and the modular definition of additional solutions and processes in PHREEQC allow inclusion of heterogeneities and various complexities within a 1D column. This example considers a column filled with small clay beads of 2 cm diameter, which act as stagnant zones. Both the first-order exchange approximation and finite differences are applied in this example, and transport of both a conservative and a retarded (by ion exchange) chemical are considered. It is furthermore shown how a heterogeneous column can be modeled by prescribing mixing factors to account for diffusion between mobile and immobile cells with keyword **
MIX**.

The example input file (table 33) is for a column with a uniform distribution of the stagnant porosity along the column. The 20 mobile cells are numbered 1-20. Each mobile cell, *
n*
, exchanges with one immobile cell, which is numbered 20 + 1 + *
n*
(cells 22-41 are immobile cells). All cells are given an identical initial solution and exchange complex, but these could be defined differently for each individual cell. It is also possible to distribute the immobile cells over the column non-uniformly, simply by omitting solutions for the stagnant cells that are not present. The connections between the mobile-zone and the stagnant-zone cells and among stagnant zone cells can be varied along the column as well, but this requires that mixing factors among the mobile and immobile cells are prescribed using the keyword **
MIX**.

As defined in table 33, the column initially contains a 1 mmol/L KNO_{
3}
solution in both the mobile and the stagnant zone (**
SOLUTION** 1-41). A NaCl-NO_{
3}
solution flows into the column (**
SOLUTION** 0). An exchange complex with 1 mmol of sites is defined for each cell (**
EXCHANGE** 1-41), and exchange coefficients are adapted to give linear retardation *
R*
= 2 for Na^{
+}
(**
EXCHANGE_SPECIES**). The first **
TRANSPORT** data block is used to define the physical and flow characteristics of the first transport simulation. The column is 2 m in length and is discretized in 20 cells (**
-cells**) of 0.1 m (**
-length**). A pulse of 5 shifts (**
-shifts**) of the infilling solution (**
SOLUTION** 0) is introduced into the column. The length of time for each shift is 3600 s (**
-timest**), which results in a velocity in the mobile pores *
v*
_{
m}
= 0.1 / 3600 = 2.78e-5 m/s. The dispersivity is set to 0.015 m for all cells (**
-disp**). The diffusion coefficient is set to 0.0 (**
-diffc**).

The stagnant/mobile interchange is defined using the first-order exchange approximation. The stagnant zone consists of spheres with radius *
r*
= 0.01 m, diffusion coefficient *
D*
_{
e}
= 3.e-10 m^{
2}
/s, and shape factor
= 0.21 according to table 1 (see "Transport in Dual Porosity Media"). These variables give an exchange factor
= 6.8e-6 s^{
-1}
. Mobile porosity is
= 0.3 (**
-stag**) and immobile porosity
= 0.1. For the first-order exchange approximation in PHREEQC, a single cell immobile zone and the parameters
,
and
are specified with **
-stag**. This stagnant zone definition causes each cell in the mobile zone (numbered 1-20) to have an associated cell in the immobile zone (numbered 22-41). The **
PRINT** data block is used to eliminate all printing to the output file.

Following the pulse of NaCl solution, 10 shifts of 1 mmol KNO_{
3}
/L (second **
SOLUTION** 0) are introduced into the column. The second **
TRANSPORT** data block does not redefine any of the column or flow characteristics, but specifies that results for cells 1 through 20 (**
-punch_cells**) be written to the selected output file after 10 shifts (**
-punch_frequency**). The data blocks **
SELECTED_OUTPUT** and **
USER_PUNCH** specify the data to be written to the selected-output file.

Table 33.**
--**Input data set for example 13A: Stagnant zone with implicitly defined mixing factors

TITLE Example 13A.--1 mmol/l NaCl/NO3 enters column with stagnant zones. Implicit definition of first-order exchange model. SOLUTION 0 # 1 mmol/l NaCl units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 Na 1.0 # Na has Retardation = 2 Cl 1.0 # Cl has Retardation = 1, stagnant exchange N(5) 1.0 # NO3 is conservative # charge imbalance is no problem ... END SOLUTION 1-41 # Column with KNO3 units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0 EXCHANGE 1-41 -equil 1 X 1.e-3 EXCHANGE_SPECIES # For linear exchange, make KX exch. coeff. equal to NaX K+ + X- = KX log_k 0.0 -gamma 3.5 0.015 END TRANSPORT -cells 20 -shifts 5 -flow_d forward -timest 3600 -bcon flux flux -diffc 0.0 -length 0.1 -disp 0.015 -stag 1 6.8e-6 0.3 0.1 # 1 stagnant layer^, ^alpha, ^theta(m), ^theta(im) PRINT -reset false END SOLUTION 0 # Original solution reenters units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0 TRANSPORT -shifts 10 -punch_frequency 10 -punch_cells 1-20 SELECTED_OUTPUT -file ex13a.sel -reset false -solution -distance true USER_PUNCH -head Cl_mmol Na_mmol 10 PUNCH TOT("Cl")*1000, TOT("Na")*1000 END

The mixing factors *
mixf*
_{
m}
and *
mixf*
_{
im}
for the first-order exchange approximation for this example are derived from equations 121 and 123 as follows:

The retardation factors *
R*
_{
m}
and *
R*
_{
im}
are not included here in the formulas for *
mixf*
_{
im}
and*
mixf*
_{
m}
because in PHREEQC the retardation is a consequence of chemical reactions. According to equations 161 and 162, for this example the mixing factors are calculated to be *
mixf*
_{
im}
= 0.20886 and *
mixf*
_{
m}
= 0.06962. The mixing factors differ for the mobile cell and the immobile cell to account for the difference in the volume of mobile and immobile water.

In PHREEQC, a mixing of mobile and stagnant water is done after each diffusion/dispersion step. This means that the time step decreases when the cells are made smaller and when more diffusive steps ("mixruns") are performed. A 20-cell model as in example 13A has one mixrun. A 100 cell model would have 3 mixruns (equation 110 requires*
n*
=3 for *
mixf*
< 1/3), and the time step for calculating *
mixf*
_{
im}
would be (3600/5) / 3 = 240 s. A time step *
t*
= 240 s leads to *
mixf*
_{
im}
= 0.01614 in the 100 cell model.

The input file with explicit mixing factors for a uniform distribution of the stagnant zones is given in table 34. The **
SOLUTION** data blocks are identical to the previous input file. One stagnant layer without further information is defined (**
-stag** 1, in the **
TRANSPORT** data block). The mobile/immobile exchange is set by the mix fraction given in the **
MIX** data blocks. The results of this input file are identical with the results from the previous input file in which the shortcut notation was used. However, the explicit definition of mix factors illustrates that a non-uniform distribution of the stagnant zones, or other physical properties of the stagnant zone, can be included in PHREEQC simulations by varying the mixing fractions which define the exchange among mobile and immobile cells.

Table 34.**
--**Input data set for example 13B: Stagnant zone with explicitly defined mixing factors

TITLE Example 13B.--1 mmol/l NaCl/NO3 enters column with stagnant zones. Explicit definition of first-order exchange factors. SOLUTION 0 # 1 mmol/l NaCl units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 Na 1.0 # Na has Retardation = 2 Cl 1.0 # Cl has Retardation = 1, stagnant exchange N(5) 1.0 # NO3 is conservative # charge imbalance is no problem ... END SOLUTION 1-41 # Column with KNO3 units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0 EXCHANGE 1-41 -equil 1 X 1.e-3 EXCHANGE_SPECIES # For linear exchange, make KX exch. coeff. equal to NaX K+ + X- = KX log_k 0.0 -gamma 3.5 0.015 END MIX 1; 1 .93038; 22 .06962 ;MIX 2; 2 .93038; 23 .06962; MIX 3; 3 .93038; 24 .06962 ;MIX 4; 4 .93038; 25 .06962; MIX 5; 5 .93038; 26 .06962 ;MIX 6; 6 .93038; 27 .06962; MIX 7; 7 .93038; 28 .06962 ;MIX 8; 8 .93038; 29 .06962; MIX 9; 9 .93038; 30 .06962 ;MIX 10; 10 .93038; 31 .06962; MIX 11; 11 .93038; 32 .06962 ;MIX 12; 12 .93038; 33 .06962; MIX 13; 13 .93038; 34 .06962 ;MIX 14; 14 .93038; 35 .06962; MIX 15; 15 .93038; 36 .06962 ;MIX 16; 16 .93038; 37 .06962; MIX 17; 17 .93038; 38 .06962 ;MIX 18; 18 .93038; 39 .06962; MIX 19; 19 .93038; 40 .06962 ;MIX 20; 20 .93038; 41 .06962; # MIX 22; 1 .20886; 22 .79114 ;MIX 23; 2 .20886; 23 .79114; MIX 24; 3 .20886; 24 .79114 ;MIX 25; 4 .20886; 25 .79114; MIX 26; 5 .20886; 26 .79114 ;MIX 27; 6 .20886; 27 .79114; MIX 28; 7 .20886; 28 .79114 ;MIX 29; 8 .20886; 29 .79114; MIX 30; 9 .20886; 30 .79114 ;MIX 31; 10 .20886; 31 .79114; MIX 32; 11 .20886; 32 .79114 ;MIX 33; 12 .20886; 33 .79114; MIX 34; 13 .20886; 34 .79114 ;MIX 35; 14 .20886; 35 .79114; MIX 36; 15 .20886; 36 .79114 ;MIX 37; 16 .20886; 37 .79114; MIX 38; 17 .20886; 38 .79114 ;MIX 39; 18 .20886; 39 .79114; MIX 40; 19 .20886; 40 .79114 ;MIX 41; 20 .20886; 41 .79114; TRANSPORT -cells 20 -shifts 5 -flow_d forward -timest 3600 -bcon flux flux -diffc 0.0 -length 0.1 -disp 0.015 -stag 1 PRINT -reset false END SOLUTION 0 # Original solution reenters units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0 TRANSPORT -shifts 10 -punch_frequency 10 -punch_cells 1-20 SELECTED_OUTPUT -file ex13b.sel -reset false -distance true -solution USER_PUNCH -head Cl_mmol Na_mmol 10 PUNCH TOT("Cl")*1000, TOT("Na")*1000 END

The stagnant zone consists of spheres with radius *
r*
= 0.01 m. Diffusion into the spheres induces radially symmetric concentration changes according to the differential equation:

The calculation in finite differences can therefore be simplified to one (radial) dimension..

The calculation follows the theory outlined in the section "Transport in Dual Porosity Media" of this manual. The stagnant zone is divided into a number of layers that mix by diffusion. In this example, the sphere is cut in 5 equidistant layers with
*
r*
= 0.002 m. Five stagnant layers are defined under keyword **
TRANSPORT** with **
-stagnant** 5 (table 36). Mixing is specified among adjacent cells in the stagnant layers with **
MIX** data blocks; the mixing factors are calculated by equations 128 and 129. For the calculation, the volume *
V*
_{
j}
of cell *
j*
(m^{
3}
) is needed, the shared surface area *
A*
_{
ij}
of cell *
i*
and *
j*
(m^{
2}
), the distance between midpoints *
h*
_{
ij}
of cells *
i*
and *
j*
(m), and the correction factor for boundary cells *
f*
_{
bc}
(dimensionless). The values for mobile cell 1 and associated immobile cells are given in table 35. The cells in the immobile layer are numbered as *
n *
+ *
l *
x*
cells *
+ 1, where *
n*
is the number of a mobile cell, *
l *
is the number of the stagnant layer, and *
cells*
is the total number of mobile cells. In this example, the boundary cell in the stagnant zone for cell number 1 is cell number 22 and the other four stagnant layers are cell numbers 42, 62, 82, and 102, with number 102 being the innermost cell of the sphere, which is connected only to one other cell (cell 82). The volume of the mobile cell (cell 1) is expressed relative to the volume of a sphere of radius 0.01 m, by multiplying this volume with
/
(4.19e-6 ∞ 0.3 / 0.1 = 1.26e-5). In table 35 the value for *
f*
_{
bc}
is 1.72 as calculated from equation 127. It may be noted that using *
f*
_{
bc}
= 1.81 slightly improves the fit to an analytical solution given in Crank (1975) for diffusion into spheres in a closed vessel (a beaker with solution and clay beads). However, changing *
f*
_{
bc}
to 1.81 has little effect on the concentration profiles for the column which are shown in figure 13

Table 36.**
--**Input data set for example 13C: Stagnant zone with diffusion calculated by finite differences (partial listing)

TITLE Example 13C.--1 mmol/l NaCl/NO3 enters column with stagnant zones. 5 layer stagnant zone with finite differences. SOLUTION 0 # 1 mmol/l NaCl units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 Na 1.0 # Na has Retardation = 2 Cl 1.0 # Cl has Retardation = 1, stagnant exchange N(5) 1.0 # NO3 is conservative # charge imbalance is no problem ... END SOLUTION 1-121 # Column with KNO3 units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0 EXCHANGE 1-121 -equil 1 X 1.e-3 EXCHANGE_SPECIES # For linear exchange, make KX exch. coeff. equal to NaX K+ + X- = KX log_k 0.0 -gamma 3.5 0.015 END MIX 1; 1 0.90712; 22 0.09288 MIX 22; 1 0.57098; 22 0.21656; 42 0.21246 MIX 42; 22 0.35027; 42 0.45270; 62 0.19703 MIX 62; 42 0.38368; 62 0.44579; 82 0.17053 MIX 82; 62 0.46286; 82 0.42143; 102 0.11571 MIX 102; 82 0.81000; 102 0.19000 # # MIX definitions omitted for mobile cells # 2-19 and associated immobile cells # MIX 20; 20 0.90712; 41 0.09288 MIX 41; 20 0.57098; 41 0.21656; 61 0.21246 MIX 61; 41 0.35027; 61 0.45270; 81 0.19703 MIX 81; 61 0.38368; 81 0.44579; 101 0.17053 MIX 101; 81 0.46286; 101 0.42143; 121 0.11571 MIX 121; 101 0.81000; 121 0.19000 TRANSPORT -cells 20 -shifts 5 -flow_direction forward -timest 3600 -tempr 3.0 -bcond flux flux -diffc 0.0 -length 0.10 -disp 0.015 -stag 5 PRINT -reset false END SOLUTION 0 # Original solution reenters units mmol/l pH 7.0 pe 13.0 O2(g) -0.7 K 1.0 N(5) 1.0 TRANSPORT -shifts 10 -punch_frequency 10 -punch_cells 1-20 SELECTED_OUTPUT -file ex13c.sel -reset false -solution -distance true USER_PUNCH -head Cl_mmol Na_mmol 10 PUNCH TOT("Cl")*1000, TOT("Na")*1000 END

Note in table 36 that 121 solutions are defined, 1-20 for the mobile cells, and the rest for the immobile cells. The input file is identical with the previous one, except for **
-stag** 5 and the mixing factors among the cells. Not all the mixing factors are shown in table 36, they are identical for subsequent cells and their neighboring stagnant cells. In this example with clay beads, only radial (1D) diffusion is considered, and only mixing among cells in different layers is defined. However, it is possible to include mixing among the immobile cells of adjacent mobile cells.

Figure 13.**
--**Results of simulations of transport with diffusion into spherical stagnant zones modeled using finite difference and first-order exchange approximations.

Figure 13 compares the concentration profiles in the mobile cells obtained with example 13A (and B) with example 13C. The basic features of the two simulations are the same. The positions of the peaks, as calculated by the two simulations, are similar. The Cl^{
-}
peak is near 1.2 m, but would be at about 1.45 m in the absence of stagnant zones. The integrated concentrations in the mobile porosity are about equal for the first-order exchange and the finite difference simulations. The exchange factor
= 0.21 for the first-order exchange approximation appears to provide adequate accuracy for this simulation. However, the first-order exchange approximation produces lower peaks and more tailing than the more exact solution obtained with finite differences. Discrepancies can also appear as deviations in the breakthrough curve (Van Genuchten, 1985). The first-order exchange model is probably least accurate when applied to simulating the transport behavior of spheres; other shapes of the stagnant area can give a better correspondence. It is clear that the linear exchange model is much easier to apply because any explicit model requires the preparation of extended lists of mixing factors (notice that a separate simulation with **
USER_PUNCH** can serve that purpose), which change when the discretization is adapted. The calculation time for a finite difference model with multiple immobile-zone layers may also be considerably longer than for the single immobile-zone layer of the first-order exchange approximation.