Population balance models (PBMs) represent a powerful modelling framework for the description of the dynamics of properties that are characterised by distributions. This distribution of properties under transient conditions has been demonstrated in many chemical engineering applications. Modelling efforts of several current and future unit processes in wastewater treatment plants could potentially benefit from this framework, especially when distributed dynamics have a significant impact on the overall unit process performance. In these cases, current models that rely on average properties cannot sufficiently capture the true behaviour and even lead to completely wrong conclusions. Examples of distributed properties are bubble size, floc size, crystal size or granule size. In these cases, PBMs can be used to develop new knowledge that can be embedded in our current models to improve their predictive capability. Hence, PBMs should be regarded as a complementary modelling framework to biokinetic models. This paper provides an overview of current applications, future potential and limitations of PBMs in the field of wastewater treatment modelling, thereby looking over the fence to other scientific disciplines.

## INTRODUCTION TO POPULATION BALANCE MODELS

*x*is termed the internal coordinate (i.e. the distributed property),

*f*

_{1}(

*x*,

*t*) is the number density function (i.e. the distribution of the distributed property), is the continuous growth term of

*x*and

*h*(

*x*,

*t*) is the PBM reaction term (through discrete events). Table 1 provides some examples in regard to WWT applications.

Process | X | f_{1}(x,t) | h(x,t) | |
---|---|---|---|---|

(De)flocculation^{a} | Floc size | Floc size distribution | Microbial growth | Aggregation, breakage |

Size/density^{b} | Size/density distribution | Aggregation | ||

Bubble dynamics | Bubble size | Bubble size distribution | Bubble growth due to pressure reduction | Coalescence, breakage |

Granulation | Granule size | Granule size distribution | Microbial growth | Granulation |

Crystallisation | Crystal size | Crystal size distribution | Crystal growth | Aggregation, breakup |

Bio P-removal | Poly-P | Poly-P distribution | Poly-P storage, release | Cell division, cell birth |

Polyhydroxybutyrate (PHB) production | PHB | PHB distribution | PHB storage, release | Cell division, cell birth |

Growth | μ_{max} | μ_{max} distribution | Growth rate gradient | Cell division, cell birth |

Affinity | Kx | Kx distribution | Affinity gradient | Cell division, cell birth |

Process | X | f_{1}(x,t) | h(x,t) | |
---|---|---|---|---|

(De)flocculation^{a} | Floc size | Floc size distribution | Microbial growth | Aggregation, breakage |

Size/density^{b} | Size/density distribution | Aggregation | ||

Bubble dynamics | Bubble size | Bubble size distribution | Bubble growth due to pressure reduction | Coalescence, breakage |

Granulation | Granule size | Granule size distribution | Microbial growth | Granulation |

Crystallisation | Crystal size | Crystal size distribution | Crystal growth | Aggregation, breakup |

Bio P-removal | Poly-P | Poly-P distribution | Poly-P storage, release | Cell division, cell birth |

Polyhydroxybutyrate (PHB) production | PHB | PHB distribution | PHB storage, release | Cell division, cell birth |

Growth | μ_{max} | μ_{max} distribution | Growth rate gradient | Cell division, cell birth |

Affinity | Kx | Kx distribution | Affinity gradient | Cell division, cell birth |

^{a}Note that this mechanism is driving the settling processes in primary and secondary sedimentation.

^{b}In this case a two dimensional (2D) PBM is obtained.

The first term on the left-hand side of Equation (1) represents the accumulation term. Distribution dynamics that can be described are either governed by continuous processes (e.g. biomass growth, crystal growth, particle drying – represented by the second term on the left-hand side of Equation (1)) or discrete processes (e.g. aggregation, breakage, coalescence and granulation – represented by the term on the right-hand side). The mathematical functions used to describe these processes are called kernels. The right-hand side term usually consists of a birth and a death term, where the birth rate describes the rate at which particles of property *x* are being formed and the death rate describes the rate at which they are being removed. In crystallisation, a nucleation term needs to be added in the smallest size class mimicking the birth of a crystal through precipitation from an oversaturated solution. The internal coordinate *x* can be either a scalar (i.e. a single independent variable) or a vector, resulting in a one-dimensional (1D) or multi-dimensional PBM, respectively. The use of multi-dimensional PBMs means that the distribution of more than one material property can be described. It should be noted that formulating (i.e. writing down kinetic expressions or kernels) and solving multi-dimensional PBMs still poses a challenge to the research community. The nature of the resulting equation depends on the presence of the and *h*(*x*,*t*) terms. If only the former is present, a partial differential equation (PDE) is obtained for which solution methods are available. *h*(*x*,*t*) usually contains integral terms expressing the interactions between members of the distribution, turning the equation into an integro-PDE (see, for example, Nopens *et al*. (2002)). Several numerical methods have been reported in the literature to solve this type of equation (Ramkrishna 2000). Applications of PBMs to WWT processes are rather scarce. The first application was introduced by Fukushi *et al*. (1995), where a PBM model was used for modelling the dissolved air flotation process in water and wastewater treatment. The authors successfully described and validated the attachment process of bubbles to flocs during the flocculation process in a turbulent flow, allowing better prediction of the flotation process. Gujer (2002) investigated the impact of lumped average cell composition versus distributed composition in the context of Activated Sludge Models 2 and 3 and concluded that adding up the behaviour of individual cells does not result in the same model prediction as when predictions are made with average cell composition. Schuler (2005) demonstrated that lumped state (= averaged) assumptions in enhanced biological phosphorus removal system performance models (e.g. storage compounds and poly-P) produced large errors due to the difference in individual residence times of organisms in different zones (aerobic, anoxic and anaerobic). This was found to be related to process hydraulics (Schuler 2006) and to affect the endogenous respiration as the latter was found to be more important when distributed models were applied (Schuler & Jassby 2007). Finally, several PBM references can be found in the field of activated sludge flocculation, ranging from very simple formulations (Parker *et al*. 1972) to more elaborate ones (Nopens *et al*. 2002; Biggs *et al*. 2003) and papers focussing on experimental validation (Nopens *et al*. 2005) and model structure analysis (Nopens *et al*. 2007; Torfs *et al*. 2012). This will be discussed further in the next section. PBMs can serve the purpose of building process understanding. The result of such a detailed modelling exercise can be included in next-generation simplified wastewater treatment plant (WWTP) models that go beyond the currently used paradigms (i.e. activated sludge model using averaged biomass behaviour combined with residence time distribution models and oversimplified aeration and settling models). Hence, PBM models should not be considered as replacement of current WWTP models, but as enhancement tools to improve the future quality of their unit process predictions. This paper intends to outline the potential of PBMs in the field of WWT through several examples of different WWTP unit processes.

## POTENTIAL APPLICATIONS OF PBMS IN THE FIELD OF WWT

### Improved flocculation to better exploit primary and secondary settling

*et al*. 2014). Moreover, the primary treatment process is often chemically enhanced, which creates an optimal dosage problem. As particle concentrations are low in primary settlers, the settling regime is not hindered but rather discrete, i.e. Stokian, and a function of particle size, shape and density, leading to a wide distribution of settling velocities as evidenced by Bachis

*et al*. (2014). The discrete settling assumption is also true for the clarification zone above the sludge blanket of a secondary settler. In discrete settling, settling velocities are directly related to size, shape and density and, hence, the particle size distribution (PSD). The PSD depends on the original flocculation state as well as on actions undertaken to improve the flocculation state (e.g. turbulent shear, coagulant addition). Note that the flocculation state also depends on the particle's history (e.g. a sludge can have the same PSD, but it can react very different if the floc strength is different, caused by a different flocculation history). A PBM can account for this either through its aggregation and breakage rates or by using a two-dimensional PBM that considers the combination of size and density as internal coordinates. Flocculation of particles is probably the most straightforward application of PBMs. Since biological growth occurs on a much longer time scale compared to aggregation and breakage it can be ignored when studying short-term flocculation behaviour of activated sludge. The PBM reaction term is then defined as in which the birth (

*B*) and death (

*D*) terms occur for both aggregation and breakage as shown in Figure 1. Through these different mechanisms flocs of any size can be formed or removed. Furthermore, it becomes clear that aggregation is a particle–particle interaction process, whereas shear-induced breakage is not. The rates of all these processes are, hence, governed by the number of flocs present (

*N*) as well as an aggregation rate (

*β*), an aggregation efficiency (

*α*) and a breakage rate (

*S*) and distribution of resulting particles (the so-called daughter size distribution). Aggregation and breakage rates are in their turn a function of the mechanisms that drive the aggregation or breakup. Traditionally, the mechanisms of shear and polymer addition are included, which is probably sufficient for the application of PBM in the context of primary and secondary settling. More details on these rates and their dependencies on shear and flocculant addition can be found elsewhere (Nopens

*et al*. 2002; Nopens

*et al*. 2005). An example of a model prediction along with measured size distribution during a batch sludge flocculation process is shown in Figure 2.

Flocculation usually takes place in the process units prior to the actual settling tank as well as in the settling tanks (if conditions are appropriate). Flocculation models as described above can be used for both. Currently, the flocculation process in primary and secondary clarifiers is not studied in detail and its effects are incorporated using rules of thumb. However, understanding the contribution of flocculation would improve their design and operation, which can significantly improve the settling performance and control (both with respect to underflow and effluent suspended solids, the latter being a source of chemical oxygen demand (COD) and nutrients and, as such, important for effluent quality). Indeed, being able to predict and control the size distribution of a population of particles arriving to either the primary or secondary settler would be a useful input to settler models that can handle a distribution of settling velocities, calculated from the size distribution derived with the PBM (Bachis *et al*. 2014). In a secondary settler, exposure of flocs to elevated shear during transport from the bioreactor to the centre well of the settler will induce reflocculation and impact the floc size distribution and the floc strength. An appropriate application of PBM for settler-induced flocculation would specify size as the internal coordinate. In primary settler applications, the availability of particle settling velocity distributions as the internal coordinate of the suspended solids has generated a more accurate prediction of the load to the secondary treatment model, reduced the need for calibration (Bachis *et al*. 2012, 2014), and produced more accurate and optimal control of chemical dosage that leads to cost savings. In addition to particle size, particle density can be included as a second internal coordinate when density varies significantly with floc size. This additional internal coordinate can be accomplished using a 2D PBM approach but comes with an increased computational and parameter estimation cost since the rate expressions need to be extended to include density, which will require a detailed investigation of the process. Another interesting route for further research is the coupling of PBMs to computational fluid dynamics (CFD) models as the latter can predict local shear which then serves as input for the PBM model. Research that couples PBM with CFD in WWT has been reported already (Griborio & McCorquodale 2006; Gong *et al*. 2011), but needs further attention. Here, again, it should be clear that coupling a 1D PBM to a CFD model is a challenging task, typically resulting in models which need very long simulation times. One strategy to reduce the computational burden is to reduce the PBM model before coupling it to a CFD model (Mortier *et al*. 2013).

### More accurate aeration modelling for better design and energy optimisation

For a long time, K_{L}a-based models have been used to capture mass transfer between the gas and adjacent phases during aeration. More recently, models taking air flow rate as input were proposed. Despite the inclusion of somewhat more complexity and the resulting improved model performance, the variability of the *α* ‘fudging factor’ in space could still not be entirely related to process variables such as sludge concentration and sludge age, i.e. a lot of unexplained variance remains. Moreover, *α* was shown to vary spatially in a reactor (Rosso *et al*. 2011). This spatial variation introduces a significant amount of uncertainty in the model prediction when a single *α* value is used, resulting in locally different dissolved oxygen concentrations and, hence, aerobic process rates. To date, a key assumption in all aeration models is assigning a single average bubble size. This assumption is not very apparent, but resides in the gas–liquid interface surface area (a) of K_{L}a and is hidden in the *α* value in oxygen transfer efficiency-based models. This constant bubble size assumption is unrealistic and can be very restrictive for the model since the bubble's interfacial area drives the oxygen transfer process. In reality, bubble size is spatially distributed (Figure 3) from the point of injection to the top of a reactor due to the process of coalescence leading to a significantly different bubble size distribution near the reactor top. Another factor that plays a role here is the pressure reduction when moving from the bottom to the top of the reactor, which will influence bubble size. Increased viscosity (due to the presence of sludge) further promotes coalescence compared to clean water (Figure 4) (Fabiyi & Novak 2008; Ratkovich *et al*. 2013). A PBM using bubble size as the internal coordinate and including coalescence as a PBM reaction process can significantly improve the local prediction of oxygen mass transfer (and hence *α*) as well as improve the design of aeration systems to maximise the oxygen transfer (in combination with CFD). It should be noted that the current work in CFD linked to aeration also uses a fixed bubble size (Fayolle *et al*. 2007). The use of PBMs in aeration systems with suspended solids has not been widely studied in WWT. PBMs applied to bubble columns are widespread in fermentation systems and the chemical engineering literature (Dhanasekharan *et al*. 2005; Sanyal *et al*. 2005; Wang 2011) and can serve as solid examples of improved benefits to these process models.

Typical mechanisms taking place in bubble breakup and coalescence are shown in Figure 5. The kernels used in a PBM describing bubble breakup and coalescence are very similar as those used in a flocculation PBM.

### Floc and granule size: distributed kinetics

Despite the fact that the size of agglomerates, be it flocs or granules, can play an important role in biological activity (e.g. different microbial consortia at different locations in the agglomerate), the simulation of such a microenvironment within a macro-scale fluid transport environment has hardly been performed. Recently, Volcke *et al*. (2012) demonstrated the significant impact of granule size distribution on the performance of an Anammox-based granular sludge reactor. The authors used a fixed size distribution for this analysis. It is clear that size distribution can be impacted by shear and, hence, will further influence the system behaviour. The absence of size and size dynamics in the currently used models indicates that the activity loss caused by particle size (causing transport limitation due to stratification, e.g. Vangsgaard *et al*. (2012)) cannot be predicted. Again, when experimental data are confronted with these models, other degrees of freedom (i.e. parameters or input variables) will be calibrated for inappropriate reasons. Sobremisana *et al*. (2011) demonstrated that including floc size can result in significantly deviating reactor performance since kinetics can be quite different depending on reactor location and the size of the biological floc. The authors used an integrated PBM–CFD approach to simulate the carbon and nitrogen removal process at both the reactor scale and internal floc scale. The effect of size was introduced by means of an effectiveness factor (i.e. ratio of rate with and without diffusional resistance) based on floc size for all different processes. For a simple baffled reactor the treatment performance deviated by 13% for COD, 10% for NH_{4} and 56% for NO_{3} compared to the same simulation not accounting for influence of size. However, further validation is required. Knowledge arising from this can be useful to (partially) decouple affinity constants in kinetic rate expressions and reduce their requirement for calibration. Understanding the interaction between size and reaction kinetics can inform researchers and engineers on how to better design and operate these processes (e.g. avoid or promote certain shear zones and account for imposed shear of mixing and aeration intensities). Moreover, it will reduce the need to adjust parameters unnecessarily to improve the model fit. Apart from size heterogeneity, incomplete mixing can lead to spatially heterogeneous concentrations in biomass and substrates that ultimately result in locally different kinetics. Integrating the effects of spatial variations in macro-scale mixing as well as the biomass and substrate concentration is another avenue for further model development. Lencastre Fernandes *et al*. (2013) demonstrated the effect of this heterogeneity for a budding yeast population using a multi-scale modelling approach that included PBM. A similar approach could be used for WWT modelling, and could be helpful in developing an improved understanding of biomass population dynamics. A complicating factor for the WWT compared to a yeast fermentation is that multiple species have to be considered to model the WWT appropriately. Models including floc size and local heterogeneities could also be helpful in developing technologies to select between wanted and unwanted microbial communities, which is a recently developing topic in view of mainstream Anammox application. Microbial selection can be done on a physical basis (size) or through selection by creating favourable growth conditions for the targeted microbial consortium (Al-Omari *et al*. 2014).

### Precipitation/crystallisation for better quality marketable products

WWTPs are transforming into water resource recovery facilities (WRRFs) leading to new modelling challenges (Vanrolleghem 2013). One important aspect in product recovery, driving their market value, will be the specifications of the recovered material. These comprise both composition and size. Crystallisation has been extensively modelled in the field of chemical engineering and pharmaceutical engineering to produce crystals with tailor-made specifications (e.g. Aamir *et al*. 2009; Nagy & Braatz 2012). A PBM with crystal size as an internal coordinate and inclusion of nucleation (function of super-saturation) and crystal growth can be used as a first approximation. If needed, more internal coordinates can be added to deal with composition or crystal shape (e.g. 2D compared to 1D; Samad *et al*. (2011)). Additional phenomena such as agglomeration and breakage can be added. Interestingly, describing crystallisation with a PBM also allows the describing of phenomena such as size-dependent crystal growth (Samad *et al*. 2011). The use of a PBM modelling framework is widely accepted when studying crystallisation processes. However, in a WRRF context, the PBM framework has not really been used thus far, with the exception of a recent manuscript by Galbraith & Schneider (2014) where a discretised PBM was used to describe the chemical precipitation of phosphorus. The main phenomena are nucleation, crystal growth and dissolution (= negative growth), agglomeration and breakage. However, it happens frequently that the PBM only considers growth and nucleation (Fujiwara *et al*. 2005), which is related to the fact that it is difficult to calibrate model parameters when additional phenomena such as agglomeration and breakage are included as well. An important variable that needs to be included in these models is the super-saturation (the driving force for crystallisation, i.e. a state of a solution that contains more of the dissolved material than could be dissolved by the solvent under normal circumstances), which will vary as a function of temperature. Description of the latter temperature dependence requires knowledge about the super-saturation curve as a function of temperature, which is usually represented by means of a polynomial.

## WHEN TO USE A PBM AND REMAINING CHALLENGES

PBMs should be regarded as an additional modelling framework that can be used complementary to other frameworks (i.e. biokinetic models, CFD). However, their inclusion into models would only increase the overall complexity and consequently their use should only be justified if pertinent to the modelling objective. A primary use of PBM will therefore be as an aid to increase process understanding. Second, the required knowledge can be translated and added to currently used models. Hence, it is obvious that in the short-term PBMs will not become the mainstream modelling paradigm in WWT modelling, but rather an addition (similar to CFD). Although PBMs have been around for 50 years, the real growth in their application came with the increased computational power in the 1990s. Therefore, PBMs still need to further develop and several challenges remain: which phenomena should be included; which kinetic expressions to select (or develop) for each phenomenon; what is the appropriate number of classes? PBM validation due to cumbersome data collection is another challenge. Despite these outstanding challenges, the time has come to start reflecting and applying PBMs in this field in order to advance their development and make them a suitable complementary modelling framework in the short term to further improve our process understanding and our strive for improved technologies.

## CONCLUSIONS

Many processes in WWTPs are governed by population dynamics of materials characterised by variation in property dynamics. These potential complexities in system behaviour are lost or significantly suppressed when only average behaviour is characterised or simulated. PBMs can deal with these process complexities and have already demonstrated their benefits in the field of (bio)chemical engineering. The majority of the models in WWT modelling that need more rigour are physical–chemical processes. Hence, more than ever we need to look over the fence and integrate available (bio)chemical engineering knowledge into WWTP models. Some examples are described in this paper, but potentially many more applications of PBM in WWT exist and can be exploited. The intention of this paper is to make WWT modellers aware of this complementary modelling framework and its potential applications, challenges and limitations.

## ACKNOWLEDGEMENTS

P.V. holds the Canada Research Chair on Water Quality Modelling. We kindly acknowledge Dr Marina Arnaldos for her critical reflections on the paper. She currently is an experienced researcher in the FP7 EU ITN-project SANITAS. This paper was presented at WWTmod2014 and the fruitful discussions are kindly acknowledged.