EXAMPLES
This example repeats the inverse modeling calculations of the chemical evolution of spring-water compositions in the Sierra Nevada that are described in a classic paper by Garrels and Mackenzie (1967). The same example is described in the manual for the inverse-modeling program NETPATH (Plummer and others, 1991 and 1994). The example uses two spring-water compositions, one from an ephemeral spring, which is assumed to be less chemically evolved, and one from a perennial spring, which is assumed to be more chemically evolved. The differences in composition between the ephemeral and perennial spring are assumed to be due to reactions between the water and the minerals and gases it contacts. The object of inverse modeling in this example is to find sets of minerals and gases that, when reacted in appropriate amounts, quantitatively account for the differences in composition between the solutions.
The analytical data for the two springs are given below:
The chemical compositions of minerals and gases postulated to react by Garrels and Mackenzie (1967) are as follows:
The keyword INVERSE_MODELING is used to define all of the characteristics of the inverse-modeling calculations, including the solutions and phases to be used, the mole-balance equations to be included, the uncertainties to be used, whether all or only "minimal" models will be printed, and whether ranges of mole transfer that are consistent with the uncertainties will be calculated. A series of identifiers (sub-keywords preceded by a hyphen) are used to specify the characteristics of the inverse model. The input data set for this example is given in table 19.
TITLE Example 11.--Inverse modeling of Sierra springs
SOLUTION 1
-units mmol/l
pH 6.2
Si 0.273
Ca 0.078
Mg 0.029
Na 0.134
K 0.028
Alkalinity 0.328
S(6) 0.010
Cl 0.014
SOLUTION 2
-units mmol/l
pH 6.8
Si 0.41
Ca 0.26
Mg 0.071
Na 0.259
K 0.04
Alkalinity 0.895
S(6) 0.025
Cl 0.03
INVERSE_MODELING 1
-solutions 1 2
-uncertainty 0.025
-range
-phases
Halite
Gypsum
Kaolinite precip
Ca-montmorillonite precip
CO2(g)
Calcite
Chalcedony precip
Biotite dissolve
Plagioclase dissolve
-balance
Ca 0.05 0.025
PHASES
Halite
NaCl = Na+ + Cl-
log_k 0.0
Biotite
KMg3AlSi3O10(OH)2 + 6H+ + 4H2O = K+ + 3Mg+2 + Al(OH)4- + 3H4SiO4
log_k 0.0
Plagioclase
Na0.62Ca0.37Al1.38Si2.625O8 + 5.5 H+ + 2.5H2O = \
0.62Na+ + 0.37Ca+2 + 1.38Al+3 + 2.625H4SiO4
log_k 0.0
END
The identifier -solutions selects the solutions to be used by solution number. Two or more solution numbers must be listed after the identifier. If only two solution numbers are given, the second solution is assumed to evolve from the first solution. If more than two solution numbers are given, the last solution listed is assumed to evolve from a mixture of the preceding solutions. The solutions to be used in inverse modeling are defined in the same way as any solutions used in PHREEQC models. Usually the analytical data are entered in a SOLUTION keyword data block, but solutions defined by reaction calculation in the current or previous simulations may also be used. The -uncertainty identifier sets the default uncertainty for each analytical datum. In this example a fractional uncertainty of 0.025 (2.5 percent) is assumed for all of the analytical data except pH. By default, the uncertainty in pH is 0.05 unit. The uncertainty for pH and any datum for any of the solutions can be set explicitly to a fractional value or an absolute value (in moles; equivalents for alkalinity) using the -balances identifier.
The phases to be used in the inverse-modeling calculations are defined with the -phases identifier. In addition, this identifier can be used to specify any phases that only dissolve or only precipitate. In this example, kaolinite, montmorillonite, and chalcedony (SiO2) are required to precipitate only. This means that kaolinite will be precipitating (negative mole transfer) in any model that contains the phase kaolinite; likewise for montmorillonite and chalcedony. Similarly, biotite and plagioclase are required to dissolve (positive mole transfer) if they are present in an inverse model.
All of the phases used in inverse modeling must be defined in PHASES or EXCHANGE_SPECIES keyword data blocks, either in the database file or the input file. Thus, all phases defined in the default database file, phreeqc.dat or wateq4f.dat, are available for use in inverse modeling. Halite, biotite, and plagioclase are not in the default database file phreeqc.dat and so they are defined explicitly in the PHASES data block in the input data set. For simplicity, the log K's are set to zero for these phases, which does not affect inverse modeling; however, the saturation indices calculated for these phases will be spurious. The formula for plagioclase has been altered slightly from that in the previous table to achieve an exactly charge-balanced reaction. All phases used in inverse modeling must have a charge-balanced reaction. This requirement is due to the inclusion of a charge balance constraint for each solution. Each solution is adjusted to charge balance for each model by adjusting the concentrations of the elements within their uncertainty limits. (If a solution can not be adjusted to charge balance using the given uncertainties, the solution will be noted in the output and no models will be found.) Because all of the solutions are charge balanced in the modeling process, phases must also be charge balanced or they will not be included in any models. Note that the reaction for plagioclase (table 19) is on two lines, but the program interprets the two lines to be a single logical line because of the backslash "\" at the end of the first of these two lines.
The -range identifier indicates that, in addition to finding all of the inverse models, each model that is found will be subjected to additional calculations to determine the range of values that each mole transfer may have, within the constraints of the uncertainties.
By default, every inverse model includes mole-balance equations for every element in any of the phases included in -phases (except hydrogen and oxygen). If mole-balance equations are needed for elements not included in the phases, that is for elements with no source or sink (conservative mixing for example), the -balances identifier can be used to include those elements in the formulation of the inverse-modeling equations (see example 12). In addition, the -balances identifier can be used to specify uncertainties for an element in each solution. For demonstration purposes in the example, the uncertainty for calcium is set to 0.05 (5 percent) in solution 1 and 0.025 (2.5 percent) in solution 2. In addition to the mole-balance equations, the following equations are included for every inverse model: charge balance for each solution, electron balance for the system, and water balance for the system.
The unknowns in these equations include the mole transfers for each phase, the mole transfers of redox reactions, and the uncertainty unknowns for each element in each solution (excluding hydrogen and oxygen). An uncertainty unknown is included for each solution for alkalinity. Finally, an uncertainty unknown is included for pH for each solution containing carbon(+4).
Results for the two inverse models found in this example are shown in table 20. The results begin with a listing of three columns for each solution that is part of the model. All columns are values in mol/kg water. The first column contains the original analytical data for the solution. The second column contains any adjustments to the analytical data calculated for the model. These adjustments must be within the specified uncertainties. The third column contains the revised analytical data for the solution, which is equal to the original data plus any adjustment.
After the listing of the solutions, the relative fractions of each solution in the inverse model are printed. With only two solutions in the model, normally the fraction for each solution will be 1.0. If more than two solutions are included in the inverse model, normally the sum of the fractions of the solutions, excluding the last solution, will equal 1.0. The fractions are actually derived from a mole-balance on water, so if hydrated minerals consume or produce significant amounts of water or evaporation is modeled (see example 12), the numbers may not sum to 1.0. The second and third column for the block giving solution fractions are the minimum and maximum fractional values that can be attained within the specified uncertainties. These two columns are nonzero only if the -range identifier is used. In this example, all fractions are identically 1.0; the amount of water from gypsum dissolution is too small to affect the four significant figures of the mixing fractions.
The next block of data in the listing contains three columns describing the mole transfers for the phases. The first column contains the inverse model that is consistent with the adjustments printed in the listing of the solutions. In this example, the adjusted solution 1 plus the mole transfers in the first column exactly equal the adjusted solution 2. Mole transfers that are positive indicate dissolution; mole transfers that are negative indicate precipitation. (Note that mole transfers of phases in reaction calculations are relative to the phase, not relative to solution: positive values mean an increase in the phase; negative values mean a decrease in the phase.) The second and third columns of mole transfers are the minimum and maximum mole transfers of each phase that are consistent with the specified uncertainties. These two columns are nonzero only if the -range identifier is used. These minima and maxima are not independent, that is, obtaining a maximum mole transfer of one phase places very strong constraints on the mole transfers for the other phases. However, the output does not show any linkages between the maximum value for one phase with maximum or minimum values for other phases.
No redox reactions occurred in this inverse model. If they had, the number of moles transferred between valence states of each element would be printed under the heading "Redox mole transfers".
Next, the sum of each uncertainty unknown divided by its uncertainty (
, a standardized sum of residuals) and the maximum percentage adjustment to any element in any solution are printed; these two values apply to the model printed in the left-hand column. Finally, if no inverse model can be found with any proper subset of the phases, the statement "Model contains minimum number of phases" is printed.
After all models are printed, a short summary of the calculations is printed, which lists the number of models found, the number of minimal models found (models with a minimum number of phases), the number of infeasible sets of phases for which inverse models were attempted but failed, and the number of calls to the inequality equations solver, cl1 (calculation time is generally proportional to the number of calls to cl1).
The results of the example show that inverse models exist using the phases suggested by Garrels and Mackenzie (1967). The main reaction is dissolution of biotite, calcite, and plagioclase, which consumes carbon dioxide; kaolinite and montmorillonite or kaolinite and chalcedony precipitate. The results of Garrels and Mackenzie (1967) fall within the range of mole transfers calculated in the first model of PHREEQC for all phases except carbon dioxide. This discrepancy is due to the fact that Garrels and Mackenzie (1967) did not account for the dissolved carbon dioxide in the spring waters. Garrels and Mackenzie (1967) also ignored a small discrepancy in the mole balance for potassium. PHREEQC accounts for this discrepancy in the adjustments to the concentrations of the elements. PHREEQC shows that by altering the concentrations within the specified uncertainty (2.5 percent) an inverse model can be found. Without making the calculations with PHREEQC including uncertainties, it is not obvious whether the discrepancy in potassium is significant. The results of PHREEQC are concordant with the results of NETPATH, except that NETPATH also must ignore the discrepancy in the potassium mole balance.
-------------------------------------------
Beginning of inverse modeling calculations.
-------------------------------------------
Solution 1:
pH 6.200e+00 + 1.242e-02 = 6.212e+00
Al 0.000e+00 + 0.000e+00 = 0.000e+00
Alkalinity 3.280e-04 + 5.500e-06 = 3.335e-04
C(-4) 0.000e+00 + 0.000e+00 = 0.000e+00
C(4) 7.825e-04 + 0.000e+00 = 7.825e-04
Ca 7.800e-05 + -3.900e-06 = 7.410e-05
Cl 1.400e-05 + 0.000e+00 = 1.400e-05
H(0) 0.000e+00 + 0.000e+00 = 0.000e+00
K 2.800e-05 + -7.000e-07 = 2.730e-05
Mg 2.900e-05 + 0.000e+00 = 2.900e-05
Na 1.340e-04 + 0.000e+00 = 1.340e-04
O(0) 0.000e+00 + 0.000e+00 = 0.000e+00
S(-2) 0.000e+00 + 0.000e+00 = 0.000e+00
S(6) 1.000e-05 + 0.000e+00 = 1.000e-05
Si 2.730e-04 + 0.000e+00 = 2.730e-04
Solution 2:
pH 6.800e+00 + -3.399e-03 = 6.797e+00
Al 0.000e+00 + 0.000e+00 = 0.000e+00
Alkalinity 8.951e-04 + -1.796e-06 = 8.933e-04
C(-4) 0.000e+00 + 0.000e+00 = 0.000e+00
C(4) 1.199e-03 + 0.000e+00 = 1.199e-03
Ca 2.600e-04 + 6.501e-06 = 2.665e-04
Cl 3.000e-05 + 0.000e+00 = 3.000e-05
H(0) 0.000e+00 + 0.000e+00 = 0.000e+00
K 4.000e-05 + 1.000e-06 = 4.100e-05
Mg 7.101e-05 + -8.979e-07 = 7.011e-05
Na 2.590e-04 + 0.000e+00 = 2.590e-04
O(0) 0.000e+00 + 0.000e+00 = 0.000e+00
S(-2) 0.000e+00 + 0.000e+00 = 0.000e+00
S(6) 2.500e-05 + 0.000e+00 = 2.500e-05
Si 4.100e-04 + 0.000e+00 = 4.100e-04
Solution fractions: Minimum Maximum
Solution 1 1.000e+00 1.000e+00 1.000e+00
Solution 2 1.000e+00 1.000e+00 1.000e+00
Phase mole transfers: Minimum Maximum
Halite 1.600e-05 1.490e-05 1.710e-05 NaCl
Gypsum 1.500e-05 1.413e-05 1.588e-05 CaSO4:2H2O
Kaolinite -3.316e-05 -5.518e-05 -1.140e-05 Al2Si2O5(OH)4
Ca-Montmorillon -8.156e-05 -1.107e-04 -5.213e-05 Ca0.165Al2.33Si3.67O10(OH)2
CO2(g) 2.909e-04 2.346e-04 3.543e-04 CO2
Calcite 1.258e-04 1.028e-04 1.326e-04 CaCO3
Biotite 1.370e-05 1.317e-05 1.370e-05 KMg3AlSi3O10(OH)2
Plagioclase 1.758e-04 1.582e-04 1.935e-04 Na0.62Ca0.37Al1.38Si2.625O8
Redox mole transfers:
Sum of residuals: 5.573e+00
Maximum fractional error in element concentration: 5.000e-02
Model contains minimum number of phases.
===============================================================================
Solution 1:
pH 6.200e+00 + 1.242e-02 = 6.212e+00
Al 0.000e+00 + 0.000e+00 = 0.000e+00
Alkalinity 3.280e-04 + 5.500e-06 = 3.335e-04
C(-4) 0.000e+00 + 0.000e+00 = 0.000e+00
C(4) 7.825e-04 + 0.000e+00 = 7.825e-04
Ca 7.800e-05 + -3.900e-06 = 7.410e-05
Cl 1.400e-05 + 0.000e+00 = 1.400e-05
H(0) 0.000e+00 + 0.000e+00 = 0.000e+00
K 2.800e-05 + -7.000e-07 = 2.730e-05
Mg 2.900e-05 + 0.000e+00 = 2.900e-05
Na 1.340e-04 + 0.000e+00 = 1.340e-04
O(0) 0.000e+00 + 0.000e+00 = 0.000e+00
S(-2) 0.000e+00 + 0.000e+00 = 0.000e+00
S(6) 1.000e-05 + 0.000e+00 = 1.000e-05
Si 2.730e-04 + 0.000e+00 = 2.730e-04
Solution 2:
pH 6.800e+00 + -3.399e-03 = 6.797e+00
Al 0.000e+00 + 0.000e+00 = 0.000e+00
Alkalinity 8.951e-04 + -1.796e-06 = 8.933e-04
C(-4) 0.000e+00 + 0.000e+00 = 0.000e+00
C(4) 1.199e-03 + 0.000e+00 = 1.199e-03
Ca 2.600e-04 + 6.501e-06 = 2.665e-04
Cl 3.000e-05 + 0.000e+00 = 3.000e-05
H(0) 0.000e+00 + 0.000e+00 = 0.000e+00
K 4.000e-05 + 1.000e-06 = 4.100e-05
Mg 7.101e-05 + -8.980e-07 = 7.011e-05
Na 2.590e-04 + 0.000e+00 = 2.590e-04
O(0) 0.000e+00 + 0.000e+00 = 0.000e+00
S(-2) 0.000e+00 + 0.000e+00 = 0.000e+00
S(6) 2.500e-05 + 0.000e+00 = 2.500e-05
Si 4.100e-04 + 0.000e+00 = 4.100e-04
Solution fractions: Minimum Maximum
Solution 1 1.000e+00 1.000e+00 1.000e+00
Solution 2 1.000e+00 1.000e+00 1.000e+00
Phase mole transfers: Minimum Maximum
Halite 1.600e-05 1.490e-05 1.710e-05 NaCl
Gypsum 1.500e-05 1.413e-05 1.588e-05 CaSO4:2H2O
Kaolinite -1.282e-04 -1.403e-04 -1.159e-04 Al2Si2O5(OH)4
CO2(g) 3.044e-04 2.474e-04 3.684e-04 CO2
Calcite 1.124e-04 8.874e-05 1.198e-04 CaCO3
Chalcedony -1.093e-04 -1.483e-04 -6.986e-05 SiO2
Biotite 1.370e-05 1.317e-05 1.370e-05 KMg3AlSi3O10(OH)2
Plagioclase 1.758e-04 1.582e-04 1.935e-04 Na0.62Ca0.37Al1.38Si2.625O8
Redox mole transfers:
Sum of residuals: 5.573e+00
Maximum fractional error in element concentration: 5.000e-02
Model contains minimum number of phases.
===============================================================================
Summary of inverse modeling:
Number of models found: 2
Number of minimal models found: 2
Number of infeasible sets of phases saved: 20
Number of calls to cl1: 62
------------------
End of simulation.
------------------
------------------------------------
Reading input data for simulation 2.
------------------------------------
-----------
End of run.
-----------