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 Maple1, 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:

\begin{equation*} \gamma_{i} = 0.5 z_{i} z_{e} \log \sqrt{I\,} \end{equation*}

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:


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

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);

# 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;

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 Psolid, Msolid and Nsolid 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/litre2. 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 distsq3, 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.

Figures 1-4

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.

3distsq DSQ in Table 2.