Thermodynamic Modeling of Phosphorus Recovery from Wastewater

Municipal wastewater is an interesting source of phosphorus and several processes for the recovery of phosphorus from this source have been described. These processes yield magnesium ammonium phosphate (MAP), a valuable fertilizer. In these processes, pH shifts and the addition of chemicals are used to influence the species distribution in the solution such as to finally obtain the desired product and to prevent the co-precipitation of salts of heavy metal ions. Elucidating these species distributions experimentally is a challenging and cumbersome task. Therefore, in the present work, a thermodynamic model was developed that can be used for predicting the species distributions in the various steps of the recovery process. The model combines the extended Debye-Hückel equation for the prediction of activity coefficients with dissociation constants and solubility product data from the literature and contains no parameters that need to be adjusted to process data. The model was successfully tested by comparison to experimental data for the Stuttgart process from the literature and used for analyzing the different process steps. Furthermore, it was demonstrated how the model can be used for optimizing the process.


Introduction
Phosphorus is a key element for every living organism. Hence, phosphorus-containing substances with a high plant availability are used as agricultural fertilizers to enhance the cultivation of plants [1]. To date, the main source for these substances is rock phosphorus. However, similar to Hubbert's peak oil predictions from 1949 [2], recent simulations indicate that rock phosphorus reserves could be exhausted within the next 100 years [3]. Making other phosphorus sources available is thus highly desirable.
The simultaneous recovery of the two nutrients phosphorus (P) and nitrogen (N) by crystallization of magnesium ammonium phosphate (MAP, MgNH 4 PO 4 ·6 H 2 O) from wastewater or sewage sludge is an attractive alternative to rock phosphorus mining [4][5][6][7]. MAP, once an undesired by-product in wastewater treatment plants [8,9], is nowadays used as a fertilizer [10]. Several processes were developed to produce MAP from wastewater [6,11], including the Stuttgart process [5], which is investigated in the present work. Although the recovery of phosphorus as MAP using digested sludge from plants with chemical P-removal requires more efforts than that from plants with enhanced biological P-removal (EBPR), chemical P-removal is the main operation mode in Germany where peak effluent values of 1 mg L −1 need to be undercut at larger waste water treatment plants and much lower values are discussed even for smaller plants.
Thus, digested sewage sludge from plants with primary sedimentation and simultaneous chemical phosphorus elimination via iron-salts serves as the feed of the Stuttgart process. Essentially, digested sewage sludge is an aqueous solution with many dissolved species loaded with a solid phase [12,13]. Hence, it is a complex multi-component, multi-phase system. Although its actual composition varies with time and location, the aqueous phase typically contains high amounts of the nonmetal ions NH 4 + and PO 4 3− [13], the metal cations Na + , K + , Mg 2+ and Ca 2+ , and the heavy metal cations Fe 3+ , Al 3+ , Pb 2+ , Cu 2+ , Zn 2+ and Ni 2+ [12]; the major anions being Cl − and sulfates. The solid phase contains a plethora of insoluble organic compounds and different phosphate salts with the cations Mg 2+ , Ca 2+ , Al 3+ and Fe 3+ . The actual composition of the sewage sludge of interest mainly depends on two factors: the source of the wastewater, i.e. whether municipal or industrial wastewater is considered-which might underlie temporal changes-and the previous process steps and operating conditions in a wastewater treatment plant. The latter encompasses e.g. measures to reduce bulky sludge via addition of aid substances and the type of thickening devices (static or mechanical) including the use of flocculants.
The production of MAP from such feeds is essentially a recrystallization process that consists of two major process steps: (i) dissolution of the metal phosphate salts from the solid phase into the aqueous phase and (ii) precipitation of phosphate together with Mg 2+ and NH 4 + . This is achieved Fig. 1 Scheme of the Stuttgart process (C 6 H 8 O 7 : citric acid, HM: heavy metal), adopted from Ref. [25] 1 3 by applying several pH shifts. A scheme of the Stuttgart process, which is investigated in the present work, is presented in Fig. 1.
In the Stuttgart process, the phosphate salts are dissolved into the aqueous phase by lowering the pH value with sulfuric acid (H 2 SO 4 ). The remaining solid is removed by a filter, and only the aqueous phase is processed further. That aqueous solution usually contains less Mg 2+ than PO 4 3− , but a stoichiometric ratio is needed for maximizing the MAP precipitation. Hence, a Mg-containing salt, typically MgCl 2 or MgO, is added to the aqueous solution until the molar ratio Mg 2+ /PO 4 3− is slightly above 1. MAP is then precipitated by increasing the pH value with NaOH. In that process step, the precipitation of undesired salts other than MAP, e.g. iron phosphates, must be suppressed. The key difference between various processes for obtaining MAP from wastewater is the design of this suppression step. In the Stuttgart process [14], the aqueous metal cations are bound in complexes with citric acid (C 6 H 8 O 7 ) [15][16][17]. Citrate anions of different valency form strong ion pairs (electrically neutral and non-neutral ones) with different cations and thereby prevent these ions from forming salts. However, the complexation of one mole of metal cations of the sewage sludge roughly requires one mole of citric acid, which is a rather expensive chemical, at least compared to H 2 SO 4 and NaOH. Furthermore, not only the unwanted heavy metal cations, but also Mg 2+ is complexed by citric acid, which is obviously unfavorable for the MAP precipitation. A thorough understanding of the species distributions and the extent of the complexation of the different ions is thus essential for gaining an improved understanding of the process. Moreover, a potential reduction of the required amount of citric acid is highly desirable.
The recrystallization of phosphate salts is determined by the species distribution of aqueous phosphates, which depends crucially on the pH value of the solution [18]. However, the effort for measuring species distributions in the multi-component aqueous solutions that are relevant here is high. Therefore, methods for predicting the species distributions based on models are important, but only few such studies have been presented in the literature. Ohlinger et al. [9] and Wu and Zhou [19] predicted the MAP formation as a function of pH in systems with few species, considering only the ionic species that make up MAP (Mg 2+ , NH 4 + , and PO 4 3− ). Jia et al. [20] applied a chemical equilibrium model for the optimization of ammonium removal from digested sewage sludge as MAP. They considered the species Na + , K + , Mg 2+ , Ca 2+ , NH 4 + and PO 4 3− and calculated the amount of NaOH required for raising the pH value of the sludge to 9. However, their model does not include aqueous heavy metal cations and the solid heavy metal phosphates that might precipitate together with MAP, which have a negative influence on the fertilizer qualities of the product [21]. Song et al. [22] studied the pH-dependence of the crystallization of calcium phosphate, which needs to be considered in MAP production from digested sewage sludge.
In the present work, a thermodynamic model of the MAP production from digested sewage sludge is presented. In contrast to previous studies in the literature, in our model we consider most of the major ionic constituents of the sewage sludge. The model is tested using the Stuttgart process as an example. To validate the model predictions, they are compared to experimental data from the literature for the dissolution step of the Stuttgart process. Unfortunately, such a comparison cannot be carried out for the other process steps due to a lack of experimental data. The model is used to predict the species distributions and the precipitating solids during the precipitation step of the Stuttgart process. Finally, systematic variations of process parameters are carried out, revealing some optimization potential of the Stuttgart process with respect to the amounts of employed chemicals (H 2 SO 4 , NaOH, and citric acid).
In the model, chemical equilibrium constants and solubility products for the precipitating salts are taken from the literature. An extended Debye-Hückel model [23] is used for calculating the activity coefficients of the dissolved species. The considered species are Na + , K + , Mg 2+ , Ca 2+ , Pb 2+ , Cu 2+ , Zn 2+ , Ni 2+ , Fe 3+ , Al 3+ , NH 4 + , Cl − as well as the different forms of phosphate ions occurring at different pH values and their many mutual combinations to salts that might precipitate from the solution. In addition, the different Table 1 Species considered in the model developed in the present work and their ion-specific parameters σ i of the extended Debye-Hückel model taken from [23], see Eq. (3) No parameters are available for the non-neutral metal-citrate ion pairs, and hence those parameters were set to zero. The parameter σ i does not apply to neutral species and hence a dash (-) is given for these

Assumptions and Considered Species
In the present work, the aqueous phase of the sewage sludge is assumed to contain the species listed in Table 1. Insoluble organic compounds are disregarded. In the aqueous solution, chemical reactions take place, which are listed in Table 2 along with their equilibrium constants. The aqueous solution can be in thermodynamic equilibrium with one or more solid phases. It is assumed that the several solids that might precipitate from the solution form individual phases of the individual pure salts. The considered salts are listed in Table 3.

Thermodynamic Model
The thermodynamic model developed in the present work is conceptually similar to a model employed in a previous work for describing the crystallization sequence during the evaporation of seawater [24]. We start the model description by briefly introducing the employed normalizations of the chemical potentials. Since we only consider processes at 298.15 K and 101.3 kPa, in all following equations we drop the temperature and pressure dependence of the chemical potentials to simplify the notation.
The chemical potentials of all solutes i are normalized to infinite dilution according to Henry's law on the molality scale as shown in Eq. (1): where 0 i is the reference chemical potential of solute i in the hypothetical ideal solution at a molality of 1 mol kg −1 , b i is the molality of solute i, the unit molality b 0 = 1 mol kg −1 makes the second term dimensionless, and γ i is the activity coefficient of solute i on the molality scale. The chemical potential of the solvent water is normalized according to Raoult as shown in Eq. (2): where 0 W is the chemical potential of pure water and a W is the activity of water.
In the present work, the activity coefficients of all ionic species are described with the extended Debye-Hückel model according to Kielland [23] as shown in Eq. (3), in which all variables are dimensionless and which only applies to 298.15 K:  Here, A = 1.1658 is the Debye-Hückel constant, z i denotes the charge of ion i, and I is the ionic strength of the solution defined in Eq. (4): In Eq. (3), σ i is the so-called "distance of closest approach" parameter, which is ion-specific. Values for σ i are taken from the literature and are compiled in Table 1. For some non-neutral ion pairs, no σ i are available. In order to keep the model as simple as possible, σ i was set to zero in these cases. Note that when setting σ i to zero, the Debye-Hückel limiting law for the species i is recovered from Eq. (3). Furthermore, for the sake of simplicity, neutral ion pairs and neutral species are considered to behave ideally, i.e. γ i = 1 for these components. These simplifying assumptions, and the use of an extended Debye-Hückel approach in the first place, are mainly motivated by the fact that the maximum ionic strength of the solutions studied here is about 0.2 mol kg −1 . For such solutions, the present approach should be a reasonable approximation to the actual behavior of the components in solution. The water activity is modeled as shown in Eq. (5) where M W is the molar mass of water and the sum runs over all solute species.
From the described activity coefficient model, the thermodynamic definition of pH shown in Eq. (6) can be evaluated straightforwardly.
The chemical equilibrium of a reaction R is characterized by its equilibrium constant K R defined in Eq. (7): Here, m are the reactants and n are the products of the reaction, b i and γ i are the molality and activity coefficient on the molality scale of component i in the solution, and ν i is the stoichiometric coefficient of component i in the reaction. The reactions considered in the present work and their equilibrium constants are listed in Table 2.
The only exception to Eq. (7) is the auto-protolysis of water, for which the equilibrium constant is defined in Eq. (8) because the chemical potential of water is normalized according to Raoult (cf. Eq. (2)).
In the different steps of a phosphorus recovery process, different salts precipitate from the solution or dissolve into the solution. The following inequality Eq. (9) is used to determine whether a solid salt S is present or not: Here, K SP,S is the solubility product and K AP,S is the so-called activity product. These properties are defined in Eqs. (10) and (11) as [38] and Here, n are the ionic species that form salt S, ν n is the stoichiometric coefficient of ion n in the salt S, and 0 n is that ion's standard chemical potential. Furthermore, b n and γ n are the molality and activity coefficient of ion n in the aqueous solution, a W is the water activity in the solution, and υ W,S is the number of moles of crystal water in one mole of S.
The activity product K AP,S of a salt S is a formal property. It is a function of the composition of the aqueous solution. By contrast, the solubility product K SP,S is a constant. Equation (9) means the following: if equality between the solubility product and the activity product holds, the aqueous solution is in equilibrium with the salt S. If the greater-than sign  holds, this is not the case, i.e. the solution is not saturated with the salt S. In the present work, metastability is not taken into account, i.e. if the solubility product is reached in the aqueous solution, the salt immediately precipitates. The salts considered in the present work and their solubility products, which were taken from the literature, are listed in Table 3.

Simulations of the Stuttgart Process
The developed model was used for predicting the species distributions as a function of pH in the different steps of the Stuttgart process. The initial composition of the sewage sludge used in the simulations is presented in Table 4. For the aqueous phase, the initial composition was obtained from the average values of digested sludge reported in Ref. [12] as this reference provides values for a wide range of components. To that composition, Cl − was added until the electroneutrality condition was met. For the solid phase, no detailed analysis is available in the literature. Hence, for carrying out the simulations it was assumed to comprise the PO 4 3− salts of Mg 2+ , Ca 2+ , Al 3+ , and Fe 3+ for the following two reasons: first, these metal ions are found to dissolve into the solution in the experiments in Ref. [25]. Second, their PO 4 3− salts are practically insoluble in water, so that they are likely to be present in the solid phase (cf. the solubility products reported in Table 3). The amounts of these salts in the solid phase of the sewage sludge were chosen such that in the simulations of the dissolution step, the maximum metal dissolution (at low pH values) agrees with the experimental data from a sewage plant in Offenburg, Germany [25].
Starting with this initial composition of the sewage sludge, simulations were carried out for investigating the Stuttgart process, cf. Figure 1 for the process scheme. Following the experiments reported in Ref. [25], the consecutive process steps were simulated as follows: for the dissolution step, H 2 SO 4 was added to the sewage sludge, and the pH value and the species distribution were calculated. The addition of H 2 SO 4 was stopped at a pH value of 4. Then, mimicking the filtering process, the solid phase was simply discarded from the simulation. In the next step, MgCl 2 was added to the solution until the molar ratio Mg 2+ /PO 4 3− was equal to 1. Then, the complexation was carried out by adding citric acid. The amount of citric acid was stoichiometric with respect to the total amount of dissolved metal cations in solution (Mg 2+ , Ca 2+ , Al 3+ and Fe 3+ ), see Eq. (12): As a consequence of the addition of citric acid, the pH value drops to about two. Finally, for the precipitation step, NaOH was added to the solution, and the pH value and the species distribution were calculated. NaOH is assumed to dissociate completely and was added until a pH value of 8.5 was reached. During this precipitation step, all the MAP that can potentially precipitate from solution does so.

Optimization Studies
Carrying out parameter variations in experiments is timeconsuming and expensive. The model from the present work provides the opportunity for systematically studying the influence of the key parameters of the Stuttgart process by simulations and to thereby gain guidance on how the process might be optimized. In the present work, two types of variations were studied. These variations pertain the influence of the amounts of the employed chemicals H 2 SO 4 , NaOH and citric acid.
First, the pH values at which the dissolution step and the precipitation step are stopped were systematically varied in the ranges 2-5 and 6-9, respectively. All other process parameters remained unaltered. A higher pH value at the end of the dissolution step means that less H 2 SO 4 is needed; a lower pH value at the end of the precipitation step means that less NaOH is needed. As the objective for optimization, the purity of the produced MAP was considered, since both chemicals are inexpensive and obtaining a higher quality product clearly outweighs the added cost if more H 2 SO 4 or NaOH are needed for that purpose. The purity of MAP is defined as the mass fraction of MAP among all precipitated solids.
Second, the amount of added citric acid was varied. A potential reduction of the citric acid demand of the process is highly desirable since citric acid is the most expensive chemical used in the process. Hence, simulations were carried out in which the ratio k between the molality of citric acid and the total molality of the dissolved metal cations in solution (Mg 2+ , Ca 2+ , Al 3+ and Fe 3+ ), defined in Eq. (13) as was varied. In the Stuttgart process as discussed in the previous section, k = 1. Again, the purity of MAP was computed, but also the total cost for the chemicals per kg of produced MAP was estimated. The cost of the individual chemicals used in these calculations are compiled in Table 5.

Implementation
The thermodynamic model described above was implemented in MATLAB. The process simulations and optimization studies were conducted with MATLAB as well, solving the sets of non-linear equations (mass balances and dissociation/reaction equilibria) with the MATLABsolver fsolve.

Simulations of the Stuttgart Process
In Fig. 2, the model predictions for the dissolution step of the Stuttgart process are compared to the experimental data from Ref. [25], where the molalities of several metal cations in the aqueous phase were measured as a function of the pH value. The experimental data from Ref. [25] are compiled in Table 6 for the reader's convenience.
For the initial composition of the sewage sludge, the present model yields a pH value of about 7.4, which compares favorably to the experimental data from Ref. [25], where the pH value of the digested sewage sludge (input to the Stuttgart process) is reported to be in the range 7.4-7.5. Considering that the initial composition of the sludge had to be taken from a different literature source [12], this favorable agreement is gratifying and indicates that both sources are commensurate. Moreover, the accuracy of prediction is similar to that presented by Zhou and Wu [19]. Figure 2 then shows that with decreasing pH as a result of adding H 2 SO 4 , the molalities of the ions in the aqueous phase increase due to the dissolution of phosphate salts from the solid phase of the sewage sludge. This can be explained by the combination of two effects: first, when lowering the pH value, the phosphorus species distribution in the liquid phase changes. Second, the solubility products of the metal HPO 4 2− salts exceed those of their metal PO 4 3− counterparts by many orders of magnitude, cf. Table 3. The combination of these two effects leads to the dissolution of the salts into the aqueous phase and, consequently, also to an increase of the metal ion concentrations. To take a specific example: the solid Ca 3 (PO 4 ) 2 dissolves as CaHPO 4 a pH value of about 5.5, which is visible as a "kink" in the molality of Ca 2+ . Similar kinks in the molalities of the other ions appear in Fig. 2 whenever the dissolution of a new salt begins. At a pH value of about 3, the phosphate salts with the cations Mg 2+ , Ca 2+ and Fe 3+ are completely dissolved in the aqueous phase and hence, the molalities of these ions remain constant.
Furthermore, Fig. 2 shows good agreement between the experimental data and the model predictions. Especially the molalities of Fe 3+ and Ca 2+ are predicted accurately by the model. At high pH values, the model overestimates the molalities of Mg 2+ , while the molalities of Al 3+ are underestimated for all pH values. Nonetheless, considering the simplicity of the model and the fact that no adjustable parameters were used and all calculations are entirely predictive, the overall agreement is favorable. Moreover, it must be borne in mind that the initial sludge composition and the ion molalities stem from two different literature sources,  Table 6. Error bars indicate the experimental uncertainties and are only shown where they exceed the symbol size Table 6 Experimental data for the ionic molalities of Fe 3+ , Al 3+ , Mg 2+ and Ca 2+ in the dissolution step of the Stuttgart process as a function of pH. The numbers shown here are averages from the series of three measurements reported in Ref. [24]  which probably also results in some discrepancies between the experimental data and the simulation results. The next step in the Stuttgart process is the combined complexation of the heavy metal ions and the precipitation of phosphorus as MAP. We start by discussing the heavy metal complexation. Figure 3 shows the molalities of the citrate complexes of Ca 2+ and Al 3+ as a function of pH. Figure 4 shows a similar plot for the citrate complexes of Mg 2+ and Fe 3+ .
In the following, we only discuss Fig. 3 (Ca 2+ and Al 3+ ) because Fig. 4  The key aspect, however, is that at high pH values, all the Al 3+ is bound in the complex with C 6 H 4 O 7 4− , but most of the Ca 2+ still remains in the free form. This distinctly different behavior is likely caused by the higher valency of Al 3+ compared to Ca 2+ : Al 3+ attracts the various citrate anions much more strongly than Ca 2+ does. As a consequence, the equilibrium constants of the formation reactions of the Al 3+ -citrate complexes are about four orders of magnitude higher than those of the Ca 2+ -complexes. Hence, at high pH values the complexation of the Al 3+ ions is pronounced, while only some of the Ca 2+ ions form citrate complexes.
A similar line of arguments holds for the comparison between Fe 3+ and Mg 2+ shown in Fig. 4. Note that the increase in the molality of Mg 2+ when raising the pH value from 2 to about 3 in Fig. 4 is caused by the addition of MgCl 2 in the same process step.
Also the other heavy metal ions (not shown in the Figures) form strong complexes with the different anions of citric acid. According to the model predictions, at a pH value of 9, there are practically no remaining free Fe 3+ and Al 3+ ions, more than 90% of the Cu 2+ , Zn 2+ and Ni 2+ ions and roughly 60% of the Pb 2+ ions are complexed. These findings compare favorably to experimental results from the literature [26].   Fig. 4 for the different Mg 2+ species). These kinks are caused by the precipitation of different Ca 2+ -containing salts, because complexation and precipitation occur simultaneously in the Stuttgart process. To discuss the precipitation in more detail, all precipitating salts are shown as a function of pH in Fig. 5.
With increasing pH value, first the solubility limits of the hydrogen phosphates of Mg 2+ and Ca 2+ are reached, and these salts precipitate from solution. However, this changes abruptly at a pH value of 6 and above, where PO 4 3− becomes the dominating phosphate species in solution. Since the salts containing PO 4 3− are only sparingly soluble in aqueous solution, they start to precipitate immediately. Most notably, the desired product MAP is the main precipitating species. Ca 3 (PO 4 ) 2 co-precipitates with MAP, but the amounts of Ca 3 (PO 4 ) 2 are almost one order of magnitude lower than those of MAP. No salts of heavy metal ions precipitate since they are effectively complexed by citric acid. Hence, the process allows for obtaining MAP with high purity, the only notable contaminant being Ca 3 (PO 4 ) 2 , which has no negative implications on the fertilizing abilities of the product. At a pH value of 8 and above, also Mg 3 (PO 4 ) 2 starts to precipitate, which comes along with a reduction of the MAP precipitation. Hence, too high pH values should be avoided. Figure 6 shows the results obtained with the present model for the purity of MAP as a function of two pH values: the pH value at the end of the dissolution step and the pH value at the end of the precipitation step.

Optimization Studies
With the operating parameters as in the Stuttgart process from Ref. [25], the maximum possible amount of MAP precipitates from solution (not shown in the Figure), and a MAP purity of about 0.72 g g −1 is achieved. The simulations suggest, however, that an increase up to about 0.80 g g −1 is possible by lowering both the pH value of the dissolution step and that of the precipitation step (but to a lesser extent). As result, more H 2 SO 4 is needed, but approximately the same amount of NaOH is required for achieving these new pH values. Since H 2 SO 4 is rather cheap, this seems to be a viable option. Figure 7 shows the influence of the stoichiometric ratio k between citric acid and the dissolved metal ions, cf. Equation (13), on the price and the purity of the produced MAP.
Both the cost for the chemicals and the purity of the produced MAP increase with increasing k. The cost increases almost linearly with k; the slightly non-linear behavior at low k is caused by the fact that in this region, the amount of produced MAP increases with increasing k. For k > 0.25, all the MAP that can potentially precipitate does so. In contrast to the cost, the MAP purity increases quickly with increasing k and remains constant from k = 0.5 onwards. This plateau of the MAP purity is intriguing: it suggests that almost the same MAP purity as in the Stuttgart process (k = 1) can be achieved with k = 0.5, i.e. with a reduction of the needed amount of citric acid by 50%. As a result, also the cost per kg of produced MAP drops by almost 50%. Hence, a reduction of the use of citric acid could significantly improve the economics of the process.

Conclusions
A thermodynamic model for sewage sludge systems was developed and used for investigating the different steps of a phosphorus recovery process for MAP production, namely the Stuttgart process. The model is based on an extended Debye-Hückel approach and literature values for equilibrium constants and solubility products. All results obtained with the model are strict predictions, which were validated by comparison to available experimental data. The model can be used to gain a physico-chemical insight into the Stuttgart process and can thereby help with finding the correct ranges for the dosage of the employed chemicals. For example, the model revealed that MAP precipitation is maximized when NaOH is added in the precipitation step up to a pH value between 6.5 and 8; at lower pH, no MAP precipitates and at higher pH, other phosphate salts start to co-precipitate with MAP. Furthermore, the model can be used for process optimization. The results of calculations with systematic variations of process parameters suggest that there is room for improvement of the process with respect to product purity and cost. Most notably, the consumption of citric acid might be reduced to about half while keeping the purity of the produced MAP almost constant. Altogether, the results of the present work show that thermodynamic modeling is viable tool in the analysis and design of phosphorus recovery processes.