Equilibrium equations and solution method
Notation and definitions are described in Nomenclature.
This page follows Bethke (2007).
Notation
Introduce the following:
: water.
: aqueous basis species. This is a label, not a quantity. for example, one of the might be Na.
: other aqueous species, which are called the "secondary" species.
: minerals in the equilibrium system of reactions that appear in the basis.
: all minerals, even those that do not exist in the equilibrium system. This might include minerals whose reactions are kinetically controlled.
: gases of known fugacity that appear in the basis.
: all gases.
: sorbed species. These do not get transported — they are surface complexes.
: sorbing sites. In the Langmuir approch and ion-exchange approach there is just one of these, and all sorbing species compete to sorb. In the "surface complexation" approach, there are an equal number of as .
Basis
Note that the index , , etc, gives meaning to the symbol that it is indexing. Although unusual in mathematics literature, this is the common convention in geochemistry. There are unknowns and the basis is Being a basis, all other other components can be expressed in terms of it.
Secondary species
All secondary species are in equilibrium with the basis. They are not sorbing, so This equation is specified in the database with the being the stoichiometric reaction coefficients. The database specifies the equilibrium constant, , and mass-action equilibrium reads Here:
is the activity of water
is the activity of the primary species
is the activity of mineral , which is assumed to be unity but is written here for clarity
is the known gas fugacity
is the molality of the secondary species
is the activity coefficient of the secondary species.
denotes molality
As discussed on the activity and fugacity pages, geochemists typically ignore the dimensional inconsistencies in this equation. Instead, it is conventional to omit dimension-full factors of "1", and simply use consistent units of moles, kg, bars and Kelvin throughout all calculations.
The above equation may be rearranged for the molality of the secondary species
Minerals, activity product and saturation index
For any mineral (including those in the basis, but the equation is trivial for them) Define the activity product which may be compared with the expression for the reaction's equilibrium constant, . Define the "saturation index": There are three possibilities:
. The mineral is at equilibrium.
. The mineral is supersaturated, which means the system will try to precipitate (perhaps slowly) the mineral.
. The mineral is undersaturated, which means the aqueous solution will (perhaps slowly) dissolve more of the mineral.
Sorption
and Feundlich approaches
Here, the basis does not include any . The equations are summarised in Chapter 9 of Bethke (2007). However, Bethke recommends against using these approaches, since adsorption may occur without limit.
Langmuir approach
Here a single sorption site is introduced into the basis () and the reaction for each sorbed species, , written in terms of the basis is Here is generally unity, but is written here for consistency of notation. Each such reaction has an associated equilibrium constant, , and a mass action that reads Here is the molality of unoccupied sites and is the molality of the secondary species.
Ion exchange
Ion exchange does not treat the sorption of a species on the surface of the porous skeleton, but the replacement of one adsorped ion with another. The development of the theory is similar to the Langmuir approach of sorption however. A new basis species representing exchange sites, is introduced. The formula for is similar to the Langmuir equation above, but involves the electric charge of the ions involved in the exchange as well as multiplicative factors.
The ion-exchange approach is not currently supported in the geochemistry module.
Surface complexation approach
Here, the basis includes an entry, for each type of surface site considered, and there is a reaction to form each surface complex : The mass action for each surface complex is (1) which sets the molality of each surface complex. Here:
[dimensionless] is the electronic charge of the surface complex
[J.C = V] is the surface potential, which is set as below
C.mol is the Faraday constant
J.K.mol is the gas constant
[K] is temperature
Computing the surface potential, , requires an additional equation, which Bethke (2007) writes as (2) Note that the right-hand side of this equation depends on through Eq. (1). The additional notation introduced here is as follows.
[m] is the surface area of material upon which the live. Usually a mineral in the system provides this surface area and usually the user specifies a specific surface area in [m/g(of mineral)]. Hence is proportional to the mole number of a mineral, , with coefficient of proportionality being the mineral's density [g.mol] and the specific surface area.
[dimensionless] is the dielectric constant: at 25C for water.
F.m is the permittivity of free space.
kg.m is the density of water.
is the ionic strength.
is the charge on the background solute. Bethke (2007) assumes
[dimensionless] is the charge of the surface complex .
Equations relating molality and mole numbers of the basis
Equations may be formed to relate the molality and mole number of the basis species to each other and the molality of the secondary species.
Let be the mass [kg] of solvent water. The secondary species also contain water: one mole of secondary species contains moles of water "inside it". The same holds for sorbed species. Therefore, so the total number of moles of water is Here is the number of moles of HO in a kg of water.
Similarly, in terms of the molality , the mole-number (total number of moles of substance) of basis species is There are of these equations.
Let be the mole number of mineral . (Note that is a mass, while has dimensions of moles: I follow the notation of Bethke (2007).) This may be unknown, since an experiment might commence by adding a certain amount of mineral, but the reactions might produce or consume some of it. Nevertheless, the total number of moles of basis mineral is There are of these equations.
It is assumed that the gas is buffered to an infinite reservoir, so for equilibrium reactions we only ned to consider the gas components that exist in the aqueous solution as part of the secondary species: There are of these equations.
The total mole number of sorbing sites, , may be written in terms of the molality of the sorbing sites, and the molality of the sorbed species, :
Full equations, assumptions, and reduced equations
In summary, the full set of equations are In addition, if surface complexation is present, Eq. (2) provides an equation for the surface potential, . In the above equations:
[mol] is the total number of moles of a substance
[mol.kg] is the molality
[kg] is the mass of solvent water
[mol] is the mole number of mineral
[dimensionless] is a stoichiometric coefficient
is an activity
is an activity coefficient
is a gas fugacity
is an equilibrium constant
is unity for the Langmuir approach to sorption, involves electonic charge for ion exchange, and involves the surface potential for the surface complexation approach
The last two equations are not dimensionally consistent: see the note above.
These equations relate quantities, as well as the addition equation for if relevant. Once the basis species are known, the molalities of the secondary and any sorbed species may be worked out simply using the final two equations.
The known conditions for the reactions can take a varity of forms:
For HO: either (total number of moles) or (number of moles acting as the solute) is known.
For each primary species: either (bulk-composition number of moles) or (molality of free primary species) is known. E.g. specifying pH specifies . Charge balance must hold.
For each mineral: either (bulk-composition number of moles) or (number of moles of basis mineral ) is known
Because an infinite buffer reservoir of gas is assumed, the initial concentration of gas is irrelevant — the user simply specifies the fixed fugacity
For each sorption site: either or is known.
To solve these equations, the following assumptions are made.
The mineral's activities are all unity. This means the secondary species' molalities and sorbed species' molalities do not depend on the mineral activities. This means that, given the secondary species' molalities and sorbed species' molalities, the equation may be used to provide in terms of known , or in terms of known , depending on the problem setup. Therefore, this equation is trivial and may be ignored during the solution process. It is simply used once the solution is found.
The gas fugacities are known and fixed. For the same reason as the previous bullet point, this means that the equation may be ignored during the solution process. It is simply used to provide (the total number of moles of gas in the aqueous solution) once the solution is found.
The activity coefficients are known functions.
If surface complexation is employed, is fixed (see Eq. (2)). If depends on , as is often the case, this breaks the "splitting into reduced equations" just described. However, the only examples I have found are when is fixed by specifying the amount of free mineral.
This means that the "reduced" equations are as well as Eq. (2) if appropriate. The problem becomes: given , find , and (and , if appropriate). The Newton-Raphson (NR) procedure (below) is used. If is known, the first equation is trivial, and needn't appear in the NR. Similarly, if any or are known (e.g. pH is fixed by the user) then that equation needn't appear in the NR.
Solution method
The NR procedure is used to find the in terms of the . Bethke (2007) describes some detail and some useful procedures that help convergence.
Initial configuration
The initial configuration for the is whichever is the most appropriate of:
set by the user;
90% of the specified values for ;
the result of the previous time step;
if the residual for basis species is extremely large, then the starting are repeatedly halved until the residual reaches a manageable size of less than .
Activity
At every step in the NR procedure, activity coefficients and as well as the activity of water, , are computed based on the current values of and . The derivatives of the activity coefficients with respect to are not included in the NR Jacobian.
Under-relaxation
Because the quantities can never be negative, the NR iteration is under-relaxed. At iteration number , define then the NR update takes the form
Charge balance
If a basis-species free molality is specified (such as ) then charge balance cannot be assured until NR has converged. To force electrical neutrality at the end of each NR iteration, the equation (where is the charge of basis species ) is used to set the bulk concentration of one of the basis species. That is, one of the equations is replaced by this charge-balance equation, to set the bulk concentration, , of a basis species in abundant concentration for which the greatest analytic uncertainty exists. This is generally Cl, but can be set by the user. See Section 4.3 and 14.2 of Bethke (2007).
Mineral undersaturation and supersaturation
After NR convergence, and substituting the result to find and , it may be found that , meaning the mineral was completely consumed (and more). Similarly, computing the saturation index of minerals not in the basis may reveal a supersaturated mineral that should be allowed to precipitate. Therefore, after NR convergence, the following procedure involving basis swaps is used.
All are computed. If one or more a negative, then the one that is the most negative is removed from the basis. It is replaced by secondary species that satisfies . The NR procedure is resolved.
The sauration index of each mineral, , that can form in a reaction model is checked for supersaturation. If one or more are supersaturated, the one with the largest saturation index is added into the basis, and something removed from the basis. This "something is": preferably, a primary species satisfying ; if no primary species appears in the reaction for then the mineral in the reaction for that satisfies . The NR procedure is re-solved (and checked for undersaturation and supersaturation, etc). Finally, due to a multitute of complexities, this procedure involving supersaturation is undesirable — the supersaturated mineral is ignored because the modeller knows it doesn't occur in reality.
Primary species with low concentration
Sometimes a species in the basis can occur at very small concentration, leading to numerical instability because making small corrections to its molality leads to large deviations in the molalities of the secondary species. In this case, it may be swapped for an abundant secondary species. Bethke (2007) doesn't give conditions for when this type of swap should occur, or whether it's supposed to be an automatic or user-controlled process
Final check
After NR convergence, and all the other checks above have been satisfied, the following checks are made:
the masses of each component are conserved
the activity product is equal to the equilibrium constant for each reaction
References
- Craig M. Bethke.
Geochemical and Biogeochemical Reaction Modeling.
Cambridge University Press, 2 edition, 2007.
doi:10.1017/CBO9780511619670.[BibTeX]