The following example demonstrates the capability of PHREEQC to calculate transient transport of heat and solutes in a column or along a 1D flowline. A column is initially filled with a dilute KCl solution at 25 °C in equilibrium with a cation exchanger. A KNO_{
3}
solution then advects into the column and establishes a new temperature of 0 °C. Subsequently, a sodium chloride solution at 24 °C is allowed to diffuse from both ends of the column, assuming no heat is lost through the column walls. At one end, a constant boundary condition is imposed, and at the other end, the final cell is filled with the sodium chloride solution and a closed boundary condition is prescribed. For the column end with a constant boundary condition, an analytical solution is compared with PHREEQC results, for unretarded Cl^{
-}
*
(R*
= 1.0) and retarded Na^{
+}
and temperature (*
R*
= 3.0). Finally, the second-order accuracy of the numerical method is verified by increasing the number of cells by a factor of three and demonstrating a decrease in the error of the numerical solution by approximately one order of magnitude relative to the analytical solution.

The input file for example 12 is shown in table 31. The EXCHANGE_SPECIES data block is used (1) to make the exchange constant for KX equal to NaX (log_k = 0.0), (2) to effectively remove the possibility of hydrogen ion exchange, and (3) to set activity coefficients for exchange species equal to their aqueous counterparts (**
-gamma**
identifier), so that the exchange between Na^{
+}
and K^{
+}
is linear and the retardation is constant. The influent, SOLUTION 0, is defined with temperature 0 °C and 24 mmol/kgw KNO_{
3}
. All solutions are defined to be in equilibrium with atmospheric oxygen partial pressure. The column is discretized in 60 cells, which are filled initially with a 1-μmol/kgw KCl solution (SOLUTION 1-60). Each cell has 48 mmol of exchange sites, which are defined to contain only potassium by the data block EXCHANGE 1-60.

The TRANSPORT data block defines cell lengths of 1 m (**
-lengths **
1), no dispersion (**
-dispersivities**
0), no diffusion (**
-diffusion_coefficient**
0), and no retardation for temperature (**
-thermal_diffusion**
1). SOLUTION 0 is shifted 60 times into the column (**
-shifts**
60, **
-flow_direction**
forward), and arrives in cell 60 at the last shift of the advective-dispersive transport simulation. The boundary conditions at the column ends are of flux type (**
-boundary_conditions**
flux flux). In this initial advective-dispersive transport simulation, no exchange occurs because the exchange sites are already completely filled with potassium, and the concentrations and temperatures in the 60 cells evolve to 24 mmol/kgw KNO_{
3}
and 0 °C. This result could be more easily achieved with SOLUTION data blocks directly, but the simulation demonstrates how transient boundary and flow conditions can be represented. The dissolved and solid compositions and the temperature of each cell of the column is automatically saved after each shift in the advective-dispersive transport simulation. The keyword PRINT is used to exclude all printing to the output file (**
-reset**
false).

In the next advective-dispersive transport simulation, diffusion is calculated from the column ends. The column composition and temperatures are initially the conditions produced by the first advective-dispersive transport calculation, except that an NaCl solution is now defined as solution 0, which diffuses into the top of the column, and as solution 60, which diffuses from the bottom of the column. The new SOLUTION 0 is defined with a temperature of 24 °C and with 24 mmol/kgw NaCl. The last cell (SOLUTION 60) also is defined to have this solution composition and temperature. The exchanger in cell 60 is defined to be in equilibrium with the new solution composition in cell 60 (EXCHANGE 60).

The TRANSPORT data block defines one diffusive transport period (**
-shifts**
1; **
-flow_direction**
diffusion). The boundary condition at the first cell is constant concentration, and at the last cell the column is closed (**
-boundary_conditions**
constant closed). The effective diffusion coefficient (**
-diffusion_coefficient**
) is set to 0.3×10^{
-9}
m^{
2}
/s, and the time step (**
-time_step**
) is defined to be 1.0×10^{
10}
s. Because only one diffusive time period is defined (**
-shifts**
1), the total time modeled is equal to the time step, 1.0×10^{
10}
s. However, the time step will automatically be divided into a number of time substeps to satisfy stability criteria for the numerical method. The heat retardation factor is set to 3.0 (**
-thermal_diffusion**
3.0). For Na^{
+}
the ratio of exchangeable concentration (maximum NaX is 48 mmol/kgw) to solute concentration (maximum Na^{
+}
= 24 mmol/kgw) is 2.0 for all concentrations, and the retardation is therefore *
R*
= 1 + *
d*
NaX/*
d*
Na^{
+}
= 3.0, which is numerically equal to the temperature retardation.

The SELECTED_OUTPUT data block specifies the name of the selected-output file to be *
ex12.sel*
. The identifier “**
-high_precision **
true” is used to obtain an increased number of digits in the printing, and the identifiers **
-dist**
and **
-temp **
specify that the distance and temperature of each cell will be printed to the file. The USER_PUNCH data block is used to print concentrations of sodium, potassium, and chloride to the selected-output file in units of mmol/kgw, and USER_GRAPH plots the data.

In the model, the temperature is calculated with the (linear) retardation formula; however, the Na^{
+}
concentration is calculated by the cation-exchange reactions. Even though the Na^{
+}
concentration and the temperature are calculated by different methods, the numerical values should be the same because the initial and the transient conditions are numerically equal and the retardation factors are the same. The temperature and the Na^{
+}
concentration are equal to at least 6 digits in the PHREEQC selected-output file, which indicates that the algorithm for the chemical transport calculations is correct for the simplified chemistry considered in this example. A further check on the accuracy is obtained by comparing simulation results with an analytical solution. For an infinite column with *
C*
_{
x}
= 0 for *
t*
= 0, and diffusion from *
x*
= 0 with *
C*
_{
x=0}
= *
C*
_{
0}
for *
t*
> 0 the analytical solution is

where De is the effective diffusion coefficient and R is the retardation factor.

The PHREEQC results are compared with the analytical solution for Cl^{
-}
and for temperature and Na^{
+}
in figure 13 and show excellent agreement. Notice that diffusion of Cl^{
-}
from the column ends has not yet “touched” in the mid-section, so that the column is still effectively infinite and the analytical solution is appropriate. Although both ends of the column started with the same temperature and concentration, at*
x*
= 0 m, the same temperature and concentrations are maintained because of the constant boundary condition. At *
x*
= 20 m (plotted at the midpoint of cell 60,*
x*
= 19.83 m), the temperature and concentrations decrease because of the closed boundary condition; no flux of heat or mass through this boundary is allowed, and the temperature and concentrations are diminishing because of diffusion into the column. The sodium concentration is not dissipating as rapidly as the chloride concentration because exchange sites must be filled with sodium along the diffusive reach.

Because this example has an analytical solution, it is possible to verify the second-order accuracy of the numerical algorithms used in PHREEQC. For a second-order method, decreasing the cell size by a factor of three should improve the results by about a factor of nine. The deviations from the analytical solution at the end of the time step are calculated at distances from 0.5 to 8.5 m in 0.5 m increments. The results are shown in table 32. As expected for a second-order method, the deviations from the analytical solution decreased by approximately an order of magnitude when the cell size was decreased by a factor of three.