[Next] [Previous] [Up] [Top]

EQUATIONS FOR SPECIATION AND FORWARD MODELING

## Equations for the Newton-Raphson Method

A series of functions, denoted by , are defined in this section. These functions describe heterogeneous equilibrium and are derived primarily by substituting the equations for the number of moles of species (derived from mass-action equations in the previous section) into mole- and charge-balance equations. Each function is presented along with the total derivative with respect to the master unknowns.

### Activity of Water

The activity of water is calculated from an approximation given by Garrels and Christ (1965, p. 65-66), which is based on Raoult's law:

(18) .

The function, , is defined as follows:

(19) ,

and the total derivative of this function is

(20) ,

The master unknown is the natural log of the activity of water.

### Ionic Strength

The ionic strength of the aqueous solution is a master unknown and is defined as follows:

(21) .

The function, , is defined as follows:

(22) ,

and the total derivative of this function is

(23) .

### Equations for Equilibrium with a Multicomponent Gas Phase

Equilibrium between a multicomponent gas phase and the aqueous phase is modeled with additional, heterogeneous mass-action equations. Only one gas phase can exist in equilibrium with the aqueous phase, but the gas phase may contain multiple components. The fugacity or activity of a gas component is assumed to be equal to its partial pressure. PHREEQC assumes the total pressure of the gas phase in equilibrium with a solution is fixed and is specified as Ptotal. If the sum of the partial pressures of the gas components in solution is less than Ptotal, the gas phase does not exist. The additional master unknown for the gas phase is the total number of moles of gas in the gas phase (including all gas components), Ngas. The number of moles of a gas component, g, in the gas phase is .

A mass-action equation is used to relate gas-component activities (fugacities) to aqueous phase activities. PHREEQC uses dissolution equations, in the sense that the gas component is assumed to be on the left-hand side of the chemical reaction. For carbon dioxide, the dissolution reaction may be written as follows:

(24) .

The Henry's law constant relates the partial pressure of the gas component to the activity of aqueous species. For carbon dioxide, the Henry's law constant is 10-1.468, and the following mass-action equation obtains at equilibrium:

(25) ,

where is the partial pressure calculated using activities in the aqueous phase. In general, the partial pressure of a gas component may be written in terms of aqueous phase activities as follows:

(26) ,

where is the partial pressure of gas component g, calculated using activities in the aqueous phase; is the Henry's law constant for the gas component; and is the stoichiometric coefficient of master species, m, in the dissolution equation. The values of may be positive or negative. For PHREEQC, terms on the left-hand side of a dissolution reaction are assigned negative coefficients and terms on the right-hand side are assigned positive coefficients.

At equilibrium, the number of moles of a gas component in the gas phase is equal to the partial pressure of the gas times the total number of moles of gas in the gas phase,

(27) .

The total derivative of the number of moles of a gas component in the gas phase is

(28) .

For mole-balance equations, the numerical model treats the gas phase components in the same way that it treats aqueous species. Thus, the terms appear in the Jacobian for the mole-balance equations for each element. The total number of moles of each element in the system includes both the number of moles in the gas phase and the number of moles in the aqueous phase.

Apart from the new terms in mole-balance equations, the one new function for the gas phase requires that the sum of the partial pressures of the component gases is equal to the total pressure, Ptotal. The function is defined as follows:

(29) .

The total derivative of with respect to the master unknowns, with the convention that positive dNgas are increases in solution concentration, is

(30) .

For data input to PHREEQC, the mass-action equations, Henry's law constant, and temperature dependence of the constant for gas phases are defined with the PHASES keyword data block. Components to include in gas-phase calculations and initial gas composition are defined with the GAS_PHASE keyword data block. (See Description of Data Input.)

### Equations for Equilibrium with Pure Phases

Equilibrium between the aqueous phase and pure phases, including single-component gas phases, is included in the model through the addition of heterogeneous mass-action equations. PHREEQC allows multiple pure phases, termed a pure-phase assemblage, to exist in equilibrium with the aqueous phase, subject to the limitations of the Gibbs' Phase Rule. The activity of a pure phase is assumed to be identically 1.0. The additional master unknown for each pure phase is the number of moles of the pure phase that is present in the system, np, where p refers to the pth phase. Terms representing the changes in the number of moles of each pure phase occur in the mole-balance equations for elements.

The new function corresponding to each of the new unknowns is a mass-action expression for each pure phase. PHREEQC uses dissolution reactions, in the sense that the pure phase is on the left-hand side of the chemical equation. For calcite, the dissolution reaction may be written as

(31) ,

and, using log K of 10-8.48 and activity of the pure solid is 1.0, the resulting mass-action expression is

(32) .

In general, pure-phase equilibria can be represented with the following equation:

(33) ,

where is the stoichiometric coefficient of master species, m, in the dissolution reaction. The values of may be positive or negative. For PHREEQC, terms on the left-hand side of a dissolution reaction are assigned negative coefficients and terms on the right-hand side are assigned positive coefficients. The saturation index for the mineral, SIp, is defined to be

(34) .

The function used for phase equilibrium in the numerical method is

(35) ,

where is a specified target saturation index for the phase (see keyword EQUILIBRIUM_PHASES) and converts base-10 log to natural log. For single-component gas phases, is equivalent to the log of the partial pressure of the gas. The total derivative with respect to the master unknowns is

(36) .

For data input to PHREEQC, the mass-action equations, equilibrium constant, and temperature dependence of the constant for pure phases are defined with the PHASES keyword data block. Initial composition of a pure-phase assemblage is defined with the EQUILIBRIUM_PHASES keyword data block. (See Description of Data Input.)

### Mole-Balance Equation for a Surface

Mole balance for a surface site is a special case of the general mole-balance equation. The total number of moles of a surface site is specified by input to the model. The sum of the moles of all of the surface species for the site must equal the total number of moles of surface sites. The following function is derived from the mole-balance relation for a surface site:

(37) ,

where the value of the function, fs, is zero when mole balance is achieved, Ts is the number of equivalents surface site s, and is the number of surface sites occupied by the surface complex. The total derivative of fs is

(38) .

For data input to PHREEQC, the number of moles of each type of surface site is defined with the SURFACE keyword data block. Surface species are defined with the SURFACE_SPECIES keyword data block. (See Description of Data Input.)

### Mole-Balance Equation for an Exchanger

Mole balance for an exchange site is a special case of the general mole-balance equation. The total number of moles of each exchange site is specified by input to the model. The sum of the moles of all of the exchange species for a site must equal the total number of moles of the exchange site. The following function is derived from the mole-balance relation for an exchange site:

(39) ,

where, the value of the function, fe, is zero when mole balance is achieved, Te is the total number of exchange sites for exchanger , and is the number of exchange sites occupied by the exchange species. The total derivative of fe is

(40) .

For data input to PHREEQC, the number of moles of exchange sites is defined in the EXCHANGE keyword data block. Exchange species are defined with the EXCHANGE_SPECIES data block. (See Description of Data Input.)

### Mole-Balance Equation for Alkalinity

The mole-balance equation for alkalinity is used only in speciation calculations and in inverse modeling. Mole balance for alkalinity is a special case of the general mole-balance equation, but special definitions of coefficients are needed. Alkalinity is defined as an element in PHREEQC and a master species is associated with this element (see SOLUTION_MASTER_SPECIES keyword). In the default databases for PHREEQC, the master species for alkalinity is . The master unknown for alkalinity is , or for the default databases, .

The total number of equivalents of alkalinity is specified by input to the model. The sum of the alkalinity contribution of each aqueous species must equal the total number of equivalents of alkalinity. The following function is derived from the alkalinity-balance equation:

(41) ,

where, the value of the function, fAlk, is zero when mole balance is achieved, TAlk is the number of equivalents of alkalinity in solution, and is alkalinity contribution of the aqueous species, eq/mol. The total derivative of fAlk is

(42) .

The value of may be positive or negative. Conceptually, a measured alkalinity differs from the alkalinity calculated by PHREEQC. In the default database files for PHREEQC the values of have been chosen such that the reference state ( ) for each element or element valence state is the predominant species at a pH of 4.5. It is assumed that all of the element or element valence state is converted to this predominant species in a theoretical alkalinity titration. However, in a real alkalinity titration, significant concentrations of species of elements and element valence states that have nonzero alkalinity contributions may exist at the endpoint of the titration, and the extent to which this occurs causes the alkalinity calculated by PHREEQC to be a different quantity than the measured alkalinity. Species that are especially susceptible to this problem are the hydroxide complexes of iron and aluminum. Thus, the alkalinity of a solution as calculated by PHREEQC, though it will be numerically equal to the measured alkalinity, is necessarily an approximation because of the assumption that a titration totally converts elements and element valence states to their reference state. In most solutions, where the alkalinity is derived predominantly from carbonate species, the approximation is valid.

For data input to PHREEQC, the alkalinity of each species is calculated from the association reaction for the species, which is defined in the SOLUTION_SPECIES keyword data block, and the alkalinity contributions of the master species, which are defined with the SOLUTION_MASTER_SPECIES keyword data block. Total alkalinity is part of the solution composition defined with the SOLUTION keyword data block. (See Description of Data Input.)

### Mole-Balance Equations for Elements

The total number of moles of an element in the system is the sum of the number of moles initially present in the pure-phase assemblage, aqueous phase, exchange assemblage, surface assemblage, gas phase, and diffuse layers of the surfaces. The following function is derived from the general mole-balance equation:

(43)

where the value of the function, fm, is zero when mole-balance is achieved. Tm is the total number of moles of the element in the system. E is the number of exchangers in the exchange assemblage, S is the number of surfaces in the surface assemblage. Np is the number of phases in the pure-phase assemblage, Naq is the number of aqueous species, Ne is the number of exchange species for exchanger e, Ns is the number of surface species for surface s, and Ng is the number of gas components. The number of moles of each entity in the system is represented by np for phases in the pure-phase assemblage, ni for aqueous species, for the exchange species of exchanger e, for surface species for surface s, ng for the gas components, and for the aqueous species in the diffuse layer of surface s. The number of moles of element, m, per mole of each entity is represented by bm, with an additional subscript to define the relevant entity; is usually, but not always, equal to (the coefficient of the master species for m in the mass-action equation), except for elements hydrogen and oxygen.

To avoid solving for small differences between large numbers, the quantity in parenthesis in the previous function is not explicitly included in the solution algorithm and the value of is never actually calculated. Instead the quantity is used in the function . Initially, is calculated from the total concentration of in the aqueous phase, the exchange assemblage, the surface assemblage, and the gas phase:

(44) .

During the iterative solution to the equations, is updated by the mole transfers of the pure phases:

(45) ,

where refers to the iteration number. It is possible for to be negative in intermediate iterations, but must be positive when equilibrium is attained.

The total derivative of the function, fm, is

(46)

For data input to PHREEQC, total moles of elements are defined initially through keyword data input and speciation, initial exchange, and initial surface calculations. Moles of elements are initially defined for an aqueous phase ( ) with the SOLUTION keyword data block, for an exchange assemblage ( ) with the EXCHANGE keyword data block, for a surface assemblage ( ) with the SURFACE keyword data block, for the gas phase ( ) with a GAS_PHASE keyword data block. The number of moles of each phase in a pure phase assemblage ( ) is defined with the EQUILIBRIUM_PHASES keyword data block. Total moles of elements and total moles of pure phases may be modified by reaction calculations. (See Description of Data Input.)

### Aqueous Charge-Balance Equation

The charge-balance equation sums the ionic charges of aqueous species and, in some cases, the charge imbalances developed on surfaces. For generality, net charge on an exchanger is also included in the derivation, though it is not justified by the theoretical framework. When specified, a charge-balance equation is used in initial solution calculations to adjust the pH or the activity of a master species (and consequently the total concentration of an element or element valence state) to produce electroneutrality in the solution. The charge-balance equation is used to calculate pH in reaction and transport simulations.

In real solutions, the sum of the equivalents of anions and cations must be zero. However, analytical errors and unanalyzed constituents in chemical analyses generally cause electrical imbalances to be calculated for solutions. If a charge imbalance is calculated for an initial solution, the pH is adjusted in subsequent reaction or transport simulations to maintain the same charge imbalance. If mixing is performed, the charge imbalance for the reaction step is the sum of the charge imbalances of each solution weighted by its mixing factor. If a surface is used in a simulation and the explicit diffuse-layer calculation is not specified, then the formation of charged surface species will result in a charged surface. Similarly, if exchange species are not electrically neutral (all exchange species in the default database are electrically neutral), the exchanger will accumulate a charge. These charge imbalances must be included in the charge-balance equation to calculate the correct pH in reaction and transport simulations. In general, the charge imbalance for a solution is calculated at the end of the initial solution calculation and at the end of each reaction and transport simulation with the following equation:

(47) ,

where zi is the charge on the aqueous species and is the charge imbalance for aqueous phase q. If charged surfaces or exchangers are not present, the charge imbalance for a solution at the end of a simulation will be the same as at the beginning of the simulation.

The charge imbalance on a surface is calculated at the end of the initial surface calculation and at the end of each reaction and transport simulation with the following equation:

(48) ,

where is the charge imbalance for the surface, and is the charge on the surface species i of surface s. If the composition of the diffuse layer is explicitly included in the calculation (-diffuse_layer in the SURFACE keyword data block), then each solution should be charge balanced using one of the charge balance options, and will equal to zero.

Normally, exchange species have no net charge, but for generality, this is not required. However, the activity of exchange species (the equivalent fraction) is not well defined if the sum of the charged species is not equal to the total number of equivalents of exchange sites (exchange capacity). If charged exchange species exist, then the charge imbalance on an exchanger is calculated at the end of the initial exchange calculation and at the end of each reaction and transport simulation with the following equation:

(49) ,

where is the charge imbalance for the exchanger, and is the charge on the exchange species i of exchanger e.

The charge imbalance for the system is defined at the beginning of each reaction or transport simulation with the following equation:

(50) ,

where is the charge imbalance for the system, Q is the number of aqueous phases that are mixed in the reaction or transport step, is the mixing fraction for aqueous phase q.

The charge-balance function is

(51) ,

where is zero when charge balance has been achieved and the double summation for surfaces is present only if the diffuse-layer composition is not explicitly calculated. The total derivative of is

(52) ,

and, again, the double summation for surfaces is present only if the diffuse-layer composition is not explicitly calculated.

For data input to PHREEQC, charge imbalance is defined by data input for SOLUTION, EXCHANGE, and SURFACE keyword data blocks combined with speciation, initial exchange, and initial surface calculations. The charge on a species is defined in the balanced chemical reaction that defines the species in SOLUTION_SPECIES, EXCHANGE_SPECIES, or SURFACE_SPECIES keyword data blocks. (See Description of Data Input.)

### Surface Charge-Potential Equation without Explicit Calculation of the Diffuse-Layer Composition

By default, PHREEQC uses the approach described by Dzombak and Morel (1990) to relate the charge accumulated on the surface with the potential at the surface, . The surface-charge density is the amount of charge per area of surface material, which can be calculated from the distribution of surface species as follows:

(53) ,

where is the charge density for surface s in coulombs per square meter (C/m2), F is the Faraday constant in coulomb per mole (96,485 C/mol), As is the specific area of the surface material (m2/g), and Ss is the mass of surface material (g). At 25oC, the surface-charge density is related to the electrical potential at the surface by the following equation involving the hyperbolic sine:

(54)

where is the valence of a symmetric electrolyte, is the ionic strength, F is the Faraday constant in kilojoules per volt-equivalent (kJ V-1 eq-1, which equals C/mol), is the potential at the surface in volts, R is the gas constant (8.314 J mol-1 oK-1), and T is in Kelvin. The following assumptions apply to equation 54: (1) Although strictly valid only at 25oC, the constant 0.1174 is used at all temperatures, and (2) the valence of the electrolyte is assumed to be 1. See the following sections, Surface Charge-Potential Equation with Explicit Calculation of the Diffuse-Layer Composition and Non-Electrostatic Surface-Complexation Modeling, for alternate formulations of surface-complexation modeling.

The charge-potential function is defined as follows:

(55) ,

and the total derivative of this function is

(56) .

For data input to PHREEQC, calculation without an explicit diffuse layer is the default. Specific surface area ( ) and mass of surface ( ) are defined in the SURFACE keyword data block. The charge on a surface species is defined in the balanced chemical reaction that defines the species in the SURFACE_SPECIES keyword data block. (See Description of Data Input.)

### Surface Charge-Balance Equation with Explicit Calculation of the Diffuse-Layer Composition

As an alternative to the previous model for the surface charge-potential relation, PHREEQC optionally will use the approach developed by Borkovec and Westall (1983). Their development solves the Poisson-Boltzmann equation to determine surface excesses of ions in the diffuse layer at the oxide-electrolyte interface. Throughout the derivation that follows, it is assumed that a volume of one liter (L) contains 1 kg of water.

The surface excess is defined to be

(57) ,

where is the surface excess in mol m-2, is the location of the outer Helmholtz plane, is concentration as a function of distance from the surface in mol m-3, and is the concentration in the bulk solution. The surface excess is related to concentration in the reference state of 1.0 kg of water,

(58) ,

where is the surface excess of aqueous species i in mol/kg water. This surface-excess concentration can be related to the concentration in the bulk solution by

(59) ,

where is a function of the potential at the surface and the concentrations and charges of all ions in the bulk solution:

(60) ,

where , is the value at the outer Helmholtz plane, and , is the dielectric constant for water, 78.5 (unitless), and is the dielectric permittivity of a vacuum, 8.854x10-12 C V-1 m-1. The value of at 25oC is 0.02931 [(L/mol)1/2 C m-2, where L is liters, mol is moles, C is coulombs, and m is meters]. The relation between the unknown (X) used by Borkovec and Westall (1983) and the master unknown used by PHREEQC is .

The development of Borkovec and Westall (1983) calculates only the total excess concentration in the diffuse layer of each aqueous species. A problem arises in reaction and transport modeling when a solution is removed from the surface, for example, in an advection simulation when the water in one cell advects into the next cell. In this case, the total number of moles that remain with the surface needs to be known. In PHREEQC, an arbitrary assumption is made that the diffuse layer is a specified thickness and that all of the surface excess resides in the diffuse layer. The total number of moles of an aqueous species in the diffuse layer is then the sum of the contribution from the surface excess plus the bulk solution in the diffuse layer:

(61) ,

where refers to the number of moles of aqueous species i that is present in the diffuse layer due to the contribution from the bulk solution, refers to the number of moles that is added to the diffuse layer due to the surface excess calculation, is the mass of water in the system excluding the diffuse layer, is the mass of water in the diffuse layer of surface s, and . The mass of water in the diffuse layer is calculated from the thickness of the diffuse layer and the surface area, assuming 1 L contains 1 kg water:

(62) ,

where is the thickness of the diffuse layer in meters.

The total derivative of the number of moles of an aqueous species in the diffuse layer is as follows:

(63)

where the second term is the partial derivative with respect to the master unknown for the potential at the surface, . The partial derivative, , is equal to the integrand from equation 60 evaluated at :

(64) ,

and the partial derivative of the function with respect to the master unknown is

(65) .

In the numerical method, it is computationally expensive to calculate the functions , so the same approach as Borkovec and Westall (1983) is used in PHREEQC to reduce the number of function evaluations. A new level of iterations is added when the diffuse layer is explicitly included in the calculations. The functions and their partial derivatives are explicitly evaluated once at the beginning of each of these diffuse-layer iterations. During the model iterations, which occur within the diffuse-layer iterations, the values of the functions are updated using the following equation:

(66) ,

where k refers to the model iteration number and is the value that is evaluated explicitly at the beginning of the diffuse-layer iteration. The model iterations end when the Newton-Raphson method has converged on a solution, however, convergence is based on the values of the functions that are estimates. Thus, diffuse-layer iterations must continue until the values of the functions have converged within specified tolerances, that is, the changes in the values of the functions are small between one diffuse-layer iteration and the next.

When explicitly calculating the composition of the diffuse layer, the function involving the sinh of the potential unknown (equation 54) is replaced with a charge-balance function that includes the surface charge and the diffuse-layer charge:

(67) ,

where, the value of the function, , is zero when charge balance is achieved. The total derivative of the function, , is

(68) .

For data input to PHREEQC, explicit calculation of the diffuse layer is invoked using the -diffuse_layer identifier in the SURFACE keyword data block. Specific surface area ( ) and mass of surface ( ) are also defined in the SURFACE keyword data block. The charge on a surface species or an aqueous species is defined in the balanced chemical reaction that defines the species in the SURFACE_SPECIES or SOLUTION_SPECIES keyword data block. (See Description of Data Input.)

### Non-Electrostatic Surface-Complexation Modeling

Davis and Kent (1990) describe a non-electrostatic surface-complexation model. In this model, the electrostatic term is ignored in the mass-action expressions for surface complexes. In addition, no surface charge balance or surface charge versus potential relation is used; only the mole-balance equation is included for each surface site.

For data input to PHREEQC, the non-electrostatic model for a surface is invoked by using the -no_edl identifier in the SURFACE keyword data block. (See Description of Data Input.)

Activity of Water
Ionic Strength
Equations for Equilibrium with a Multicomponent Gas Phase
Equations for Equilibrium with Pure Phases
Mole-Balance Equation for a Surface
Mole-Balance Equation for an Exchanger
Mole-Balance Equation for Alkalinity
Mole-Balance Equations for Elements
Aqueous Charge-Balance Equation
Surface Charge-Potential Equation without Explicit Calculation of the Diffuse-Layer Composition
Surface Charge-Balance Equation with Explicit Calculation of the Diffuse-Layer Composition
Non-Electrostatic Surface-Complexation Modeling

User's Guide to PHREEQC - 07 MAY 96
[Next] [Previous] [Up] [Top]