Resource for Students and Scientists

William D. Scott and T. H. Venkitachalam

Division of Environmental Science, Murdoch University

Perth, Western Australia, 6150, AUSTRALIA

Scott et al. (1991) and Wrigley et al. (1992) have produced a computer model of the complex equilibria involved in the formation of Struvite. This paper uses the symbolic codes to calculate the amount of Struvite that is expected to precipitate in a simplest inorganic system. A weighted interpolation and a double quadratic fit are used to close the calculation.

Struvite precipitation; waste; symbolic manipulators; phosphorus; magnesium; nitrogen; slow release fertilizer; kidney stones;

Waste systems and kidney stones often contain significant amounts of
struvite ($\ce{MgNH4PO4 * 6H2O}$). Taylor *et al.* (1963) and
Webb (1990) present experimental data on the system at 25°C and 30°C,
respectively. The modelling of such systems is becoming clear. Linder and
Little (1986) produced computer models for precipitation in such systems,
including calcium oxalate monohydrate as well as dicalcium phosphate
dihydrate but were unsuccessful in modelling the struvite system. The
solubility of struvite was computed by Verplaetse *et al.* (1985) but
they only allowed for the 1:1:1 ratios of the constituents, without the full
ramifications of electroneutrality, mass balance and common ion effects.
The full system was first considered by Scott *et al.* (1991) and
Wrigley *et al.* (1992). The calculation allows total flexibility
in the complex ionic system that results when struvite precipitates from
aqueous solution. This article reports the specifics; the generation of
symbolic codes with the symbolic manipulator Maple^{1}, and the numerical
solution in Fortran. Examples of the design calculation are presented in
tabular and graphical form.

The equilibrium calculation assumes that:

(1) Species are limited to those that can be derived from water, magnesium,
ammonia, and phosphate. In the liquid, species included are $\ce{Mg^{2}+}$, $\ce{Mg(OH)+}$,
$\ce{NH3}$, $\ce{NH4+}$, $\ce{PO4^{3}-}$, $\ce{HPO4^{2}-}$, $\ce{H2PO4-}$, $\ce{H3PO4}$,
$\ce{MgNH4PO4}$, $\ce{MgHPO4}$, $\ce{MgPO4-}$, $\ce{MgH2PO4+}$, $\ce{NH4PO4^{2}-}$,
$\ce{(NH4)2PO4-}$, $\ce{NH4HPO4-}$, as well as $\ce{OH-}$, $\ce{H+}$, and $\ce{H2O}$
itself. There is an allowance for the presence of excess non-reactive positive or
negative ions, *e. g.*, $\ce{Na+}$ or $\ce{Cl-}$.

(2) The effect of ionic strength $I$ on activity coefficients is cared for by
an *ad hoc* form that allows the calculation of the activities of
individual ions:

where $z_{e} = \prod z_{j}^{y_{j}}$ is the equivalent charge of the
environmental field of ion $i$, $y_{j}=\dfrac {m_{j}}{\sum m_{j}}$ is the
fraction of charge $j$ to the ‘pool‘ of ions of charge opposite in sign
to ion $i$. $\prod $ is the overall product from ions that have the opposite
charge to $i$. $\sum m_{j}$ is the sum over all ions of charge $j$ (Wrigley
*et al.*,1992).

(3) Only one solid phase is formed with all solid species mutually soluble in one another. The mole fractions of these species are identical to their activities.

The nature of this algorithm evolves from the fact that two independent input variables are necessary to fully describe the system. This is apparent from Gibb’s Phase Rule for the case of a constant temperature and pressure, with four independently variable components and two phases. It is found that the solutions are polynomial functions to high powers in $\ce{H}$ and $\ce{Mg}$ (at least 11th power in $\ce{H}$); hence these two ion concentrations are appropriate independent variables, input to the algorithm.

The calculation proceeds using the mass balance and electroneutrality equations in an input file to Maple, as presented in Table 1:

TABLE 1 MAPLE INPUT FILE MNP := K1*Mg*NH4*PO4; #Struvite NH4 := K2*NH3*H #Ammonium # statement not included if NH4 is used as a primary calculating variable H3P := K3*H2P*H; #Phosphoric acid H2P := K4*HP*H; #Dihydrogen phosphate HP := K5*PO4*H; #Hydrogen phosphate MP8 := K6*Mg^3*PO4^2; #Magnesium phosphate with 8 H2O MP22 := K7*Mg^3*PO4^2; #Magnesium phosphate with 22 H2O MHP := K10*Mg*HP; #Magnesium hydrogen phosphate MHPd := K12*MHP; #dissolved Magnesium hydrogen phosphate MP := K13*Mg*PO4; #Magnesium phosphate ion MH2P := K14*Mg*H2P; #Magnesium dihydrogen phosphate ion #MNPd := K15*MNP; #dissolved Struvite #NP := K16*NH4*PO4; #associated ammonium phosphate ion #N2P := K17*NH4^2*PO4; #associated diammonium phosphate ion #NHP := K18*NH4*HP; #associated #NH4 := NLQD/(1+1/(H*K2)); #allows Ammonium ion to be replaced with NLQD # OH, MgOH and MgOH2 are defined later in the Fortran programme #The electroneutrality equation, Exions cares for added positive ions ffff:=Exions+H+NH4+2*Mg+MgOH+MH2P-OH-3*PO4-2*HP-H2P-MP; #ffff:=Exions+H+NH4+2*Mg+MgOH+MH2P-OH-3*PO4-2*HP-H2P-MP-2*NP-N2P-NHP; PO4 := solve(ffff,PO4); readlib(fortran): #Activates the Fortran conversion library fortran(PO4,optimized,filename=`PO4.for`); #produces fortran codes in file gggg := MNP+MP8+MP22+MHP+MgOH2-1; #Overall balance in solid phase ggg := normal(gggg); gg := op(1,ggg); #Need to check that this is the right term g:=collect(gg,NH3); a:=coeff(g,NH3,2); b:=coeff(g,NH3,1); c:=coeff(g,NH3,0); fortran(a,b,c,optimized,filename=`abc.for`); quit; # auxiliary and alternative expressions: #NH4 := K2*NH3*H; #Ammonium #OH := K11/H; #Water dissociation #MgOH := K8*Mg*OH; #Magnesium hydroxide ion #MgOH2 := K9*Mg*OH^2; #Magnesium hydroxide # for the fortran program depending on the primary calculating variables # Mass balance expressions: #NLQD :=NH4+NH3; #Nitrogen balance in the liquid #PLQD :=PO4+HP+H3P+H2P+MHPd+MP+MH2P; #Phosphorus balance in the liquid #MLQD :=Mg+MgOH+MHPd+MP+MH2P; #Magnesium balance in the liquid # collection of the quartic terms using NH4: #g:=collect(gg,NH4); a:=coeff(g,NH4,4); b:=coeff(g,NH4,3); #c:=coeff(g,NH4,2); d:=coeff(g,NH4,1); e:=coeff(g,NH4,0); #fnct := a*x^4 + b*x^3 + c*x^2 + d*x + e; #fortran(a,b,c,d,e,optimized,filename=`abcde.for`);

The details of the Maple programme depend on the
number of species that are recognised and the primary calculating
variable, here NH3 or NH4. Inclusion of 6 additional species
(see Wrigley *et al.*, 1992 or Wrigley, 1993) results in a quartic
expression (with 5 constants, {a,b,c,d,e}, see the commented statements)
and a Newton-Raphson trial and error solution. The version presented (up to the
quit statement) is the simplest; it was used to produce the symbolic codes
in the Fortran programme. Note that the 12 equilibria presented fit the
experimental 25°C data of Taylor *et al.* (1963) with a Coefficient of
Determination of 0.76 and a Standard Relative Error of 26%
(Wrigley *et al.*, 1992).

The equation $gg = 0$ is quadratic in NH3 (the moles of $\ce{NH3}$ or $\ce{NH4OH}$ per litre). The expressions for the constants, {a,b,c} are given on page 8 of the Fortran programme, as listed in the APPENDIX. Each constant is a complex function of H and Mg as well as the equilibrium constants and auxiliary, known concentrations such as OH, MgOH. The structure of the calculation evolves so as to assign algebraic values to appropriate, unknown species concentrations. The electroneutrality equation is used to calculate PO4 and the solid mass balance, to calculate NH3. Though the expressions for the constants are complex, Maple outputs them as Fortran codes. The codes may be copied by hand but the chance of an error is great. It is recommended that the Maple programme be run to generate them.

The design calculation evolves from the mass balance relations:

\begin{equation*} P = P_{lqd} + f \cdot P_{solid}\qquad M = M_{lqd} + f \cdot M_{solid}\qquad N = N_{lqd} + f \cdot N_{solid} \end{equation*}where *P*, *M* and *N* are, respectively, the total Phosphorus, Magnesium
and Nitrogen in the mixed liquid/solid system in gram-atoms/litre and
*P*_{solid}, *M*_{solid} and *N*_{solid} are, respectively, the ratios of
the gram-atoms of Phosphorus, Magnesium and Nitrogen in the solid to the total
gram atoms in the solid. The symbol *f* is the total gram atoms of the
three components in the solid per litre; the unknown that specifies how much
solid will precipitate from the solution. Noting that two arguments
(here *H* and *Mg*) are necessary to specify the concentrations in the
solution and in the solid, the *lqd* and *solid* concentrations in the above
relations are specified. In practice, however, *P*, *M* and *N* are
specified as inputs. The task is to use the above
relations to calculate *f* (or *pf*) and the required values of *H* (or *pH*)
*Mg* (or *pMg*). This ‘inversion’ is an interactive trial and
error approach.

The APPENDIX lists the Fortran program. It relies on approximately 1000 ‘data
points’ to initiate the interpolation/extrapolation using do loop 8.
The range of the loop should allow for all reasonable values of
the arguments, *H* and *Mg*. The heart of the calculation is subroutine
state which calculates some 21 concentrations. The main program
gathers these data in a matrix (a listing of records) that is augmented,
as necessary, to include ‘points’ in regions of interest.

First the total concentrations in the mixed liquor sample, (*P*, *M* and *N*),
are entered. This includes all solid
and liquid forms, in mg/litre^{2}.
It is not always possible to ‘tag’ the calculation to *P* and compute
$f = (P - P_{lqd})/P_{solid}$. Hence a select option is provided; flag variable *ichoose* sets the variable
for the calculation of *f* with (1, 2, and 3) corresponding to (*P*, *M* and *N*);
*ichoose* = 2 corresponds to $f = (M - M_{lqd})/M_{solid}$; *ichoose* = 3,
$f = (N - N_{lqd})/N_{solid}$, ‘tagging’ the calculation to nitrogen. Convergence
seems to be related to the dominance of a particular element in the domain of
interest.

Next, the excess amount of positive and negative non-reactive ions in the solution are supplied. The electrical sum, Exions, enters into the electroneutrality expression. The two, individual excess ion values affect the activity coefficients calculations through the ionic strength.

With values of *H* and *Mg*, p-values of the total magnesium, nitrogen and phosphorus
(*pM*, *pN* and *pP*) are produced. From these ‘points’, the square of the ‘distance’
from the required values is calculated as *distsq*^{3},
which is used as an inverse weighting factor for the *pf*, *pH* and *pMg*
values for each ‘point’. These estimates are lodged in vector *pint*
after statement 20. During the (interpolation) procedure the 6 points with the
smallest value of *distsq* are found and ranked.

From here an (extrapolation) from the chosen domain is allowed by fitting
the 6 closest points to a double quadratic function (see do loop 30), using
a matrix inversion scheme. The inverse interpolation technique simply uses
the errors in the two arguments (*M* and *N*, *P* and *N*, or *M* and *N*
depending on *ichoose*). After the calculation, several possible input values
of *pH* and *pMg* are available on the screen through format statement 36.
The next guess is made following statement 50. If the guess is good,
it may be added to the data set, and the procedure repeats. After
convergence, final solid calculations are completed, including the struvite
content of the solid in weight % (statement 133).

The calculation proceeds on a rather complex, rough surface and it requires skill to gain convergence. However, if the solution is covered by the initial ‘transects’ of the solution space, the 6 points collected will be reasonably close to the solution. If the DSQ value is larger than about 0.01, there is little likelihood of a solution. That is, the subroutine state contains checks to assure that the solution obeys the mass and electroneutrality criteria, though it is possible that the concentrations become too high for validity or that the iterations with activity coefficients calculate unreasonable values. Accepting the limitations of the codes and that the algorithm requires that there be a solid phase present, the lack of a solution implies that no solid exists.

At this point, four different estimates can be used with each cycle of
the calculation. One picks *pH* and *pMg* values that are close
to the nearest values. Another uses a weighted estimate (interpolated).
This gives an estimate if the answer is within the selected set and
it is only useful during a single cycle. A third is to use the value
from the quadratic fit (extrapolated); this is only effective when close to
the correct value and may be an over-estimate. The last is
an average of the last two estimates. In different nutrient
or excess non-reactive ion concentrations, different estimates were
found to be useful. Fast convergence occurs when close to a solution.

Table 2 illustrates the convergence for the case when approximately equal molar
quantities of *P*, *Mg* and *N* are added at 0.01 molar (*P* = 310 mg/litre,
*M* = 243 mg/litre, and *N* = 140 mg/litre), with no non-reactive ions:

Table 2 An Example of the Convergence of the Algorithm Guess | Calculated Values comment | pH pMg | pf pH pML pPL pNL pMg pP pM pN DSQ 1.998 8.600 2.198 2.186 2.164 4.600 2.000 2.004 2.008 0.00007029 initial point 8.667 4.663 1.990 8.667 2.197 2.190 2.150 4.663 2.000 2.002 1.994 0.00004689 guess intpol 8.635 4.630 1.994 8.635 2.195 2.188 2.160 4.630 2.000 2.001 2.003 0.00000713 guess the avg 8.673 4.667 1.989 8.673 2.196 2.190 2.151 4.667 2.000 2.000 1.994 0.00004258 guess the avg 8.647 4.642 1.992 8.647 2.195 2.189 2.157 4.642 2.000 2.001 2.000 0.00000101 guess the avg 8.648 4.642 1.992 8.648 2.194 2.189 2.158 4.642 2.000 2.000 2.000 0.00000003 guess the avg -- Msolid,Psolid,Nsolid 0.354 0.347 0.299 Total Solids -- 870.373 mg/litre 87.235% w/w struvite

Another practical example with large nutrient levels is presented in Figs. 1-4:
Here a solution with 500 mg/litre total *P* and 500 mg/litre total *N* was
subjected to additions of total *Mg* from 100 mg/litre to 1000 mg/litre
(Fig. 1 and Fig. 2). A solution with 500 mg/litre total *N* and 500 mg/litre
total *Mg*, to additions of total Phosphorus from 500 mg/litre to
1000 mg/litre (Fig. 3). A solution with all three components at 500 mg/litre
was subjected to additions of excess chloride (Fig. 4). Except for Fig. 4,
all the solutions have background concentrations of *Na* and *Cl* of
500 mg/litre or equivalent. The solutions in Fig. 4 have excess
chloride levels above 500 mg/litre as given on the abscissa.

These examples have purposely used large concentrations as precipitation of
solid at dilute concentrations is difficult, despite the presence of many
common ions. There seems to be a defined solubility. The graphs show the
formation of large amounts of solid as either the Magnesium levels or the
Phosphorus levels are raised – around 3000 mg/litre solids. Note that all
the solids contain large amounts of water of hydration – as many as 22 with
$\ce{Mg3(PO4)2}$. However, large amounts of excess chloride can
redissolve this solid (Fig. 4). This corresponds to pH reduction (as might
occur with addition of *HCl*, for example) and relatively small amounts of
struvite in the solid (30%). Phosphorus addition (again – conceptually as
phosphoric acid) has an opposite effect on the fraction of struvite in the
solid but the solid precipitate still increases (Fig. 3). The Magnesium
addition (perhaps as $\ce{Mg(OH)2}$) increases both the solid amount and the
amount of struvite in the solid (Fig. 1). Along with this precipitation, the
*pH* rises to high levels (about 9.5); the total magnesium in the liquid
approaches 10^{–5} molar or 0.24 mg/litre; the total phosphorus in the
liquid approaches 10^{–3} molar or 31 mg/litre though the nitrogen
levels in the liquid hardly change, remaining around 10^{–1.7} molar or
280 mg/litre. The additions envisioned would best be considered as single
cation or anion additions, with *OH* or *H* as a balance. Of course,
in nature, all combinations are possible within the restraints of mass balance
and electroneutrality.

Presented are design calculations to predict complex equilibrium calculations for waste water treatment, soil science, biology and the canning industry. It is possible to complete these calculations realistically with struvite and similar compounds. Of course, this is still a rather ideal model for its purpose; it does not include $\ce{CO2}$ chemistry, $\ce{Ca}$ compounds or organics.

Linder, P., and J. Little (1986). Prediction of computer modelling of the precipitation of stone forming solids from urine. *Inorg. Chim. Acta*, *3*, 137-45.

Scott, W. D., T. J. Wrigley and K. M. Webb (1991). A computer model of struvite solution chemistry. *Talanta*, *38*, 889-95.

Taylor, A. W., A. W. Frazier and E. L. Gurney (1953). *Trans. Faraday Soc.*, *59*, 1580-84.

Verplaetse, H., R. M. H. Verbeeck, H. Minnaert, and W. Oosterlinck (1985). Solubility of inorganic kidney stone components in the presence of acid-base sensitive complexing agents. *Eur. Urol.*, *11*, 44-51.

Webb, K. M. (1988). *The solubility of struvite and its application to a piggery effluent problem.* Honours Thesis, Environmental Science, Murdoch University.

Wrigley, T. J. (1993). *Water quality improvement of piggery effluent*, Ph.D. Thesis, Murdoch University.

Wrigley, T. J., W. D. Scott and K. M. Webb (1992). An improved computer model of struvite solution chemistry. *Talanta*, *39*, 1597-603.

The Fortran programme strevl.f produced the above results (see Scott *et al.*,1991 and Wrigley *et al.*, 1992).

1 …Maple Maple is a registered trademark of Waterloo Maple Software, Waterloo, Ontario

2 …mg/litre The programme uses molar, mole fractions or gram-atoms and p-values.

3 …*distsq* DSQ in Table 2.