The aqueous geochemistry sub-model was developed to simulate aqueous speciation and solubility equilibrium control of pure-phase minerals using an approach similar to PHREEQC and as reported previousy in Salmon et al (in review). Each dissolved component, X, and mineral phase, PP, is a state variable within the model and is mass-conserved. In addition to external loading from inflows, the model includes a configurable dissolved sediment flux term for the dissolved species and a sedimentation term for the particulates. If we remove the hydrodynamic terms (e.g., inflows/outflows, mixing) for brevity, the differential equation is defined for the dissolved state variables as:
Aqueous speciation and solubility equilibrium control
The aqueous geochemistry model is conceptually similar to other equilbrium codes (e.g. ... ). The model defines geochemical reactions in terms of components. An aqueous species C is created as a product of reaction between components A and B, as shown in the following example:
aA+bB⟺cC
where a, b and c are the stoichiometric coefficients. The reaction shown above is able to proceed back and forward depending on the prevailing conditions of the aqueous solution. The equilibrium of the reaction is a function of the standard free energy, and from this the familiar mass-action expression is defined:
KC=[C]c[A]a[B]b
where KC is denoted the equilibrium constant. In any natural aquatic system, numerous components are present and there are therefore hundreds of reactions similar to Eq. above that will be close to equilibrium. The many reactions form a complex and interdependent set of simultaneous expressions which can be solved numerically, discussed next.
For each component, Xx, (generally components are equivalent to elements, except for the case of elements with multiple redox states, in which case each redox state is assigned as a unique component), a set of M components is defined, X∈{X1,X2,...,XM}. Combination of all the components into an aqueous solution results in numerous reactions taking place of the form of Eq. .., and we define the resultant set of aqueous species that result as: C∈{C1,C2...,CN} , where N is the total number produced.
The `activity' of an aqueous species is not equivalent to its molality due to the nonideality of aqueous solutions; a species activity, ˜C is related to its molality according to the activity coefficient, γ:
[Ci]=~Ci=γiCi
Numerous methods exist for estimating a species activity coefficient. The first is known as the Davies equation:
logγi=−c1z2i(√μ1+√μ−0.3μ)
and the other one used in the AED2 solution is the Debye-Huckel equation:
logγi=−c1z2i√μ1+c2c3√μ+c4μ
where zi is the ionic charge of the ith species, c1−4 are constants that depend upon properties such as the temperature of the solution. The activity of a species is dependent on the ionic strength, μ, a solution property defined as:
μ=12N∑iz2iCiWaq
where Waq is mass of water in the aqueous phase.
To numerically solve the system of equations, the total number of moles of any component is estimated by adding up the amount contained in each species:
Tj=N∑i=1aijCij=1..M
where ))aix is the stoichiometric coefficient of component x in species i. The activity of each species is calculated using the mass-action equations, which are defined more generically as:
˜Ci=KiM∏j=1˜Xaijjj=1..N
Taking the logarithm of both sides of Equation above leaves:
log˜Ci=logKi+M∑j=1aijlog˜Xjj=1..N
which can be defined in matrix notation as:
˜C∗=K∗+A˜X∗
where ∗ denotes a vector quantity. The solution is achieved by iterating using a Newton-Rhapson scheme, which aims to minimize the residual, R, between the estimated component molalities (T) and the known values:
Rx=N∑i=1aixCi−Txi=1..N
In matrix notation:
R=A−1C−T
The Newton-Rhapson technique iterates towards a solution by employing a derivative (termed Jacobian) matrix, J:
R=JΔX
which is solved for ΔX at each iteration. The new estimates for X are then updated according to:
Xn+1=Xn+ΔX
and the whole scheme proceeds until the solution has sufficiently converged, R<ϵ, where ϵ is a pre-defined convergence criterion (10−10 in this code). The Jacobian matrix is defined as:
J=dRxdXk
which can be expanded to:
dRxdXk=N∑i=1aikdCidXk−dTxdXk
using Eq. .... The derivatives are then defined according to:
dCidXk=aikCiXk
and
dTxdXk=0
which are the substituted into Eq. .. to leave:
dRxdXk=N∑i=1aikaikCiXk.
Therefore for a given total number of moles of each element in each computational cell, using this scheme the model solves for the activity of each aqueous species, ionic strength and pH. The charge balance variable (denoted ubalchg ) is also subject to advection and mixing as all other state variables, and importantly, an estimate of ubalchg must be provided at any inflow boundaries. If this is not done, then the charge balance will be compromised, and will manifest in a poor pH prediction.
The Newton-Rhapson scheme described above is implemented by making use of the Simplex numerical solver (Barrodale and Roberts, 1980; Parkhurst and Appelo, 1999). This is because AED2 also allows for simulation of pure phase (i.e. non-aqueous phase) minerals. When minerals can precipitate and dissolve, properties such as the ionic strength ( MU ), activity of water ( XH2O ) and charge balance ( CB ) may be non-conservative in the aqueous solution. These are therefore included as unknowns in the simplex solver, which uses the Newton-Rhapson technique to solve the following matrix:
(RPPRX1⋮RXNRMURH2ORCBRWaqXPP)=(\pdRPPln(˜X1)…\pdRPPln(˜XN)\pdRPPdμ\pdRX1ln(˜X1)⋮\pdRWaqXPP1)(dln(˜X1)RX1⋮dln(˜XN)dμdln(˜XH2O)dln(˜XH+)dln(Waq)dXPP)
The operation of the simplex solver is complicated and the reader is referred to Barrodale and Roberts (1980) for a detailed account. The aqueous species used in the simulations is based on those from Nordstrom et al. (1990), termed the
WATEQ4F database.
Redox transformations
...
Other geochemistry sources/sinks (optional)