Development, implementation, and validation of a generic nutrient recovery model (NRM) library

Abstract The reported research developed a generic nutrient recovery model (NRM) library based on detailed chemical solution speciation and reaction kinetics, with focus on fertilizer quality and quantity as model outputs. Dynamic physicochemical three-phase process models for precipitation/crystallization, stripping and acidic air scrubbing as key unit processes were developed. In addition, a compatible biological-physicochemical anaerobic digester model was built. The latter includes sulfurgenesis, biological N/P/K/S release/uptake, interactions with organics, among other relevant processes, such as precipitation, ion pairing and liquid-gas transfer. Using a systematic database reduction procedure, a 3- to 5-fold improvement of model simulation speeds was obtained as compared to using full standard thermodynamic databases. Missing components and reactions in existing standard databases were discovered. Hence, a generic nutrient recovery database was created for future applications. The models were verified and validated against a range of experimental results. Their functionality in terms of increased process understanding and optimization was demonstrated.


Introduction
In the transition from waste (water) treatment plants (WWTPs) to water resource recovery facilities (WRRFs), mathematical models are becoming important tools to hasten nutrient recovery process implementation and optimization . Indeed, models may aid in technology development, process operation, optimization, and scale-up in a cost-effective way (Rieger et al., 2012;Yu et al., 2011). Although to date many processes for the recovery of nutrients from waste(water) have been proposed and applied to varying degrees, for example, struvite precipitation, ammonia stripping and acidic air scrubbing (Vaneeckhaute et al., 2017a), no generic models for nutrient recovery aiming at the construction and optimization of treatment trains for resource recovery are currently available (Vaneeckhaute et al., 2017b). Moreover, existing model libraries for WWTPs, e.g., activated sludge models (ASMs) provided by the International Water Association (IWA) (Henze et al., 2000;Gernaey et al., 2004), do not allow the integration of nutrient recovery unit processes and/or the coupling of a nutrient recovery treatment train. This is due to the omission of key fundamental physicochemical components and transformations that are essential to describe nutrient recovery unit processes (Batstone et al., 2012;Brouckaert et al., 2010;Flores-Alsina et al., 2015). Critical elements to be dealt with include accurate descriptions of acid-base reactions, slow precipitation kinetics, liquid-gas exchange and sorption/desorption in the complex mixture of chemical species that the resource recovery systems in place deal with (Batstone et al., 2012). Consequently, the potential to use models to adequately put together an optimal treatment train of unit processes and set the operating conditions that maximize nutrient recovery and fertilizer quality is missing.
Over the last 10 years, important progress was made towards the development and integration of a physicochemical modelling framework compatible with the current more biological processoriented modelling frameworks provided by IWA Hauduc et al., 2015;Lizarralde et al., 2015;Mbamba et al., 2016). However, the scope of the existing studies stops at the anaerobic digestion of WWTP sludge, where it mainly aims at the prediction of uncontrolled struvite precipitation during digestion through phosphorus modelling. At the start of this research, no generic models were available that may allow to predict, optimize and control under dynamic conditions the recovered product quality (e.g., macronutrient content, particle size, density), yield and process performance of a series of nutrient recovery technologies following digestion of various waste(water) flows (manure, sludge, food waste, etc.).
The reported research aimed at developing a library of generic integrated biological-physicochemical three-phase mathematical process models for the most established nutrient recovery systems currently available as selected in Vaneeckhaute et al. (2017a), as well as a compatible model for anaerobic digestion. The models are based on detailed solution speciation and reaction kinetics, as brought forward in Vaneeckhaute et al. (2017b). This nutrient recovery model (NRM) library is a synthesis of the large body of knowledge on nutrient recovery processes that is currently available from research studies and operational experience. In contrast to existing model libraries for waste(water) treatment, e.g., the ASM library (Henze et al., 2000;Gernaey et al., 2004), the scope of the proposed NRM library starts at the anaerobic digester and focusses on the nutrient recovery treatment train following the digester.
In addition to the development of a generic physicochemical modelling framework, a critical and challenging step when combining biological and physicochemical differential equations is their numerical solution (Lizarralde et al., 2015). This is due to the stiffness that arises when considering reactions with very different conversion rates, i.e. the range of system time constants is large (Batstone et al., 2012;Brouckaert et al., 2010;Lizarralde et al., 2015;Musvoto et al., 2000;Rosen and Jeppsson, 2006;Sotemann et al., 2005). Previous attempts towards inclusion of a physicochemical modelling framework in existing WWTP models (e.g., Flores-Alsina et al., 2016;Hauduc et al., 2015;Mbamba et al., 2016;Tak acs et al., 2006) applied a limited literature-based selection of chemical species and reactions for self-implementation in the modelling software and used self-coded numerical solvers that have shown difficulties with convergence . Moreover, when one wants to extend these models with new species and reactions, time-consuming and complicated programming work is required. Model flexibility is, however, particularly important for modelling of WRRFs due to the variability of waste (water) flows in time and between different facilities . Lizarralde et al. (2015) proposed the coupling of an external existing geochemical software tool for inclusion of some basic speciation calculations in dynamic process models. The use of an external geochemical software tool with designated thermodynamic databases is interesting for accurate calculation of chemical speciation and pH. Software tools as PHREEQC and MINTEQ are generally accepted tools for equilibrium water quality modelling and have a dedicated and proven solver for chemical speciation calculations (Allison et al., 1991;Parkhurst and Appelo, 2013). However, simulation times using the full PHREEQC/ MINTEQ thermodynamic databases for chemical speciation may be longer than when dedicated program code is used (Lizarralde et al., 2015). Hence, an important challenge exists in the development of an efficient methodology for including and solving the stiff equations related to the chemical speciation submodel in nutrient recovery models. A compromise should be found between model accuracy and simulation times.
The present paper describes the specifications, the development methodology and the implementation of the generic NRM framework. A systematic procedure to allow for inclusion of accurate chemical speciation in dynamic nutrient recovery process models at minimal computational effort is proposed. Model functionality in terms of increased process understanding and optimization is demonstrated through testing and validation. Recommendations for further experimental research required to fully calibrate the model dynamics, as well as case-specific potential model extensions, are provided.

Nutrient recovery model (NRM) building methodology
The methodology used for development of the NRM library can be represented by six steps, shown in Fig. 1 and described in detail in the sections below.
The proposed generic models are based on mass balances to describe physicochemical and biochemical transformation and transport processes, as well as an accurate calculation of water chemistry in order to correctly define solution speciation and driving forces for component transformation. Two key features of the models should be stressed. First, a dynamic modelling approach, i.e. one that accounts for time-dependent changes in the state of the system, was applied, because the models should be applicable to time-varying situations and variable operating conditions, such as i) periodical load variations, e.g., truck loads of waste, sludge treatment during working hours only, and seasonal variations, ii) individual disturbances, e.g., rain events and incorrect manipulations, and iii) systems that are operated intermittently or cyclically as is the case for multiple nutrient recovery processes, e.g., intermittent aeration in stripping systems and (semi-) batch processes to obtain target fertilizer specifications, e.g., a predefined ammonium sulfate (AmS) concentration in an acidic air scrubber. Second, the geochemical software tool PHREEQC was used for two purposes in the development of the NRM library: 1) PHREEQC for NRM building (Section 2.2.1), which involves the selection of species and reactions to be included in the models, the preparation of a reduced PHREEQC model database, and the definition of PHREEQC selected outputs; 2) PHREEQC for NRM simulation (Section 2.3), which involves the tight coupling of the reduced PHREEQC model to a kinetic and dynamic mass balance model in order to accurately and efficiently calculate speciation and driving forces for component transformations at each time step during the model simulations.
As opposed to previously used speciation modelling methodologies (e.g., Flores-Alsina et al., 2016;Hauduc et al., 2015;Mbamba et al., 2016;Tak acs et al., 2006), the proposed methodology guarantees convergence and flexibility (Lizarralde et al., 2015). In order to reduce simulation times, a systematic procedure for thermodynamic model database reduction was proposed. Finally, it should be noted that in the following sections, variables will be defined with their dimension given in straight brackets: ½M for mass, ½L for length, and ½T for time.

2.1.
Step I: definition of modelling objectives 2.1.1. Selection of considered/included unit processes and input waste streams A literature review on nutrient recovery technologies (Vaneeckhaute et al., 2017a) was conducted in order to select the best available technologies as key unit processes for modelling (Table 1: four key units). The selection was made based on the economic feasibility, full-scale application at this stage, and the potential to produce marketable end products for agricultural applications. With the purpose of modelling treatment trains, four ancillary units were additionally selected ( Table 1).
As mentioned above, in contrast to existing studies, the scope of the present research starts at the anaerobic digestion unit and focusses on the nutrient recovery treatment train following the digester. No recycle flows to upstream facilities in the WRRF, e.g., to an activated sludge (AS) system, were currently considered. In later stages, the proposed NRM models could be coupled to activated sludge models (ASMs), if a generic physicochemical framework is also integrated in the ASMs.
As input waste stream to the digester, manure and sludge (primary and secondary sludge, and mixtures of these) from WWTPs removing nitrogen (N) and chemical oxygen demand (COD) were considered. Digestate, the remaining product after digestion, was considered as input stream to the key units for controlled nutrient recovery following the digester. Next to manure, WWTP sludge was selected since the current most advanced models for anaerobic digestion originate from the municipal wastewater and sludge treatment sector. Nevertheless, for future applications, the generic NRM-AD implementation allows easy extension to co-digestion of other organic-biological wastes, e.g., using the general integrated solid waste co-digestion (GISCOD) modelling tool proposed by Zaher et al. (2009). The NRM-AD model can also be extended to allow for specific reactions occurring during the treatment of sludge from enhanced biological phosphorus (P) removal (EBPR) as, e.g., in Ikumi (2011) or Wang et al. (2016), but this was considered to be outside the scope of this paper.

Specification of model outputs and influencing factors
In order to develop valuable tools for process optimization, the desired model outputs and factors that may affect these outputs were defined for each NRM key unit on the basis of a detailed literature review (Vaneeckhaute et al., 2017b).
Obviously, the total content of principal macronutrients, N, P, and potassium (K), in the fertilizer product and the amount of biogas produced are important model outputs, so as to quantitatively and qualitatively determine the overall resource recovery. Next to the three principal macronutrients, N, P, and K, previous studies have shown the relevance of the secondary macronutrient, sulfur (S), in the context of nutrient recovery . Some motivations for inclusion of S in the models were: i) the demand for S fertilization in agriculture is increasing, hence its recovery deserves Step I: definition of modelling objectives; Step II: theoretical model development; Step III: model implementation and numerical solution; Step IV: data collection and identification of data needs; Step V: model application and validation; Step VI: scenario analyses and process optimization; MSL ¼ model specification language.
attention (Till, 2010), ii) S may precipitate with iron (Fe), making Fe less available for P precipitation, iii) sulfate reducing bacteria (SRBs) compete with methane (CH 4 ) producing bacteria for the same substrate, hence CH 4 production may be reduced at high S concentrations (Oyekola et al., 2007), iv) hydrogen sulfide (H 2 S) is an important inhibitor of CH 4 producing bacteria (Oyekola et al., 2007), and v) high biogas H 2 S values cause important concerns (toxicity, corrosion, biogas pollution), e.g., in the paper industry (Reiter and Piccot, 2004). Calcium (Ca) and magnesium (Mg) are also of importance, mainly for their soil improving properties and their interaction with P (Vaneeckhaute et al., 2016).
For all nutrient recovery systems, the percentage recovery of the target nutrient is a key performance measure. It was calculated using Equation (1): in which S in i and S out i are the in-and outgoing liquid flow concentrations for component i ½M L À3 , and Q in and Q out are the in-and outgoing flow rates ½L 3 T À1 .
Furthermore, the macronutrient use efficiency (N, P, K, S) in the fertilizer end products is an important factor in determining the agronomic potential and sustainability of the produced fertilizers. It was evaluated as the percentage available or mineralized nutrient content over the total nutrient content, e.g., NH 4 -N/total N and ortho-P/total P. This percentage can be obtained by means of a chemical solution speciation calculation (Section 2.2.1). Next, the fertilizer pH and salt content are of important concern as they may impact soil quality. The pH was directly calculated from solution speciation. Salts were characterized using the sodium adsorption ratio (SAR), i.e. the relative amount of available sodium (Na) over divalent cations, Ca and Mg (Hillel, 2008).
Factors that may additionally determine the value of the recovered product are the particle size (for solid fertilizers), the density (for liquid fertilizers), and the product purity. In this work the particle size was evaluated as mean particle diameter (Section 2.2.2), but in future research one may be interested in particle size distributions (PSDs) (Nopens et al., 2014;Perez et al., 2008).
For the NRM-Prec unit, product purity was evaluated by calculating the fraction of precipitated target mineral(s) over the total product collected, taking in account the presence of multiple competing and concurrent precipitation reactions. To this end, also the precipitation of principal micronutrients occurring in waste(water) treatment, such as Fe and aluminium (Al), were evaluated, since these precipitates may negatively impact the fertilizer P release in the soil (Vaneeckhaute et al., 2016). Moreover, pollution with organics was accounted for (see Section 2.2.1). For the NRM-Strip/NRM-Scrub units, purity was evaluated by calculating the amount of volatile target component(s) captured over the total amount of gas/liquid captured.
Finally, the formation of scale within the treatment unit is an important operational bottleneck for multiple nutrient recovery technologies. Especially calcium carbonate (CaCO 3 ) and magnesium carbonate (MgCO 3 ) formation in the stripping and scrubbing units are of concern. To determine scale formation, the amount of CaCO 3 and MgCO 3 precipitates formed were evaluated, next to other relevant precipitation reactions. The scaling potential was then examined by using the scaling criteria of the Ryzner Index (Tchobanoglous et al., 2003).

2.2.
Step II: theoretical model development The dynamic mathematical model of each unit process was built using (Fig. 2): i) the definition of a chemical speciation model by means of geochemical modelling software (PHREEQC for model building, Section 2.2.1), ii) the description of a kinetic physicochemical and biochemical transformation model tailored to the models developed in the first step (Section 2.2.2), and iii) the selection of a reactor mass balance model to describe the (timedependent) process conditions (Section 2.2.3).

Chemical speciation model: PHREEQC for NRM building
In order to describe the water chemistry in each system, first the potentially present chemical components and species were identified (step 1), and the possible heterogeneous physicochemical Fig. 2. Development of combined physicochemical-biological three-phase (liquid-solid-gas) process models. COD ¼ chemical oxygen demand; G ¼ gas; P ¼ precipitate; Q_gas ¼ gas flow rate; Q_liq ¼ liquid flow rate; Q_prec ¼ precipitate extraction rate (for NRM-Prec); S ¼ soluble; X ¼ biological particulate COD. transformation reactions (gas transfer, precipitation) were selected using generally accepted geochemical software for equilibrium water quality modelling, PHREEQC 3.0.6 (Parkhurst and Appelo, 2013). Visual MINTEQ 3.1 was used as a control (Allison et al., 1991). Since the involved homogeneous reactions (acid-base, ion pairing) in a speciation calculation are very rapid compared to heterogeneous physicochemical reactions and biological reactions, instantaneous equilibrium can be assumed adequate for solving water chemistry in NRMs (Batstone et al., 2012).
In order to compromise between model accuracy and simulation times when coupling the speciation model to the dynamic mass balance model, a reduced PHREEQC database and input script with defined selected model outputs were developed for each key unit (Section 2.3.1). The four-step procedure proposed for NRM building, involving the selection of the relevant species/reactions and the preparation of the reduced PHREEQC chemical speciation model, is presented in Fig. 3 and further described below. 2.2.1.1.
Step 1 -selection of relevant components for each unit process. Based on literature, collected experimental data and prior knowledge, the most important physicochemical dissolved components to include in models for nutrient recovery from both (digested) manure and sludge were selected for each key unit process (Table 2). In line with the selected model outputs (Section 2.1.2), it was aimed to represent five important component classes: 1) All important macronutrients for recovery in line with the findings in Vaneeckhaute et al. (2014) (~impact recovery efficiency and fertilizer value); 2) Gaseous compounds (~impact biogas production, volatilization, odors, greenhouse gas emissions, among other); 3) Salts (~impact ionic strength and soil quality); 4) Inorganic and organic carbon compounds (~impact biogas production, product purity, and scaling); 5) Micronutrients that may occur in large quantities in waste(water) treatment, e.g., Fe and Al as a result of coagulation/flocculation practices (~impact product purity and recovery potential).
Since redox reactions were also considered, components that exist in more than one valence state in solution were identified by their component name followed by their valence. For instance, i) the component S_C_4_ (carbon þIV) constitutes CO 3 2À plus HCO 3 À plus H 2 CO 3 (or CO 2,aq ) plus various other carbonate complexes present in the solution, such as MgCO 3 and CaHCO 3 þ , and ii) S_N_min3_ (nitrogen eIII) constitutes both NH 4 þ and dissolved NH 3 , as well as its various complexes. Only for Fe, the two valence states, Fe (þII) and Fe (þIII), were lumped together into one component for total Fe, since the measurement of its valence is complicated and generally not provided in practice in WRRFs, nor in literature. Yet, in the speciation calculation, the Fe (þII) /Fe (þIII) redox equilibrium was considered, as calculated from the occurring redox potential. The input Fe redox states, e.g., Fe(þIII)Cl 3 and Fe(þII)SO 4 , can optionally be specified, if such data are available. As it is well-known that the presence of organic compounds may influence the purity of recovered products (Kozic et al., 2011), relevant interactions between inorganic and organic components were also accounted for. Among the organic biological components considered (see Section 2.2.3), volatile fatty acids (VFAs) up to valerate were included as individual components in the physicochemical models. Oh and Martin (2010) indeed emphasized the particular importance of their physicochemical behaviour in WRRFs. The remaining soluble organic chemical oxygen demand (COD) fractions (see Section 2.2.3) were lumped into one component, i.e. dissolved organic carbon (DOC; 1 g DOC z 0.33 Â g COD). For DOC, the complexation with metals (Ca, Mg) was computed using a competitive Gaussian model for dissolved organic matter (DOM; 1 mol DOC z 8.6 Â 10 À2 mol DOM; USEPA, 1999). This simplified approach may be further refined for future applications, if more insights in the physicochemical behaviour of each particular COD fraction become available.
Finally, it should be remarked that heavy metals, such as cadmium (Cd), copper (Cu), and zinc (Zn), were not yet included in the speciation models. Nevertheless, heavy metals and the corresponding reactions are available in PHREEQC. Hence, the generic approach used for chemical speciation allows easy extension of the models to incorporate heavy metals for future applications.

2.2.1.2.
Step 2 -addition of relevant components/species/reactions to generic geochemical databases. To verify completeness, the generic PHREEQC (Phreeqc.dat) and MINTEQ (minteq.v4.dat) databases were compared with each other, as well as with prior knowledge and with literature. Two observations were made: 1) the generic MINTEQ database is more complete than the PHREEQC one in view of WRRF modelling, 2) some important components, species, and reactions that can be expected in WRRFs are not included in either database. Hence, the generic database files were extended prior to use for speciation calculation (Table 3). The corresponding acidbase constants, ion pairing constants, solubility products, and other thermodynamics were taken from literature or other model libraries, as indicated in Table 3.
It should be noted that in the context of nutrient recovery from waste(water) flows as fertilizer products, the database extensions provided concern a fundamental contribution to the field. For example, K-struvite is, next to N-struvite, an interesting fertilizer, though its precipitation reaction is not included in the standard databases. Also precipitation of aluminium phosphate (AlPO 4 ) is highly important in waste(water) treatment since Al-salts are often dosed for sludge conditioning, whereas the precipitation reaction of ammonium sulfate, (NH 4 ) 2 SO 4 , is essential for description of the scrubbing process. Noteworthy is also the clear impact of the omission of the species monosodium phosphate, i.e. NaH 2 PO 4 (aq), on the simulation results, as was observed during model validation of the NRM-Prec (see Section 3.4.1). The generic extended database in view of nutrient recovery was named 'Nutricover.dat' and will be made available for inclusion in future PHREEQC and MINTEQ software packages.

2.2.1.3.
Step 3 -setting up the speciation submodel ¼ selection of relevant species and reactions. The following methodology was used for selection of the relevant species and reactions: A. Specification of input scenarios (components þ operational conditions); B. Run PHREEQC under the various conditions defined in A; C. Select relevant species and reactions based on the PHREEQC outputs; D. Verify the selection of species and reactions with literature.
A. Specification of input scenarios: Realistic ranges for the input component concentrations and operational conditions (e.g., pH and temperature) for the speciation calculations were adopted from literature and experimental data as described in Section 2.4, as well as through contact with technology providers. The operational conditions and input streams tested for each key unit process are the following: ➢Anaerobic digestion: no oxygen, pH: 5e8.5, temperature: 20e55 C, input: sludge and manure; ➢Precipitation unit: pH: 7e11, temperature: 20e50 C, with and without Ca(OH) 2 , CaO, MgCl 2 , Mg(OH) 2 , or MgO dosing (0e500 mol m À3 ), input: digestate; ➢Stripping unit: pH: 7e11, temperature: 20e70 C, with and without NaOH, Ca(OH) 2 , CaO, Mg(OH) 2 , or MgO dosing for pHincrease (0e500 mol m À3 ), input: digestate; ➢Air scrubber: H 2 SO 4 -solution at pH: 1e4 and temperature: 15e25 C, input: stripped air. Values between brackets represent the use of air instead of chemicals for pH-adjustment.
PHREEQC makes calculations using an input script in which the problem is specified via 'KEYWORDS' and associated data blocks. First, all possible realistic scenarios were introduced using the maximum/minimum values of all considered operational factors and input variables for each unit separately. Next, for each unit the composition of 20 different possible input flows (from literature: Astals et al., 2013;Bhuiyan et al., 2007;Cesur and Albertson, 2005;Martin, 2003;Mattocks et al., 2002;Tchobanoglous et al., 2003;Vaneeckhaute et al., 2012Vaneeckhaute et al., , 2013aVaneeckhaute et al., , 2014Vlaco, 2012;Zaher et al., 2009) was used for simulation under variable operating conditions. To this end, a PHREEQC input script was developed for each unit, involving the identification of the input waste flows (PHREEQC data blocks: 'SOLUTION' and/or 'GAS'). A batch reaction calculation was also coded in case there is both a gas and liquid input, i.e. for the stripper and scrubbing unit (PHREEQC data block: 'REACTION'). Then, one factor at a time was allowed to increase within its range (e.g., PHREEQC code: REACTION_TEMPERATURE 20.0e70.0 in 51 steps), while the other factors were kept fixed. As such, a broad range of input scenarios was screened. Note that currently no alternative strategy is available in PHREEQC for selection of the various simulation scenarios (Parkhurst D. personal communication 2014). Yet, the development of an adequate, but more timeefficient, procedure to go through a multidimensional set of factors will be aspect of further research.
B. Run PHREEQC: Speciation calculations in PHREEQC/MINTEQ are made using designated thermodynamic databases that include a wide range of data for mineral phases and compounds. The calculations are based on three types of equations: 1) equilibrium relationships, 2) concentration conditions or mass balances (one per component), and 3) electro-neutrality conditions or charge balances (Chapra, 2008;Stumm and Morgan, 1996). By inclusion of oxidation/reduction reactions in the database, also the components' redox states were defined in the speciation calculations. The pH may be defined or adjusted according to the charge balance. The Davies equation was selected for ion activity correction in the NRMs, similar to Ali and Schneider (2008), Galbraith et al. (2014), Lizarralde et al. (2015), Ohlinger et al. (1998) and Flores-Alsina et al. (2016). The Davies ion activity correction was also recommended by Hafner and Bisogni (2009) above other relevant approaches, such as the Pitzer ion interaction approach. Moreover, the Peng-Robinson equation of state, which corrects for the non-ideal behaviour of gases, was used for calculating partial pressures ðpÞ and solubility (Parkhurst and Appelo, 2013). Furthermore, the temperature dependency of the thermodynamic equilibrium coefficients was expressed by means of the Van't Hoff relationship (Zumdahl, 2005), while the value of the water dissociation constant (K w ) at different temperatures (other than 25 C) was computed using the equation of Harned and Hamer (1933).
C/D. Selection criteria þ verification: From the speciation calculations the distribution of aqueous species (¼ ion activities) and saturation indices (SI) for phases (¼ driving forces for precipitation and gas transfer) were obtained. Soluble species with an insignificantly low activity, i.e. less than 0.01% of the total component activity in all scenarios, were excluded from the NRMs. Solids that may potentially precipitate (SI ! j0j) as well as gases that may volatilize (partial pressure (p) > 0) in the different units were selected. Conditions (pH, temperature) and rates for precipitation of the various forms of the selected minerals were also researched in the literature. The aim was to confirm the exclusion of the selected insignificant species and precipitates, while further identifying potential species and reactions that should be included in the database for each unit. The number of species and reactions that were found to be relevant according the speciation calculations and that were included in each NRM are presented in Table 4a. The list of species involved and the transformation reactions included in each model are presented in Appendix 1 (Table A1.1 and  Tables A1.2e1.6, respectively).

2.2.1.4.
Step 4 -building of a reduced model. Knowing that the generic geochemical model databases contain more than 3000 species (Allison et al., 1991), it was expected that the elimination of irrelevant species and reactions can have a significant impact on the simulation speed. As such, with the purpose of reducing model complexity and simulation times when coupling PHREEQC for NRM simulation (Section 2.3), a new PHREEQC database file including only the selected reactions and species was set up for each unit process. Moreover, a 'SELECTED_OUTPUT' data block was coded in the input script for each unit in order to transcribe only the appointed species and driving forces to the resulting output file. The latter is required for efficient coupling of the selected outputs to the kinetic and mass balance model (Section 2.3).
Finally, simulation results and speeds using the reduced model were compared with results and speeds obtained by running the developed chemical speciation scripts using the full Phreeqc.dat (P) and minteq.v4.dat (M) databases available in the PHREEQC 3.0.6 release.

Physicochemical transformation model
Heterogeneous physicochemical reactions, such as liquid-gas transfer and precipitation, occur much slower than the homogeneous reactions involved in the speciation calculations presented above. Hence, a kinetic approach was applied in order to allow for dynamic variation of the constituents.
Gas exchange processes in resource recovery systems can occur passively, i.e. without intensive gas bubbling (NRM-AD), or actively, i.e. with gas bubbling driven by an external air flow (NRM-Strip, NRM-Scrub). In each case similar kinetic gas exchange formulations based on the concentration driving force between the liquid and gas phases apply (Eq. (2)): where S liq;i is the liquid phase activity of component i ½M L À3 , p gas;i is the partial pressure in the gas phase of component i (atm), H T;i is the temperature-dependent Henry coefficient ½M L À3 atm À1 , H T;i :p gas;i represents the saturation concentration of gas component i in the liquid, K L=G;i is the overall liquid-gas mass transfer coefficient ½L T À1 , and a is the specific surface of the gas bubbles per reactor volume ½L À1 . Temperature dependency of H was described by a Van't Hoff relationship (Powers et al., 1987), while temperature dependency of K L=G;i a was described using the Arrhenius equation (Chapra, 2008). Through the coupling with PHREEQC (Section 2.3.1), both S liq;i and p gas;i can be calculated at every time step during the simulations. The total gas phase pressure was computed using Dalton's law of partial pressures (Stumm and Morgan, 1996). For calculation of K L=G;i a, a distinction was made between active and passive systems, since the values may differ significantly in practice (Chapra, 2008;Sotemann et al., 2006;Tchobanoglous et al., 2003). Moreover, a second distinction was made depending on the solubility of the gas considered, which determines whether mass transfer is liquid film controlled (for low to moderate soluble gases: H > 0.55, i.e. for CH 4 , CO 2 , H 2 , H 2 S, N 2 , O 2 ¼ all gases considered in the NRMs, except for NH 3 ) or gas film controlled (for very soluble gases: H < 0.55, e.g., for NH 3 ). As such, four potential mass transfer scenarios were considered, which are described in detail in Appendix 2: 1) Active liquid-gas/gas-liquid transfer (NRM-Strip, NRM-Scrub) of low to moderately soluble gases; 2) Active liquidgas/gas-liquid transfer (NRM-Strip, NRM-Scrub) of very soluble gases; 3) Passive liquid-gas/gas-liquid transfer (NRM-AD) of low to moderately soluble gases; 4) Passive liquid-gas/gas-liquid transfer (NRM-AD) of very soluble gases. The kinetic liquid-solid/solid-liquid transfer mechanisms described in all NRMs are nucleation (¼ birth of crystals), crystal growth, and redissolution. All reactions were represented by an empirical power law (Eq. (3)) using relative supersaturation (S À 1) as driving force (Ali and Schneider, 2008;Galbraith et al., 2014;Harrison et al., 2011;: in which S is the saturation ratio (¼

IAP Ks
, v refers to the stoichiometric precipitation coefficient which represents the total number of species involved in the precipitation reaction, IAP is the ion activity product ½M L À3 , K s is the solubility product ½M L À3 , k T is the temperature dependent transfer coefficient ½M L À3 T À1 , and n is the reaction order. The value of S was directly derived from the saturation index, ¼ log IAP Ks , which is calculated by PHREEQC at every time step during model simulations. The temperature dependency of the reaction rate was modelled by means of the Arrhenius equation (Greenberg and Tomson, 1992;. Using literature values for the molecular weight (MW) and density of the different precipitates, the total volume ðV fertilizer Þ, total mass/moles ðM fertilizer Þ; and MW ðMW fertilizer Þ of the recovered fertilizer product (composed of the various precipitates) was calculated at every time step. The time-dependent number of particles ðN part Þ was then determined using the Avogadro constant (N A ¼ 6.022 Â 10 23 mol -1 ). The mean particle diameter ðd p Þ of the precipitates was calculated assuming spherical particles using Equation (4): The kinetic precipitation/dissolution coefficient k T and the reaction order n in Equation (3) were adjusted according to the liquid-solid/solid-liquid transfer mechanism occurring: k G;T and n G for growth, k B;T and n B for nucleation, k D;T and n D for dissolution. The prevalent mechanism depends on the value of S and the amount of seed material in the reactor. Hence, these values were checked at every time step. As such, four possible scenarios were considered, which are described in detail in Appendix 2: 1) Supersaturation occurs (S > 1; SI > 0) and seed material is available; 2) Supersaturation occurs (S > 1; SI > 0), but no seed material is available and/or the crystal size is not large enough to have any influence on the process, i.e. the induction time is not exceeded; 3) The solution is undersaturated (S < 1; SI < 0) and precipitate is present in the system; 4) Equilibrium occurs (S ¼ 1; SI ¼ 0).
Finally, for the NRM-Prec, a generic mechanism for agglomeration and floc break-up through the effect of mixing was included using the spherical particle model for macroscale flocculation (Crittenden et al., 2012, Appendix 2). A time-dependent agglomerate number balance was also provided (Section 2.2.4). By division of the total fertilizer volume by the number of agglomerates, the agglomerate volume was obtained. The mean agglomerate diameter can then be computed in the same way as the particle diameter (Eq. (4)).
It should be remarked that mixing energy may also have to be included in Equation (3). Growth can be assumed surface integrated controlled when the system is well mixed, so the mixing effect can be neglected for the growth equations in unit processes with proper mixing Rahaman et al., 2014). However, mixing may affect the nucleation mechanism and induction time through microscale flocculation (Ohlinger et al., 1998). This mechanism is very site and species specific, hence it was considered out of the scope of the present generic model development study. However, by selecting a generic empirical equation based on S (Eq. (3)), the models could easily be extended to include mixing effects Perez et al., 2008), if appropriate parameter correlations are available. As mentioned above, future extensions may also involve particle size distributions (PSDs) (Nopens et al., 2014;Perez et al., 2008).

Biochemical transformation model
Biochemical processes and state variables are clearly important for the NRM-AD model. The description, stoichiometry, and kinetics of biochemical transformations that may be expected in the NRM-AD were based on the Anaerobic Digestion Model No. 1 (ADM1; Batstone et al., 2002), resulting in a total of 19 processes (Appendix 3: Table A3.1). pH, H 2 , and NH 3 inhibition expressions were taken from Batstone et al. (2002). Over the last ten years, various WRRF modellers (e.g. Flores-Alsina et al., 2016;Mbamba et al., 2016;Solon et al., 2015;Wang et al., 2016) have developed extensions of ADM1, mainly focused on the inclusion of a limited selection of chemical species and reactions to predict unwanted struvite precipitation and S inhibition in the digester. Since pH plays a critical role in anaerobic digestion modelling (Batstone et al., 2012;Solon et al., 2015;Vaneeckhaute et al., 2017b), inclusion of a more accurate and complete chemical speciation calculation, with associated efficient numerical solution procedure, to predict pH and driving forces for physicochemical and biochemical transformations is highly relevant.
In this study, ADM1 was for the first time extended with all essential physicochemical components and processes (acid-base reactions, ion pairing, liquid-solid transfer, liquid-gas transfer, redox transformations) that significantly impact anaerobic digester performance and digestate quality, selected in Section 2.2.1 (Appendix 3: Table A3.1, Extension 1). Ion pairing of cations with VFAs was also accounted for. On top of being important for predicting anaerobic digestion pH and performance, inclusion of such detailed physicochemical framework is essential for predicting process performance and product quality of physicochemical nutrient recovery unit processes that follow the digester. Indeed, the output digestate characteristics from the anaerobic digestion model should be sufficiently specified and compatible with the required input to the nutrient recovery unit process models. Similar as in Lizarralde et al. (2010), biological sulfate reduction (¼ sulfurgenesis) was incorporated based on the model proposed by Knobel and Lewis (Appendix 3: Table A3.1, Extension 2). An inhibition term for H 2 S was incorporated in the appropriate biokinetics (I H2S ), and its transfer to the gas phase was included as described in Section 2.2.2. The decay of SRBs was included in the same way as the decay of other organisms described in the ADM1 model . N, P, K, and S release from biomass, as well as nutrient uptake by growing biomass was accounted for (Appendix 3: Table A3.1, Extension 3). Modelling of EBPR sludge was considered beyond the scope of this study (Section 2.1.1), but for future applications the NRM-AD could be further extended using equations from, e.g., Ikumi (2011) Table A3.1, Potential Extension 4). Finally, N and P release through disintegration of complex particulates, P release from lipid hydrolysis, N release from protein degradation and amino acid uptake, as well as the N and P content of soluble and particulate inerts were also included. The detailed stoichiometric matrix and kinetic transformation equations proposed can be found in Appendix 3 (Tables A3.2-A3.4).
In this study, the biological solids leaving the digester were supposed to end up mainly in the solid fraction after solid-liquid separation of the digestate. Hence, in the subsequent key units for nutrient recovery, it was assumed that biochemical particulate transformations do not play a significant role. Nevertheless, in order to allow coupling of NRMs to activated sludge models (ASMs) in a later stage (through return liquors, for instance), the biological state variables were integrated in all NRMs. Note that the physicochemical interactions with the remaining soluble COD components were included in all models (Section 2.2.1).

Reactor model
The used reactor design and the default specifications and features for each unit process are compiled in Appendix 4. For each unit process, a mass balance was written, not only for all components in the liquid phase (S), e.g., Equation (5), but also for all components in the gas phase (G), all precipitated components (P), and all particulate biological solids (X), including both a transport term (based on in-and outgoing flow rates) and a transformation term (involving liquid-gas/gas-liquid transfer, liquid-solid/solidliquid transfer, and biochemical transformations): where P j¼1:n r j $v i;j is the summation of the specific kinetic process rates for process j ðr j ; ½M L À3 T À1 Þ multiplied by the stoichiometric coefficient for component i on process j ðv i;j; ½M M À1 Þ, Q liqin and Q liqout are the in-and outgoing liquid flow rates ½L 3 T À1 , V liq is the bulk reactor volume ½L 3 , and S liqin;i and S liq;i refer to the activities of the in-and outgoing liquid components ½M L À3 .
In addition, a mass balance for the seed material in the reactor was included, similar as Equation (5). The mass of seed material was adjusted in time according to the mass of precipitates present in the reactor and the liquid volume. Hence, it was assumed that newly formed crystals act as seed material for precipitation, similar as was experimentally discovered by Le Corre et al. (2007). External seed material can also be added.
For the precipitation unit (NRM-Prec), also particle and agglomerate number balances were implemented. The number of free precipitated particles was assumed to reduce according to the agglomerates formed, as in Crittenden et al. (2012). Note that agglomeration was only accounted for when mixing is present in the reactor (Section 2.2.2).

2.3.
Step III: model implementation and numerical solution

Model coding and state vector definition
The main coding language used in this study was Modelica, which is a high-level, declarative, and object-oriented modelling language (Claeys et al., 2006;Elmqvist et al., 1999). It is similar to the model specification language (MSL), which is currently used in Tornado/WEST (mikebydhi.com; Vanhooren et al., 2003), one of the most common software packages used in waste (water) quality modelling. However, Modelica has a better readability and expressiveness, and because of the more important industrial use (Audi, Ford, Siemens, etc.) of Modelica compared to MSL, the modelling community using Tornado/WEST recently decided to convert all conventional models for waste (water) treatment from MSL to the more powerful and more widely supported Modelica coding language. Tornado/WEST supports the use of both MSL and Modelica languages (Claeys et al., 2006).
As mentioned above (Section 2.2.1), a PHREEQC script was written for each unit process separately in order to include water chemistry. A 'SELECTED_OUTPUT' statement involving the selected species activities, saturation indices (SI's), partial pressures (p's), as well as the pH, temperature, alkalinity, and ionic strength was defined. The obtained SI's and p's are then used as driving forces for precipitation and gas transfer in the Modelica code describing the slow transformation processes (Eqs. (2) and (3)).
Since only small differences exist between the selected components for the different NRMs (Table 2), it was decided to define one generic component state vector for each different phase. As such, five different NRM component state vectors were enumerated (Appendix 5: Table A5.1): 1) Components_S1: the components in the liquid phase, i.e. the main waste flow; 2) Components_S2: the components in the H 2 SO 4 -solution used in NRM-Scrub; 3) Com-ponents_G: the components in the gas phase; 4) Components_P: the components in the precipitated phase; 5) Components_X: the particulate biological solids. The Components_S1 state vector was further split into a Components_S1_PC and a Components_S1_Bio state vector in order to describe physicochemical transformations and biological COD transformations separately. All state variable quantities involved in the physicochemical calculations (Compo-nents_G, Components_P, Components_S1_PC) were expressed on a molar base, whereas the state variables only involved in biological transformations (Components_X, Components_S1_Bio) were expressed on a COD-base. Moreover, for each model separately, a species state vector was enumerated referring to the PHREEQC selected output, which is different for each unit process.
Parameters and equations for the (slow) physicochemical and biochemical transformations, and mass balances for all total components were implemented in Modelica using a multi-matrix structure. The Tableau method matrix implementation of Morel and Herring (1993) was used as generic method for linking total soluble component activities to species activities and total precipitated component concentrations to precipitate concentrations in the NRMs, whereas the Gujer (2008) matrix implementation was used to describe the biochemical reactions involved.

Numerical solution and model execution procedure
To overcome problems related to the numerical solution of stiff systems (see Section 1), the slower reactions (Sections 2.2.2e2.2.4) and mass balances (Section 2.2.5) were represented by ordinary differential equations (ODE) coded in Modelica, while the fast reactions (Section 2.2.1) were assumed to reach steady state instantaneously and were calculated algebraically by use of algebraic equations (AE) at each iteration step using the software tool PHREEQC (Parkhurst and Appelo, 2013). In contrast to other WRRF modellers (e.g., Flores-Alsina et al.,2016;Hauduc et al.,2015;Mbamba et al., 2016;Tak acs et al., 2006) that implemented their own water chemistry module, the use of PREEQC to solve water chemistry was brought forward in this study (see Section 2; Vaneeckhaute et al., 2017b). PHREEQC has a dedicated and proven solver (Newton Raphson-based) for the complex set of implicit non-linear equilibrium equations involved. PHREEQC was preferred over other geochemical models (e.g., MINTEQ, WHAM, and WATEQ4F), because of its ease of integration with diverse scripting languages and other model libraries, next to its more precise methodology for precipitation calculations . Recently, a C-callable API (Application Programming Interface) for the PHREEQC engine has become available under the name IPhreeqc. It allows for easily coupling of the PHREEQC engine to software developed in other programming languages. The API provides direct access to the geochemical processes in the PHREEQC library, as well as support for new PHREEQC specification keywords that allow for easier manipulation of PHREEQC input and output data .
The models coded in the Modelica language, with invocations of the PHREEQC engine for speciation calculation, were then executed through the Tornado/WEST framework for modelling and virtual experimentation on the basis of sets of complex ODEs and AEs. A generic mechanism for calling PHREEQC from Modelica-specified models using Tornado was developed (Fig. 4). It consists of a Tornado-specific PHREEQC wrapper library containing only a predefined set of methods to be used in Tornado/WEST, as well as a reduced PHREEQC database and a PHREEQC script with selected outputs (Section 2.2.1). Any PHREEQC code can now be run, using input data supplied by Tornado/WEST and providing output data to be used by Tornado/WEST, in a flexible manner without the need for any case-specific C/Cþþ code modifications by the user. As a result, the combined kinetic-equilibrium models can now be used for simulation and other tasks such as parameter estimation, optimization, scenario analysis, Monte Carlo simulation, sensitivity analysis, and steady-state analysis, through the Tornado CUI (Command-line User Interface) tool, the user-friendly Tornado Experimenter GUI (Graphical User Interface), or WEST (Fig. 4).
Finally, for numerical solution in Tornado, two different solvers, RK4ASC (Runge Kutta 4 Adaptive Step size Control integration algorithm; Press et al., 1992) and VODE (Variable-coefficient Ordinary Differential Equation solver; Brown et al., 1989), were compared. The RK4ASC algorithm was retained, since simulation times were much faster and results more stable. This is likely related to its higher ability to solve models with certain discontinuities (i.e. sharp switches in behaviour, e.g., transitions in precipitation mechanisms as function of the saturation index) and dynamic inputs/disturbances (Claeys, 2008).

PHREEQC-Tornado interface
In order to connect state vectors used by PHREEQC (C code) and Tornado (Modelica code), a PHREEQC-Tornado interface was developed (Fig. 5). The interface makes special use of the data defined by the 'SELECTED_OUTPUT' data blocks (Section 2.3.1), and allows this array of data to be returned to Tornado without the necessity to read or write files. Hence, the data can be transferred between PHREEQC and Tornado through internal computer memory. This method of tight model coupling has significant merits with respect to calculation time and programming: a PHREEQC instance is only created once and is subsequently reused, preserving its internal state. In general, an order of magnitude decrease in run times is obtained compared to a loosely-coupled model, which requires starting PHREEQC as an external process for each time step (Müller et al., 2011). On top of that comes the gain in simulation time by using the developed reduced PHREEQC databases and scripts instead of full PHREEQC (Section 2.2.1). Hence, a reduction of execution time is obtained at two critical points during model simulations: i) the uploading and reading of database and input files, and ii) the transfer of data between PHREEQC and Tornado.

Model verification and debugging
After implementation, the models were subjected to a battery of tests to ensure implementation correctness, also referred to as model verification (Dochain and Vanrolleghem, 2001). A generic six-step procedure for model verification of NRMs was developed and applied to each unit process separately: 1. Verification of the PHREEQC-Tornado interface: Comparison of speciation calculations in Tornado through tight coupling to reduced PHREEQC, with simulation results from the independent full PHREEQC engine; 2. Verification of the physicochemical transformation model: Implementation of slow physicochemical transformations in Modelica code, execution in Tornado, and mass balance check; 3. Verification of the biochemical transformation model: Implementation of slow biochemical reactions in Modelica code, execution in Tornado, and i) mass balance check, ii) check against independent implementations, e.g., ADM1   As such, typing errors, inconsistencies, gaps, and conceptual errors were eliminated, while software bugs were discovered and dealt with.

Step IV: dataset collection and identification of data needs
One of the issues in the development of new models is the necessity to provide data for the estimation of model parameters and for the input variables. The different types of data required for each key NRM and the datasets that were used are provided in Appendix 6 (Table A6.1).
First, a thorough review of literature and existing models was conducted to provide default values for the different parameters involved (Appendix 6: Table A6.2-A6.5). Physicochemical stoichiometry and thermodynamic parameters are incorporated in the PHREEQC and Visual MINTEQ modelling software, where they are mainly taken from the National Institute of Standards and Technology (NIST, 2001) database. Default values for the kinetic precipitation coefficients were taken from literature, while default values for biomass kinetic coefficients were taken from the ADM1 model , except for the SRB kinetics for which the parameters were taken from Knobel and Lewis (2002) and Lizarralde et al. (2010).
Next to literature studies, also new experimental data aiming at NRM validation were collected through lab/pilot-scale testing and contact with industry. For NRM-AD, full-scale data at steady state from an anaerobic reactor treating S-rich paper mill primary sludge located at the WRRF Holmen Paper, Madrid, Spain has been obtained from the Center of Studies and Technical Research (CEIT, San Sebastian, Spain; Appendix 6: Table A6.6). An input fractionation was conducted following the procedure proposed by Grau et al. (2007).
For validation of the NRM-Prec, lab tests were conducted for P recovery from digestate under different operating conditions, i.e. different Mg:P-ratios and contact time notably. For this purpose, two different digestates were sampled at the full-scale biogas plants of SAP Eneco Energy (Houthulst, Belgium) and Wittevrongel Eneco Energy (Aalter, Belgium), which both treat agricultural wastes, mainly manure. A detailed input characterization was performed prior to the experiment (Appendix 6: Table A6.7). The precipitate was separated from the effluent by means of a centrifuge (5 min at 2000 rpm; Heraeus megafuge 1.0, Kendro Laboratory Products, Hanau, Germany), after which both fractions were also physicochemically analyzed. The P recovery efficiency (%) was then Fig. 5. Overview of the PHREEQC-Tornado interface coupling chemical speciation calculations to slow physicochemical and biochemical dynamic transformations at every time step. AE ¼ algebraic equations; ODE ¼ ordinary differential equations; X(0) ¼ initial state of the system; X(t) ¼ state of the system at time t. calculated using the P recovery of a control (no Mg addition) as a reference. For detailed methodology and experimental results, reference is made to De Corte (2012).
To obtain data for the NRM-Strip/NRM-Scrub, a technical and financial survey for a case treating 2000 m 3 d À1 of digestate at 200 mol NH 4 -N m À3 (more details: Appendix 6: Table A6.8) was carried out with various key suppliers in the field. As such, insights in the variability of the processes available to date were obtained, e.g., different target ammonium sulfate concentrations, operational pH and temperature, consumables, among other. These detailed data provided by the suppliers were used for further model refining and validation.
Finally, it should be stated that during model development new data needs appeared for which to date literature references are lacking. Such data gaps were identified and recommendations for future experiments and data collection are provided further in this paper (Section 3.5).

Step V: model validation
Model validation was performed in four different ways: i) validation against prior knowledge, ii) validation against existing models, iii) validation against literature or technical inquiries, and iv) validation against collected experimental results. In all cases, the default stoichiometric and kinetic parameter values determined in Section 2.4 were used. Input waste stream compositions, design data, and operational conditions were taken from the dataset involved. During the validation procedure, attention was given to the reduced PHREEQC database used in order to assure that all required species and reactions are included in the calculations. If required, an additional evaluation was conducted using the full PHREEQC and/or MINTEQ database, and missing species/reactions were additionally added to the reduced database.

Step VI: scenario analyses and process optimization
To gain more insight into the results and to further explore the model outcomes, scenario analyses were performed in Tornado/ WEST (Claeys, 2008). Moreover, the applicability of the models for process optimization was demonstrated by running optimization experiments in Tornado/WEST (Claeys, 2008).

Results and discussion
The implementation of the models developed in Section 2 was verified and validated. First, simulation times for the reduced models are evaluated in Section 3.1. General verification results and verification examples showing the correctness of the PHREEQC-Tornado interface are presented in Section 3.2. An example of model validation against experimental results, including scenario analyses and/or process optimization, is given for each NRM in Sections 3.3e3.5. Finally, recommendations for further research are provided in Section 3.6.

Evaluation of reduced PHREEQC simulation times
A comparison of simulation times of the developed scripts for each unit process model using the full databases and the corresponding reduced database is presented in Table 4b.
A 3-to 5-fold average improvement of model simulation speeds was obtained using the reduced database as compared to full Phreeqc.dat and minteq.v4.dat, respectively. The observed deviation in simulation times between PHREEQC and MINTEQ shows again the higher completeness of the MINTEQ database. Note that the presented simulation times in Table 4b concern the chemical speciation model only, so without the coupling to the kinetic and mass balance model. Yet, this model reduction is clearly relevant for simulation of WRRFs, since the speciation model is run at every time step during NRM model simulations (Section 2.3.3). As such, running a complete digestate treatment train under dynamic conditions for one year would take approximately 15 min (depending on the operating conditions and input characterization) using the reduced PHREEQC model, whereas it would take 45 min using the full PHREEQC model, both with tight model coupling (Section 2.3.3) to the kinetic model. As mentioned above, it is important for model validation to keep in mind that a model reduction was performed. As such, for example, it was discovered during initial validation of the NRM-Prec model that the species NaH 2 PO 4 (aq) was lacking, though essential for correct prediction of P recovery (Section 3.4.1).

General results and issues
During model verification, various software bugs were discovered and communicated to DHI, Merelbeke, Belgium, who successfully resolved the issues. As such, this research also contributed to the development of the Tornado/WEST software kernel.
Each step in the verification procedure was completed successfully. First, the PHREEQC-Tornado interface was found to be effective (see Section 3.2.2). Next, the mass balance check provided good results for each NRM. The step-by-step comparison of the Gujer matrix with other digester model implementations showed that the biochemical reactions were correctly implemented. Tests performed to check the ability of the models to realistically respond to model inputs, both under steady state and dynamic conditions, allowed eliminating small implementation errors. Some examples of tests and effects performed for model verification/validation can be found in Appendix 7. Finally, simulation results obtained from the two different independent implementations of the same model for each unit process, i.e. one using individual equations and one using a multi-matrix structure, were identical.
During model verification, three important general issues were observed to which future WRRF model developers must pay attention. First, it was found that some components, species, and precipitates that are highly important for modelling of WRRFs are not yet included in the generic PHREEQC and/or MINTEQ databases (Section 2.2.1: Table 3). Hence, for each new nutrient recovery model, the chemical speciation calculation should be verified with Table 4b Simulation times (s) and speed-up factor using the reduced PHREEQC database as compared to the full Phreeqc.dat (P) /minteq.v4.dat (M) databases for simulation of the chemical speciation scripts developed for each key unit in the nutrient recovery model (NRM) library. AD ¼ anaerobic digestion; Prec ¼ precipitation/crystallization; Strip ¼ stripping; Scrub ¼ scrubbing. multiple software packages, with literature, and with prior knowledge in order to comprehensively select which components, species, and precipitates should be included in the model and which ones can be excluded. Secondly, if an input to PHREEQC is set to 0 or if a species is not defined or not present in the calculation, then a value of À999.999 is printed as output for this component's species distribution and the corresponding saturation indices and partial pressures. In the Modelica code, these outputs are then used as driving forces for slow transformations, leading to incorrect calculations. This issue was solved by introduction of an if-thenelse statement in the PHREEQC-Tornado interface. Finally, attention should be paid to the use of units for input and output variables. Input concentrations in PHREEQC are expressed by default as mole m À3 , whereas the outputs are given by default as kmole m À3 .
Deviations from these standard units should be declared in the PHREEQC script.

Verification of PHREEQC-Tornado interface
When comparing simulation results using the stand-alone full PHREEQC engine and Tornado (with tight coupling to reduced PHREEQC), identical model outputs were obtained for all NRMs. As an example, the results for the NRM-Scrub are given in Table 5. An initial gas phase flow with high NH 3 load (coming from the NRM-Strip) was used as input to the NRM-Scrub and brought into contact with a sulfuric acid solution for NH 3 absorption. The outputs, i.e. the logarithm of the partial pressures (log(p), atm) in the purified gas phase and the activities (mole m À3 ) of some species in the ammonium sulfate solution after gas-liquid exchange, obtained with both the stand-alone PHREEQC engine and Tornado-PHREEQC are presented. It can be concluded that the implementation of the PHREEQC-Tornado interface and the PHREEQC invocation in Modelica are correct.

NRM-AD validation
3.3.1. Case study anaerobic reactor at Holmen Paper Madrid (Spain) The NRM-AD model was validated using experimental data collected under steady state conditions from an anaerobic digester for treatment of S-rich paper mill primary sludge at a full-scale WRRF (Holmen Paper, Madrid, Spain). The same case was previously used for validation of the Lizarralde et al. (2010) model for anaerobic S reduction. The input sludge characteristics, design parameters, initial reactor state variables, and operating conditions are given in Appendix 6 (Table A6.6). Kinetic and stoichiometric parameters were set at default (Section 2.4). A comparison of experimental and simulation results using the NRM-AD and the model proposed by Lizarralde et al. (2010) is given in Table 6.
Simulation results using the NRM-AD show good agreement with the experimental results for COD removal and biogas CH 4 and CO 2 composition at a particular pH. The model also seems to give a very good prediction of the digestate pH and P content, and a relatively good prediction for NH 4 -N in the digestate. The slightly higher digestate nutrient value for NH 4 -N obtained with the NRM-AD may be attributed to losses of NH 3 during digestate sampling and analysis, although potential model deficiencies may not be excluded.
The NRM-AD seems to underpredict the biological SO 4 removal and corresponding H 2 S production by SRBs, as will be explored below. However, from a pure validation perspective (note: no parameters were calibrated!), when comparing with the Lizarralde et al. (2010) model, overall the performance of the NRM-AD is significantly better, very probably due to the underlying more detailed chemical speciation and the inclusion of multiple competing physicochemical transformation reactions.

Exploration of hypotheses regarding S cycle measurements
Through model scenario analyses, four potential hypotheses were tested to explore the underestimation of biological SO 4 removal in the above case study. First, it was observed that the biogas H 2 S concentration was very sensitive to variations in pH (cfr. Al-Zuhair et al., 2008). Model simulations were carried out at the digestate pH (7.21). However, the input pH was significantly lower (6.66) and the digestate pH may be influenced through contact with air. Hence, there exists some uncertainty about the actual reactor pH.
To explore this hypothesis, a scenario analysis was conducted in order to evaluate the effect of pH (variable) on the % CH 4 , CO 2 , and H 2 S in the biogas at fixed waste input COD:SO 4 -ratio. Assuming that the pH in the reactor ranged from 6.66 (waste input pH) to 7.21 (digestate pH), the biogas composition varied from 61% CH 4 , 34% CO 2 , 2.94% H 2 S to 80% CH 4 , 16% CO 2 , 1.90% H 2 S. Hence, with the present implementation, it was not possible to obtain 6% H 2 S in the biogas at a pH in that range.
It should be remarked that the experimentally obtained biogas H 2 S content of 6% is extremely high compared to literature values. Typical biogas H 2 S values for similar concentrated sulfurous streams from the paper industry range between 1 and 2% H 2 S (Reiter and Piccot, 2004). Hence, a second reason for the uncertainty may be related to the H 2 S analysis itself, conducted by the operators.
A third explanation may be the exclusion of lactate in the present NRM-AD implementation. Lactate is a preferred substrate for sulfate reducing bacteria and would thus aid in increasing SO 4 removal and H 2 S production (Oyekola et al., 2007). This may explain the slight overestimation of biogas CO 2 production and underestimation of H 2 S production. In the present case, no lactate measurements were available, but future research should consider this component.
Furthermore, the non-consideration of reactions (precipitation/ ion pairing) with Al and Fe, due to lack of input Al/Fe measurements at the WRRF, may explain the lower SO 4 removal found through simulation (cfr. Zhang et al., 2013). This can also explain why model predictions for COD removal and CH 4 production were good, while additional COD would be required for additional SO 4 removal by SRBs. Based on a similar reasoning, Lizarralde et al. (2010) assigned potential sulfate precipitation (which was not considered in their model) to the highly overestimated H 2 S production found with their model.
An attempt to calibrate input Al in the present case study showed that a reactor concentration of 276 mol Al m À3 resulted in a SO 4 removal of 78% (¼ experimental value in Table 6) and a biogas H 2 S concentration of 3%. However, in this scenario the pH lowered to a value of 6.26. The higher SO 4 removal found through addition of Al was likely the result of a combination of multiple effects. It was, for example, observed that the addition of Al affected the amount of Ca/Mg sulfates and Ca/Mg precipitates formed. The addition of Fe resulted in a lower H 2 S production because of FeS precipitation, but it did not aid in SO 4 removal. Finally, other model gaps can of course not be ruled out and one should bear in mind that the above validation is based on a onetime test.
It can be concluded that more detailed waste(water) input characterizations, including all selected components for the NRM-AD unit process (Section 2.2.1: Table 2), as well as instantaneous pH measurements in the reactor, are required in order to correctly calibrate the model for biological S removal. Nevertheless, clearly, exploration using the NRM-AD leads to increased insights and better understanding of the various interacting processes occurring in digesters.

Phosphorus precipitation at different Mg:P-ratios
For validation of the NRM-Prec model, batch experiments were carried out in the lab for P recovery from two different crude digestates (Section 2.4; Appendix 6: Table A6.7). Different Mg:Pratios obtained through dosing of MgCl 2 :6H 2 O were applied aiming at the production of N-struvite (MgNH 4 PO 4 :6H 2 O or MAP) or Kstruvite (MgKPO 4 :6H 2 O or MKP) fertilizer. Initial simulation results showed a large deviation from the experimental results (Table 7). After evaluation using the full PHREEQC and MINTEQ databases, this deviation could be attributed to ion pair formation of NaH 2 PO 4 , a species that was initially not included in the reduced PHREEQC database, nor in the generic PHREEQC database (Table 3). Indeed, due to the high Na concentration of both digestates, Na paired with P, making it less available for precipitation. When NaH 2 PO 4 was added as species to the reduced database, a very good agreement between the simulation and the experimental results was obtained for P recovery at steady state (after 12 h; Table 7).
This finding is in line with the results obtained by Li et al. (2012), who found a ± five times higher residual effluent P concentration when NaH 2 PO 4 þ MgCl 2 :6H 2 O were dosed for struvite precipitation, compared to the dosing of H 3 PO 4 þ MgCl 2 :6H 2 O. Moreover, recently Chauhan and Joshi (2014) found that at high Na:NH 4 -ratios, NaH 2 PO 4 is formed instead of, or next to, NH 4 H 2 PO 4 , the precursor for MAP precipitation. In turn, this compound may be transformed into Na-struvite through the following reaction: The formation of Na-struvite was not yet included in the NRM-Prec model due to lack of knowledge on the existence, the stoichiometry, and the kinetics of this precipitation reaction. However, knowing that current practice often involves the addition of NaOH for pH-increase prior to struvite crystallization, the case study above clearly shows the relevance of further research on Na-P ion pair formation and Na-struvite precipitation kinetics in waste(waters). The phenomenon may not only impact the effluent quality, but also the quality of the resulting recovered fertilizer product, i.e. a potential mixture of N/K-and Na-struvite may appear.

Exploration for process understanding and optimization
Two questions arise from the experimental and simulation results presented above (Table 7): 1. Why is the P recovery efficiency rather low for both digestates? 2. Why does increasing the Mg dose not improve the P recovery efficiency?
The ability of the models to find an answer to such questions is presented below.
First, it was observed experimentally and through simulations that the main precipitated components, next to P, were Al, Ca, Fe, K, Mg, and N(eIII). Hence, the product recovered was definitely not pure MAP or MKP. A scenario analysis evaluating these components was conducted for both digestates in order to obtain more insights in the results (Fig. 6). The two tested digestate compositions under study are marked as stars in Fig. 6.
The maximum achievable P recovery as function of the input Mg and Ca content was 56.2% for digestate 1 (Fig. 6A), whereas it amounted to 90.7% for digestate 2 (Fig. 6B). This discrepancy can be attributed to the higher concentration of Fe and Al in digestate 1  Fig. 6). Up to ±110 mol m À3 of input Ca (the margin in which the digestates under study are situated), mainly ion pairing of CaHPO 4 (aq) and CaPO 4 À was observed, which decreased the amount of P available for precipitation (cfr. Lin, 2012). Above a value of ±110 mol m À3 , calcium phosphates became oversaturated, precipitation occurred, and P recovery increased. This effect of Ca inhibition observed through model simulations is in agreement with the experimental findings of Huchzermeier and Wendong (2012). The latter concluded that struvite purity decreased because of the formation of calcium phosphates when the Ca:P activity ratio was greater than 0.5 to 1. Secondly, the fact that the P recovery efficiency in the presented experiment was not much influenced by increasing Mg:P-ratios can, according to the model, be attributed to the formation of dolomite (CaMg(CO 3 ) 2 ), as well as Mg(OH) 2 and Mg 2 (CO 3 )(OH) 2 :3H 2 O at higher Ca and Mg concentrations. Indeed, higher Ca and Mg doses are associated with a pH-increase, which favours carbonate and hydroxide precipitation (Zumdahl, 2005). When the input Ca concentration would be 0, one can see an increase in P recovery with increasing Mg dose due to the formation of MKP (note: high K-concentration in the input waste stream) and Mg-phosphates. This competitive effect between Mg, Ca, and P found through NRM-Prec simulations is in agreement with the findings of Lin (2012), who obtained a precipitate mixture of struvite, dolomite, Mg(OH) 2 , calcium phosphates, and CaCO 3 in experiments on P recovery from digested swine manure.
Based on the above-mentioned findings, two optimizations of the process can be proposed if the aim would be to produce high purity struvite: 1. Removal of CaCO 3 through precipitation prior to the experiment, e.g., using a filtration system as in Huchzermeier and Wendong (2012); 2. Elimination or reduction of the use of Fe and Al in the WRRF processes upstream of the precipitation unit, e.g., for improved sludge dewatering. This measure could also be assessed by locating the struvite precipitation unit upstream in the WRRF, e.g., immediately after the activated sludge (AS) system (cfr. combined use of the WASSTRIP and Pearl process for improved P release and struvite recovery; Cullen et al., 2013). In fact, the AS system itself could also (partially) be replaced by a strip/scrub system.
When applying these proposed measures in a treatment train for digestate 1, the maximum achievable P recovery through simulation became 91%, consisting of MKP, Mg(OH) 2 , and Mg 3 (PO 4 ) 2 . Hence, a pure Mg/P/K fertilizer would be obtained (Fig. 7). Remark that the main precipitate found, MKP, is not included in the standard PHREEQC/MINTEQ databases. Hence, again, the extensions provided to the database are clearly relevant (Section 2.2.1).
It should also be noted that in Fig. 7, the Mg dose was allowed to change within the range of 0e500 mol m À3 (so no point measurements are presented). Hence, the abrupt changes in slope are related to actual changes in precipitation mechanisms, which could, e.g., involve transitions from nucleation to particle growth and/or agglomeration and/or redissolution, and changes in precipitated species due to changes in saturation indices. Moreover, an interesting observation made through model simulations was that, without any addition of Mg, a high P recovery efficiency of 72% could be obtained. This could be appointed to the precipitation of K 2 NH 4 PO 4 :6H 2 O (¼ pure N/P/K fertilizer) due to the high amounts of K available in the digestate (Appendix 6: Table A6.7). In this case, an economic analysis is recommended to select a target fertilizer, thereby taking into account local fertilizer market demands, and environmental and fertilizer regulations. On the one hand, the use  of chemical Mg may increase the operational costs of P recovery, while on the other hand a higher recovery efficiency can be obtained, with larger mean particle diameter of the recovered precipitates (mainly MKP), as predicted with the NRM-Prec. Larger particles generally increase the revenues from fertilizer sales (Vaneeckhaute et al., 2017a).

NH 3 recovery at different operating conditions
During validation of the NRM-Strip and NRM-Scrub models, NH 3 stripping was found to be very sensitive to the total and relative input concentration of carbonates, Ca, and Na, since these determine the input alkalinity and pH. Since operators usually only measure NH 3 and pH (þsometimes total alkalinity), an identifiability problem arises. For example, when using the design parameters and input flow characterization (S_N_min3_, pH) of Collivignarelli et al. (1998), a good agreement was obtained between experimental and simulation results for NH 3 recovery (Table 8).
However, due to lack of some fundamental input flow characteristics for pH calculation using the NRM-Strip model, the input composition had to be calibrated in order to approximate the operational pH. Evidently, there are multiple ion combinations possible to obtain the specified pH, and the choice of the combination may influence the model outputs. Hence, in order to effectively use the NRM-Strip/NRM-Scrub models for process optimization, the initial waste flow composition should be characterized in more detail than is usually done at WRRFs today. Irrespective thereof, it can be seen in Table 8 that the model responded correctly to disturbances/operational decisions, such as an increase in pH, temperature, and air flow rate (cfr. Collivignarelli et al., 1998).

Treatment train for NH 3 recovery
In order to overcome the above-mentioned identifiability issue, a technical survey was sent out to key suppliers of strip/scrub units for the treatment of a particular digestate flow (Section 2.4). Using the predefined input characteristics (Appendix 6: Table A6.8), as well as the dimensions, operating conditions, effluent quality and stripping performance offered by the different suppliers, the models were again validated for the different configurations received. To this end, first a treatment train consisting of NRM-Chem, NRM-Strip, and NRM-Scrub was built to reflect a full-scale installation. Then, model simulations using the design data were conducted and scenario analyses were performed in order to compare simulation results with the experimental data obtained from the suppliers.
The data of the supplier who provided the most detailed experimental results (RVT Process Equipment GmbH) are presented below as an example. An NH 3 recovery efficiency of ±90% at 55 C was guaranteed by the supplier, when increasing the pH to a value of 10.3 by addition of 102.5 mol m À3 NaOH d À1 under the design conditions provided in Appendix 6 (Table A6.8). The same results were obtained through treatment train simulation (Table 9).
Finally, technology providers also advised to remove excess input carbonate buffer capacity prior to treatment, e.g., through CO 2 stripping, in order to minimize NaOH consumption for pH-increase as well as CaCO 3 precipitation in the reactor (P erez, 2002). This   (Collivignarelli et al., 1998)   recommendation could be confirmed using the NRM-Strip model: Fig. 8 shows the decreasing NH 3 recovery efficiency as function of carbonate buffer capacity, if the NaOH consumption and other operating conditions would not be adjusted. Hence, the more carbonate is stripped off, the higher the reactor pH and the higher the NH 3 recovery efficiency. Note that, based on this principle, some technology suppliers provide an integrated CO 2 and NH 3 stripping process without using NaOH for pH-increase (Vaneeckhaute et al., 2017a).

Recommendations for further experimental research
The results show that the performance of all resource recovery systems under study is very sensitive to the input waste stream composition, e.g., through its direct effect on the pH. Therefore, in order to obtain good model predictions for a particular waste flow, the input flow should be characterized in more detail than is usually done at WRRFs today. This observation is similar to activated sludge modelling in which influent characterization is considered as the most important step for achieving accurate results (Rieger et al., 2012). It is clear that a better characterization of the input composition may help to adjust the use of consumables (e.g., chemical dose and air requirements) to a minimum, thereby reducing the operational costs. As such, the models can be used as an invaluable tool for process optimization. New experimental results, including detailed input characterizations, are currently being collected at pilot/full scale under dynamic conditions in order to further calibrate and validate the proposed NRMs.
A second issue observed is that values for the precipitation kinetics (k T ) and gas transfer coefficients (K L=G a) used from literature are commonly determined under ideal conditions, i.e. gas transfer in clean water and precipitation in a synthetic solution containing only the target species involved in the reaction, e.g., Mg, NH 4 , and P for struvite precipitation. However, the actual value of these parameters may be highly influenced by the complex matrix of the waste streams involved, e.g., through ion pairing (Section 3.4.1), concurrent and competing precipitation reactions (Section 3.4.2), and the presence of seed material. Studies evaluating kinetic rates under actual process conditions are lacking in literature, but should be focus of further research in order to correctly calibrate these parameters in the NRMs. Moreover, rates and mechanisms for nucleation, agglomeration and dissolution of various precipitates are still unknown and should be further studied. In this perspective, the use of the simple empirical equation (Eq. (3)) for liquid-solid/ solid-liquid transfer in the NRMs is interesting compared to previously used approaches in wastewater treatment (e.g., Hauduc et al., 2015;Lizarralde et al., 2015;Musvoto et al., 2000).
Another important complication is related to the characterization of the precipitates formed. X-ray diffraction is the commonly used technique to characterize precipitates in pure solutions. However, it generally requires pure crystals of high regularity to solve the structure of a complicated arrangement of atoms. Also, the results usually represent a very local microstructure, and it requires a lot of work to obtain a certain statistical reliability on the results (Tanigawa et al., 2003). More research is required on the development of a generic and cost-effective experimental method to accurately characterize the different precipitated species from a complex waste matrix. Such a procedure may not only be used to determine the precipitated species in precipitation units and hence the recovered product purity , but also, for example, the precipitates in the digestate leaving the digester. The latter is relevant as these precipitates may act as seed material for precipitation downstream. Finally, interesting model extensions have been identified. They lead to the inclusion of: i) Lactate as specific substrate for biological sulfate removal in the NRM-AD, e.g., as in UCT (2007); ii) A transformer tool in the NRM-AD to allow for co-digestion of multiple input streams, e.g., the GISCOD tool (Zaher et al., 2009); iii) Biochemical transformations of EBPR sludge in the NRM-AD, e.g., as in Flores-Alsina et al. (2016), Ikumi (2011), or Wang et al. (2016; iv) Microscale flocculation in the NRM-Prec, e.g., as in Crittenden et al. (2012); v) Particle size distributions in the NRM-Prec, e.g., as in Perez et al. (2008); vi) Differential settling in the NRM-Settle and (if relevant) in the NRM-Prec, e.g., using the Stokes equation (Crittenden et al., 2012) or a particle settling velocity distribution (Bachis et al., 2015); vii) Heavy metals (and other contaminants) in all NRM models. These extensions will of course lead to further experimental data requirements.

Conclusions and future perspectives
The first available generic nutrient recovery model (NRM) library including dynamic mathematical models based on both detailed chemical solution speciation calculations, as well as physicochemical and biochemical reaction kinetics, was developed through PHREEQC-Tornado/WEST coupling, and successfully validated at steady state; Implementation correctness was verified under steady state and dynamic conditions using a 6-step procedure, including, e.g., a comparison of the simulation outputs of two independent model implementations for each unit process: one based on individual equations and one compact matrix-based implementation; Using a systematic procedure for PHREEQC database reduction a 3-to 5-fold improvement of model simulation speed was obtained as compared to the use of standard thermodynamic databases, on top of the improvement obtained through tight model coupling; Because of gaps in existing standard thermodynamic databases, an extended database with the purpose of nutrient recovery was made available, named 'Nutricover.dat'; Detailed input characterization was found to be most critical for accurate prediction of resource recovery process performance; Simulation results showed the high potential of the NRM library to increase understanding of nutrient recovery process interactions and to optimize integrated nutrient and energy recovery systems; The use of the NRM library by the various stakeholders in the field to facilitate the implementation, operation, and optimization of nutrient recovery technologies can stimulate the transition from WWTPs to more sustainable WRRFs.

APPENDIX 1
Physicochemical species and reactions included in the nutrient recovery models (NRM) Values between brackets represent the use of air instead of chemicals for pH-adjustment.
Values between brackets represent the use of air instead of chemicals for pH-adjustment.

APPENDIX 2
Detailed description of mass transfer scenarios for a) liquid-gas/gasliquid transfer and b) liquid-solid/solid-liquid transfer a) Liquid-gas/gas-liquid transfer If the resistance to mass transfer is on the liquid side, the overall liquid mass transfer coefficient, K L;i , can be perfectly adequate, while the overall gaseous mass transfer coefficient, K G;i , provides a good estimation if the resistance is on the gas side. The relationship between the two coefficients can be represented by Equation (A1) (Chapra, 2008;Tchobanoglous et al., 2003): in which R is the universal gas law constant (0.082 l atm mol À1 K À1 ) and T the temperature (K). It should be noted that the abovementionned overall mass transfer coefficients are actually derived from the individual mass transfer coefficients by Equation (A2) (combined with Eq. (A1) for K G;i ): in which k L;i and k G;i are the individual mass transfer coefficients that depend on the conditions at the interface and the bulk of the liquid and gas phase, respectively (Chapra, 2008;Tchobanoglous et al., 2003). Nevertheless, since the concentrations at the interface are difficult to measure, the overall mass transfer coefficient is generally used for practical purposes. As such, four potential mass transfer scenarios were considered: 1) Active liquid-gas/gas-liquid transfer (NRM-Strip, NRM-Scrub) of low to moderately soluble gases. In this case, the penetration theory of Higbie (1935) was used to calculate the liquid mass transfer coefficient, K L a ½T À1 . It states that diffusion is a non-steady state process and that the molecules of the solute are in constant random motion. Clusters of these molecules arrive at the interface, remain there for a fixed period of time, and some of them penetrate while the rest mixes back into the bulk of the phase. The transfer velocity was then formulated in terms of the average contact time of a gas bubble at the interface (Eq. (A3); Chapra, 2008;Gujer, 2008): in which d is the average gas bubble diameter (default ¼ 3 mm; Gujer, 2008), u is the rise velocity of the gas bubbles The obtained D l values using Equation (A4) showed good equivalence with D l values found in literature for wastewater systems (Chapra, 2008;Gujer, 2008;Tchobanoglous et al., 2003).
3) Passive liquid-gas/gas-liquid transfer (NRM-AD) of low to moderately soluble gases. In this case, the mass transfer rate needs to be calibrated based on experimental results, e.g. as in Tourlousse and Ahmad (2007), because the rise velocity of gas bubbles is usually not measurable or very difficult to measure. For convenience, the Table A1.6 Gas-liquid / liquid-gas exchange reactions (GL) included in each nutrient recovery model (NRM). AD ¼ anaerobic digestion; Prec ¼ precipitation/crystallization; Strip ¼ stripper; Scrub ¼ scrubber.

No.
Gas-liquid /liquid-gas exchange reaction AD Prec Strip Scrub CO 2 (aq) 4 CO 2 (g) a Values between brackets represent the use of air instead of chemicals for pHadjustment.
K L a is usually calculated from the K L a of oxygen gas (O 2 ) as a reference compound, since the latter is easy to deduce from experimental data and rate constants for volatile solutes can be assumed proportional to each other (Chapra, 2008;Ikumi, 2011;Mackay and Yeun, 1983;Munz and Roberts, 1989;Musvoto et al., 2000). However, the use of O 2 as a reference compound, as selected by Musvoto et al. (2000), is quite odd for anaerobic digestion, because normally no O 2 is present in such reactors. Therefore, in the NRM-AD model, H 2 was used as volatile reference compound occurring in digesters, similar as in (Pauss et al. (1990); Eq. (A5)): 4) Passive liquid-gas/gas-liquid transfer (NRM-AD) of very soluble gases. In this case, the mass transfer rate should be determined independently of the low to moderately soluble gases above (Sotemann et al., 2005). If no experimental data are available, the K G a value for NH 3 in anaerobic digestion is usually set to a very low value ranging from 1.92 to 3.2 d À1 (default in NRM-AD ¼ 3.2 d À1 ; Ikumi, 2011;Musvoto et al., 2000;Sotemann et al., 2005). This is to ensure an extremely low loss from the liquid phase through stripping. However, as the transfer rate depends much on design, operating conditions, and characteristics of the waste flow to be treated, it is advised to determine the K G;NH3 a under actual environmental conditions, as e.g. in Arogo et al. (1999). This is especially important for the stripper and scrubber unit processes.
b) liquid-solid/solid-liquid transfer 1) Supersaturation occurs (S > 1; SI > 0) and seed material is available. In this case, the crystallization of sparingly soluble salts in WRRFs is mainly controlled by surface spiral growth. This means that the integration of the cations into crystal lattice positions at kinks in the surface is the rate-determining molecular mechanism Hauduc et al., 2015;Koutsoukos et al., 1980;Musvoto et al., 2000;. The kinetic precipitation coefficient (Eq. (A6)) was then assumed to be proportional to the available seed material (cfr. Koutsoukos et al., 1980;Parkhurst and Appelo, 2013): in which k G;T is the temperature dependent growth rate coefficient ½M L À2 T À1 ; a seed is the specific area of surface per gram of seed material before the seed crystals start to grow in the crystallizing solution ½ L 2 M À1 (default ¼ 600 m 2 g À1 ; Parkhurst and Appelo, 2013), and M seed is the time-dependent mass of seed material in the reactor ½M (default initial mass ¼ 0.0005 kg; Parkhurst and Appelo, 2013). The latter is calculated at every time step by means of mass balances on the seed material for each precipitate (Section 2.2.4), taking in account the mass of newly formed precipitates and redissolution. The default reaction order for surface controlled growth (n G ) was set at 2, which generally provides a good approximation to represent precipitation in WRRFs (Bouropoulos and Koutsoukos, 2000;Mehta and Batstone, 2013;Musvoto et al., 2000;. 2) Supersaturation occurs (S > 1; SI > 0), but no seed material is available and/or the crystal size is not large enough to have any influence on the process, i.e. the induction time is not exceeded.
In this case, primary nucleation occurs, which was often not accounted for in previous studies Schneider et al., 2013), though very relevant . The value of k T and n in Equation (3) are then switched to the nucleation rate, k B;T (default ¼ 10 6 nuclei L À1 T À1 ; Mehta and Batstone, 2013), and the nucleation reaction order, n B . The latter is usually higher for nucleation than for growth (3e4; default ¼ 3; Tavare, 1995). The induction time is inversely proportional to the logarithm of S, and should be estimated experimentally for each precipitate Mehta and Batstone, 2013).
3) The solution is undersaturated (S < 1; SI < 0) and precipitate is present in the system. In this case, the NRMs allow for precipitate redissolution until equilibrium is reached using the reverse reaction of Equation (3) (Morse and Arvidson, 2002). However, the kinetic dissolution rateðk D;T Þ and the reaction order for dissolution ðn D Þ may be different than those for precipitation. Significantly more work is needed to better understand the dissolution behaviour of the various precipitates in complex waste(water) matrices (Greenberg and Tomson, 1992;Morse and Arvidson, 2002).
4) Equilibrium occurs (S ¼ 1; SI ¼ 0).. In this case, the liquid-solid /solid-liquid transfer rate is set at 0. Finally, for the NRM-Prec, a generic mechanism for agglomeration and floc break-up through the effect of mixing was included using the spherical particle model for macroscale flocculation (Crittenden et al., 2012). The net rate of floc appearance (Eq. (A7)) was written as: in which K a ½À is the aggregation constant (¼ 4a=p for laminar flow where a is the collision efficiency factor; default for turbulent flow ¼ 5 Â 10 À4 ), K b ½T dÀ1 :L À3 is the floc break-up constant (¼ 0 for laminar flow; default for turbulent flow ¼ 10 À7 ; Crittenden et al., 2012), G is the root mean square velocity gradient ½T À1 which depends on the power input (Camp and Stein, 1943), and d is the turbulence constant. Under turbulent conditions, the values of K a and K b should be determined empirically in laboratory or pilot-scale tests (Argaman, 1971;Parker et al., 1972). Note that when the G value is set to 0, it is assumed that no agglomeration occurs.

APPENDIX 3
Biochemical processes and Gujer matrix included in the nutrient recovery model for the anaerobic digester (NRM-AD) Table A3.2 Stoichiometry of the biochemical (BC) Gujer matrix incorporated in the nutrient recovery model for the anaerobic digester (NRM-AD). For process description: see    Sulfurgenesis (Knobel and Lewis, 2002;Lizarralde et al., 2010) Release/uptake of P, K, S from bacterial cells and other biochemical components EBPR sludge (Ikumi, 2011) NRM-AD Extension 1: Acid-base systems: Table A1.2 Redox reactions: Table A1.3 Ion pairing reactions: Table A1.4 Solid-liquid transfer: Table A1.5 Gas-liquid exchange:     Yield of biomass on substrate kg COD X kg À1 COD S -Variable volume as function of retained precipitant volume; -Precipitate flow rate (Q_prec) extracts fraction of the precipitates continuously or at specific times when selected specifications are reached, e.g. target particle diameter, purity, etc.; -Allows to study the effect of mixing power and reactor seeding on, e.g., the mean particle/aggregate diameter; -Optional: use of gas flow instead of chemicals for pH-increase in the reactor; -Potential extension: inclusion of particle (differential) settling velocity (Crittenden et al., 2012).

NRM-Strip
Stirred tank for active liquid-gas exchange (based on Gujer, 2008) -Continuous in-and outgoing liquid and gas flows; -Newly formed gas bubble enters the reactor at an initial gas phase concentration; -Model parameters averaged over all bubbles; -Heterogenous gas transfer throughout the reactor height; -User-selectable number of liquid layers to represent spatially dependent liquid transfer. a NRM-Scrub Stirred tank for active gas-liquid exchange (Gujer, 2008) -Similar specifications as NRM-Strip, but: i) Default use of sulfuric acid solution at pH 1.3 for NH 3 absorption; ii) AmS recycle flow (Q_rec) with extraction as fertilizer flow when user-selected AmS specifications (usually 25e40% AmS concentration) are reached (cfr. semi-batch process).
NRM-Chem Point mixer -Closed tank to avoid NH 3 emissions through pH-increase; -Allows addition of the most important amendments in WRRFs: i) MgCl 2 , Mg(OH) 2 , and Ca(OH) 2 prior to P precipitation, ii) NaOH and Ca(OH) 2 prior to stripping (goal ¼ pH [, CaCO 3 scaling Y); -Usually followed by NRM-Prec to allow for species precipitation and flocculation.

NRM-Heat
Point heater -Colder fluid gaining heat from a hot gas/steam flow or a hot liquid flow; -Generic equation based on the specific heat of the fluid, the surface area of the heat exchanger, and the overall heat transfer coefficient (AIC, 2014); -Application prior to NRM-AD and NRM-Strip. a Some literature studies show that hydraulic levels and reactor design have no effect on the NH 3 recovery efficiency as equilibrium conditions are reached in a very small time interval (Arogo et al., 1999;Collivignarelli et al., 1998;Gujer, 2008;Powers et al., 1987). However, other studies believe that liquid transfer should be modelled heterogeneously, i.e. spatially dependent (Yu et al., 2011). Because of this discussion, an option was included in the NRM-Strip and NRM-Scrub to calculate NH 3 removal and absorption for a user-selectable number of liquid layers. The Gujer (2008) model is based on homogeneous liquid transfer.

APPENDIX 4
Reactor design and the default specifications and features for each unit process

APPENDIX 5
State vectors used in the nutrient recovery model (NRM) library

NRM-AD:
Reducing the input alkalinity to the digester results in a (delayed) pH decrease (less carbonate buffer) because of volatile fatty acid accumulation. Methanogenic bacteria are very sensitive to pH decreases (Vanrolleghem and Lee, 2003). Hence, a reduction of the biogas production is observed. Obviously, the output alkalinity decreases as well.
Increasing the input pH results in an increased formation of carbonate precipitates in the digester, whereas decreasing the pH stimulates the stripping of CO 2 (see carbonate equilibria as function of pH; Zumdahl, 2005). Setting the pH inhibition level of sulfate reducing bacteria (SRBs) at 5, but for the other bacteria at 6, leads to increased H 2 S production if the pH in the digester becomes lower than 6.
Hence, the other bacteria are inhibited, whereas the SRBs still work at pH values lower than 6. Increasing the temperature in the digester stimulates the production of biogas. The increased temperatures facilitate faster reaction rates, and thus more biogas can be produced from the organic matter in an equal amount of time (Tchobanoglous et al., 2003).

NRM-Prec
Decreasing the P concentration in the input waste flow reduces the potential for struvite (MgNH 4 PO 4 :6H 2 O) precipitation. Decreasing the Mg concentration in the input waste flow decreases the pH in the reactor, which is obvious as a Mg source is  often added to induce P precipitation (Le Corre et al., 2007). Hence, less Mg-P precipitates are formed, the effluent P concentration increases, while the P recovery efficiency decreases. Increasing the P concentration in the input waste flow at a particular (neutral to high) pH increases the amount of P precipitates formed (precipitation is driven by supersaturation). Decreasing the pH by decreasing the concentration of nutrients, such as Mg and Ca, in the input waste flow reduces the resulting fertilizer density and molecular weight (fewer and less heavy P precipitates).

NRM-Strip
Decreasing the reactor height has no influence on the N recovery efficiency because the NH 3 -NH 4 þ equilibrium between a gas bubble and the surrounding water is reached in a very small time interval (Gujer, 2008). Increasing the temperature increases the NH 3 stripping performance (Wang et al., 2007). Hence, lower effluent NH 4 -N concentrations and higher NH 3 partial pressures in the gas phase are found. The more NH 3 is stripped out, the lower the effluent pH.
Increasing the liquid flow rate, reduces the residence time in the system. As such, the (slow) formation of CaCO 3 precipitates in the reactor is reduced, and thus also the scaling potential.

NRM-Scrub
Decreasing the reactor height has no influence on the N recovery efficiency because the NH 3 -NH 4 þ equilibrium between a gas bubble and the surrounding water is reached in a very small time interval (Gujer, 2008). Increasing the partial pressure of NH 3 in the incoming gas phase (coming from the stripper) decreases the fertilizer alkalinity (through NH 2 COO À formation) and increases the N concentration in the resulting ammonium sulfate solution. Hence, more N can be recovered in an equal amount of time.