Multiscale modelling of the start-up process of anammox-based granular reactors

1.   Introduction

In the last decades, granular sludge technologies have completely revolutionized the treatment and valorization of industrial and municipal wastewater as they can be applied for the simultaneous removal of organic, nitrogen and phosphorus compounds and the production of bioenergy ^[1]. Granular sludge reactors are biofilm systems where biomass grows arranged in granules, dense and compact aggregates with an approximately spherical shape ^[2]. In contrast to the traditional biofilm systems, where biofilms develop on solid surfaces, biofilm formation in granular sludge reactors occurs due to the self-immobilization of cells without the involvement of a surface, in a process known as granulation. Such process can be initiated from an inoculum in granular or suspended form. In the latter case, the process is named de novo granulation. Conversely to conventional wastewater systems, high amounts of extracellular polymeric substances (EPS) constitute biofilm granules which, consequently, have higher biomass densities, more regular shapes and stronger structures ^[3]. These characteristics also provide protection for sensitive microbial species which have difficulties to develop in suspended form ^[4]. Moreover, due to the high settling velocity of granular sludge, solid-liquid separation is facilitated ^[5] and high biomass concentrations can be achieved in the system ^[6]. Additionally, the geometry and free movement of granules limit external boundary layer resistances and promote the mass transport of substrates towards the various granule microbial layers ^[1]. All these features contribute to high removal efficiencies and reduced-footprint systems, and, consequently, specific granular biofilms, such as aerobic, anaerobic and anammox granules, have been adapted for various wastewater treatment processes.

Although the traditional process of nitrification/denitrification (N/D) is commonly used to remove nitrogen compounds from wastewater, it is energy-intensive as it requires dissolved oxygen supply for ammonium oxidation. Moreover, an external carbon source is necessary for the heterotrophic metabolism of denitrifying bacteria, in the case of low carbon-to-nitrogen ratio wastewater. Therefore, in recent years the combination of partial nitritation and anammox processes (PN/A) has been increasingly studied for the treatment of nitrogenous wastewater. Such process allows the conversion of ammonium into molecular nitrogen: in the first step the partial nitritation takes place and about half of the ammonium present in the wastewater is converted to nitrite by aerobic ammonia-oxidizing bacteria ($AOB$); then, in the subsequent anammox step, the nitrite produced and the remaining ammonium are simultaneously converted into nitrogen gas and small amounts of nitrate by anammox bacteria ($AMX$). In the case of suspended biomass, the two processes occur in separate reactors arranged in series as they require aerobic and anoxic conditions. The granular biofilm technology represents a cost-effective and alternative solution, since both processes can be carried out simultaneously in a granular sludge reactor. Indeed, the formation of two distinct zones inside the granules is induced by providing a constant and appropriately low oxygen level in the reactor: an external zone where oxygen necessary for partial nitritation is guaranteed and an internal zone where oxygen is not present and ideal anoxic conditions for anammox processes occur.

Compared to the traditional N/D process, PN/A granular process results in lower aeration costs, $CO_2$ emissions and sludge production. Additionally, due to the autotrophic metabolism of anammox bacteria, no external carbon addition is required. For these reasons, PN/A granular process is considered a promising technology for N-removal. Nevertheless, anammox bacteria are very sensitive to environmental and operating conditions, and are characterized by very low growth rates and cellular yields ^[7,8,10,12]. Consequently, the start-up of anammox granular sludge reactors is a long and complex process, which represents the main drawback of this technology ^[8,11] and deeply impacts the operating strategies and procedures ^[13].

The start-up process of anammox granular sludge systems can be divided into four phases: cell lysis phase, lag phase, activity elevation phase and stationary phase ^[8,9,10,15]. The lysis phase occurs frequently when microbial species constituting the granular sludge inoculum find new and unknown environmental conditions and carry out processes of microbial autolysis leading to the disintegration of biofilm granules. A transition period occurs in the successive lag phase, when the biomass begins to adapt to the reactor conditions and new biofilm granules begin to develop. In the elevation phase, the well-adapted biomass grows and a continuous and increasing ammonium and nitrite removal is observed. Finally, in the stationary phase, red mature granules are detected, dominated by anammox bacteria, with the presence of ammonia-oxidizing bacteria and denitrifying heterotrophic bacteria. In this phase, an optimum and stable N-removal efficiency is achieved ^[8,14].

Many studies deal with the start-up of PN/A granular sludge systems by focusing on factors which govern the biological processes involved and the growth of anammox bacteria, such as inoculum sludge, hydraulic retention time, dissolved oxygen, temperature, pH, wastewater composition and reactor configuration. Among these, the inoculation procedure appears to be a key element to reduce the start-up period ^[14,16]. In experimental works reported in literature, different selected sludges have been used as inocula of lab-, pilot- and full-scale bioreactors, such as anaerobic granular sludge ^[7,8,9,10], flocculant nitrification sludge ^[7,9,13], flocculant denitrification sludge ^[7], activated sludge with or without addition of anammox sludge ^[15]. Anaerobic granular sludge can be a competitive alternative solution, as the biomass is already in granular form and acclimatized to anaerobic conditions. Moreover, the conventional flocculant sludge represents a further alternative, due to its greater availability. The latter two types of inoculum can also be enriched with anammox sludge. In this context, the search for new procedures and strategies that can significantly accelerate the process start-up remains an interesting challenge.

The study of the anammox process and start-up strategies through experimental activities requires high cost and long time, mainly due to very low growth rates of anammox bacteria. Consequently, mathematical modelling appears to be an attractive alternative solution for the description and optimization of PN/A granular sludge reactors, for the understanding and investigation of microbial dynamics which govern the growth of anammox granules, and for testing a wide range of environmental and operational conditions which could influence the process. In the recent years, numerous mathematical models have been proposed to describe anaerobic ^{[17,18,19,20]}, aerobic ^[21,22,23] and anammox processes ^{[24,25,26,27,28]} in granular-based systems, by considering the evolution of biofilm granules. Depending on the approach used to model the development and the structure of granular biofilms, two types of models can be distinguished: continuum and discrete models. The first ones describe granular biofilms as spherical continuum domains, through a quantitative and deterministic approach ^{[18,19,21,25]}, while discrete models, such as individual-based models ^[23,29], consider microbial cells as discrete entities and introduce elements of randomness and stochastic effects in the solution. Most of continuous models have been formulated as spherical free boundary problems with radial symmetry. Some of them take into account attachment ^[17,30,31] and invasion ^[31] processes. Moreover, the initial formation of biofilm granules has been modelled in Tenore et al. ^[31] by setting a zero initial granule radius. Continuum models frequently assume one single granule size class ^[25,27,28], while someone takes into account the size distribution of granules within the reactor by considering more size classes ^[18,19,26]. However, some works ^[19,26] demonstrate that one single size class allows to correctly describe the global treatment process, while the granule size distribution could be required to investigate more specific aspects, such as the microbial composition and the solute exchange between granules of different sizes. Although most models consider all the biomass in the granular form, a few include both the sessile biomass and the planktonic biomass present in the reactor and model the microbial mass fluxes between granular biofilms and liquid medium ^[28,31]. Many free boundary models describe the granules evolution by fixing the steady-state dimension ^{[25,26,27,28]}, while in other works the steady-state dimension is supposed to be a function of microbial metabolic activities and operating conditions of the system ^[19,32]. Except for a few cases ^[33], almost all models on granular biofilms consider the perfect retention of biofilm granules and, consequently, assume the number of granules in the system as a constant ^[21,25].

Among the models on granular biofilms, some focus on processes of partial nitritation and anammox ^[25,26,27], by taking into account the dynamics of main microbial species involved, such as aerobic ammonia-oxidizing bacteria, nitrite-oxidizing bacteria ($NOB$), anammox bacteria and main soluble substrates, such as ammonium, nitrite and nitrate. However, someone includes also the conversion of organic compounds by heterotrophic bacteria ($HB$) ^[34]. Indeed, when sufficient amounts of organic carbon occur in the system, heterotrophic bacteria can play a negative role in the N-removal process, by proliferating within the granules and competing with $AOB$ and $AMX$.

Models on PN/A granular systems propose some interesting numerical studies, aimed at describing the PN/A process and looking into the effects of some key factors on the reactor performance and microbial composition of granules. Specifically, the effects of the following factors have been investigated: bulk oxygen concentration ^[24,25] and aeration pattern ^[35], granule size ^{[24,25,26,34]} and granule size distribution ^[26], influent concentration ammonium ^[24,25], influent concentration COD and heterotrophic growth ^[28,34]. Furthermore, the impact of the coexistence of microbial flocs and granular biofilm on the reactor performance has been studied by Hubaux et al. ^[28]. The emission of nitrous oxide and nitric oxide occurring in the treatment process has been investigated by Vangsgaard et al. ^[27], exploring the effects of ammonium load, granule size and temperature. Finally, a model which propose the integration of methane removal in PN/A granular sludge reactors is reported by Castro-Barros et al. ^[36]. Nevertheless, all models related to anammox granular biofilms completely neglect the attachment process in the initial phase of biofilm formation, and the initial data that prescribe biofilm domain size and species composition are arbitrarily assigned.

In this context, the present work proposes a mathematical model aimed at describing and investigating aspects of the PN/A granular sludge bioreactor that have never been addressed in literature, such as the initial formation of anammox granules (anammox de novo granulation) and the start-up process. The general framework has been introduced in Tenore et al. ^[31] for anaerobic biofilm granules and has been applied here to PN/A processes. The model is formulated as a spherical free boundary value problem under the assumption of radial symmetry. The de novo granulation process is modelled by assuming a vanishing initial granule size ^[37,38]. This means that all the biomass present in the system at the beginning of the process is supposed in planktonic form (flocculant sludge inoculum). Then, the granules formation is initiated by the attachment phenomena. This first stage in biofilm formation process is followed by a subsequent colonization of anammox bacteria that benefit from the protective environment characterizing the anoxic core of granular biofilms. By using a continuum approach ^[39,40], the model takes into account the dynamics of soluble substrates and biomasses in planktonic and sessile form. In particular, processes of microbial growth, substrates conversion, microbial invasion, attachment and detachment are included in the model.

This model has been integrated numerically by developing an original code in the MATLAB platform and numerical studies have been carried out for the following purposes: (i) test the model behaviour; (ii) examine the fundamental role of anammox bacteria invasion in the de novo granulation process; (iii) explore the formation, evolution and ecology of anammox granules in PN/A granular sludge systems; (iv) study the autotrophic nitrogen removal through PN/A processes; (v) investigate and optimize the process start-up of these bioreactors. In particular, numerical simulations have been carried out to investigate how the size and the addition time of the anammox inoculum can affect the invasion process and optimize the process start-up. The numerical results refer to both the individual biofilm granule and the macroscale reactor performance and include the distribution and relative abundance of active sessile biomasses within the granule, the evolution of granule dimension, and the profiles of soluble substrates and planktonic biomasses within the reactor. This work is organized as follows. The mathematical model is reported in Section 2, while the biological context is described in Section 3, where model variables and kinetic rate equations are introduced. Numerical studies are presented and discussed in Section 4. Finally, the conclusions are outlined in Section 5.

2.   Mathematical model

As mentioned in Section 1, the mathematical formulation of granular biofilm reactors presented in Tenore et al. ^[31] for anaerobic digestion is applied here to partial nitritation/anammox granular processes. In this Section, the model equations are recalled for convenience.

A granular-based reactor is an extremely complex multiphase biological system which requires the introduction of some assumptions to be modelled. Such a system is supposed here to consist of two distinct components: the granular biofilm phase and the bulk liquid phase. These components influence each other through continuous mass exchanges involving sessile and planktonic biomass and soluble substrates. The granular biofilm phase is represented by a fixed number of biofilm granules ($N_G$) immersed within the bulk liquid and assumed to have identical properties at any instant of time. Specifically, biofilm granules are modelled as spherical free boundary domains with radial symmetry and with a vanishing initial radius. The attachment flux of planktonic biomass from the bulk liquid is accounted to initiate the granulation process. Biofilm granules evolves over time as a result of various processes such as metabolic activities, detachment and transport of sessile biomass, diffusion and conversion of soluble substrates and invasion of planktonic species. The term attachment is used here to indicate the aggregation of planktonic cells, which contributes to the genesis and growth of biofilm granules. Invasion phenomena consist in the colonization of a pre-existing granule mediated by planktonic motile cells living in the surrounding environment, which can penetrate the porous matrix of the biofilm and convert to sessile biomass. Detachment phenomena lead to sessile biomass losses, induced by external shear forces, substrates depletion and biomass decay.

In order to model all these processes, the following model variables have been considered within the granular biofilm domain:

● radius of the biofilm granule: $R(t)$,

● concentration of $n$ sessile species: $X_i(r, t), \ i = 1, ..., n$,

● concentration of $n$ planktonic invading species: $\psi_i(r, t),$ $\ i = 1, ..., n$,

● concentration of $m$ dissolved substrates: $S_j(r, t), \ j = 1, ..., m$.

The free boundary domain is described by the time-dependent granule radius R(t), while all other variables are expressed as functions of time $t$ and radial coordinate $r$, where $r = 0$ identifies the granule center. The liquid present in the voids of granules is not included as a model variable, since it is supposed to not play a limiting role in the microbial metabolic activities. Planktonic cells are supposed to have a negligible spatial occupancy due to the small particle size. Assuming that all sessile species are characterized by the same constant density $\rho$, the biofilm volume fraction of each individual species $f_i$ can be calculated by dividing $X_i$ by $\rho$. Furthermore, it is assumed that the sum of the biomass volume fractions is equal to one at each location and time, $\sum_{i = 1}^{n} f_i = 1$ ^[41]. Since $X_i$ and $f_i$ are mutually dependent variables, only $f_i$ is included among the model unknowns.

The reactor is modelled as a completely mixed continuous system. Therefore, the properties of the bulk liquid are the same at every point and change over time due to the conversion processes of the planktonic biomass and soluble substrates present in the bulk liquid and due to mass exchanges with the biofilm granules. In order to take into account these aspects, the following model variables have been considered within the bulk liquid:

● concentration of $n$ planktonic biomasses: $\psi^*_i(t), \ i = 1, ..., n$,

● concentration of $m$ dissolved substrates: $S^*_j(t), \ j = 1, ..., m$.

In the following lines, all model equations and boundary and initial conditions related to the biofilm domain and the bulk liquid domain are reported.

The growth and the transport of the $i^{th}$ sessile species across the granular biofilm is governed by the following hyperbolic partial differential equations (PDEs):

where $u(r, t)$ is the biomass velocity, $r_{M, i}(r, t, \bf X, \bf S)$ and $r_i(r, t, {\boldsymbol{\psi }}, \bf S)$ are the specific growth rates due to sessile and planktonic species, respectively, ${\bf X} = (X_1, ..., X_n)$, ${\bf S} = (S_1, ..., S_m)$, ${\boldsymbol{\psi }} = (\psi_1, ..., \psi_n)$, $v_{a, i}$ is the attachment velocity of the $i^{th}$ planktonic biomass and $\psi^*_i(t)$ is the concentration of the $i^{th}$ planktonic biomass in the bulk liquid.

When the attachment flux from bulk liquid to granule $\sigma_{a}(t)$ is higher than detachment flux from granule to bulk liquid $\sigma_{d}(t)$, the free boundary is a space-like line and the condition (2.2) at the interface granule-bulk liquid $r = R(t)$ is required. Conversely, when $\sigma_{a}(t)$ is lower than $\sigma_{d}(t)$, the free boundary is a time-like line and the condition (2.2) is not needed.

The function $u(r, t)$ satisfies the following problem:

where $G(r, t, \textbf{f}, \textbf{S}, {\boldsymbol{\psi }}) = \sum_{i = 1}^{n}(r_{M, i}(r, t, \textbf{f}, \textbf{S})+r_i(r, t, {\boldsymbol{\psi }}, \textbf{S}))$ and ${\bf f} = (f_1, ..., f_n)$.

By considering Eq. (2.3), Eqs. (2.1) and (2.2) can be rewritten as follows:

The free boundary evolution is described by the granule radius $R(t)$ and depends on processes of sessile metabolic growth, detachment and attachment. Attachment phenomena dominate the granulation process, while detachment phenomena become predominant as the granule dimension increases. The variation of $R(t)$ is governed by the following ordinary differential equation (ODE), derived from the mass balance on the granule volume ^[31]:

In particular, attachment process is modelled through a continuous flux from bulk liquid to granule, given by the sum of the attachment fluxes of each planktonic species $\sigma_{a, i}(t)$, which are linearly dependent on the concentration of planktonic biomasses within the bulk liquid $\psi^*_i(t)$ ^[37,38]:

Meanwhile, the detachment process is modelled through a continuous flux from granule to bulk liquid, which is a quadratic function of the granule radius $R(t)$ ^[42]:

where $\lambda$ is the detachment coefficient and is supposed to be equal for all microbial species.

The diffusion and conversion of planktonic cells and soluble substrates within the biofilm granule are governed by the following parabolic PDEs:

where $r_{\psi, i}(r, t, {\boldsymbol{\psi }}, {\bf S})$ represents the loss rate for the invading species; $r_{S, j}(r, t, {\bf f}, {\bf S})$ represents the substrate production or consumption rate due to microbial metabolism; $D_{\psi, i}$ and $D_{S, j}$ denote the diffusion coefficients within the biofilm for the $i^{th}$ planktonic species and $j^{th}$ dissolved substrate, $\psi^*_i(t)$ and $S^*_j(t)$ denote the concentrations of planktonic cells and dissolved substrates within the bulk liquid, respectively. All equations which refer to the biofilm domain do not require initial conditions, since the extension of the biofilm domain is set zero at $t = 0$.

$\psi^*_i(t)$ and $S^*_j(t)$ represent the solutions of the following ordinary differential equations, which describe the dynamics of planktonic biomass and soluble substrates within the bulk liquid, respectively, and are derived from mass balances on the bulk liquid volume:

where $V$ is the volume of the bulk liquid, assumed equal to the reactor volume, $Q$ is the continuous flow rate, $A(t)$ is the area of the granule and is equal to $4 \pi R^2(t)$, $\psi^{in}_i$ is the concentration of the planktonic species $i$ in the influent, $S^{in}_j$ is the concentration of the substrate $j$ in the influent, $r^*_{\psi, i}(r, t, {\boldsymbol{\psi }}^*, {\bf S^*})$ and $r^*_{S, j}(r, t, {\boldsymbol{\psi }}^*, {\bf S^*})$ are the conversion rates for $\psi^*_i$ and $S^*_i$, $\psi^*_{i, 0}$ is the initial concentrations of the $i^{th}$ planktonic species within the bulk liquid, ${\bf S^*} = (S^*_1, ..., S^*_m)$, ${\boldsymbol{\psi }}^* = (\psi^*_1, ..., \psi^*_n)$.

No contribution by detachment to planktonic or detached biomass is considered in this model. Indeed, the detached biomass has different characteristics from both sessile and planktonic biomass, and several hours are required for its conversion into the planktonic form ^[43,44,45]. Moreover, granular-based reactors are typically characterized by high selection pressures (low $HRT$ and high velocities) to promote the biomass aggregation ^[2]. Under these high selection pressures, granules are retained in the reactor, while planktonic and detached cells are rapidly washed out ^[4].

The mass of the $i^{th}$ sessile species within the granule can be derived from:

while, the total mass can be calculated as follow:

3.   Modelling de novo anammox granulation

The mathematical model described in the previous Section has been applied to the partial nitritation/anammox process, with the aim of describing the dynamics of anammox granules and investigating the start-up of combined partial nitritation-anammox reactors.

The anaerobic ammonia oxidation process (anammox process) allows to remove nitrogen from wastewater via anaerobic pathways of specific autotrophic microbial species, known as anammox bacteria. Such bacteria use ammonium as electron donor to convert nitrite into nitrogen gas and small fractions of nitrate. However, since nitrite are not commonly present in nitrogenous wastewater, the anammox process is preceded by a partial nitritation step, where the necessary amount of nitrite is produced. As mentioned in Section 1, under appropriate operating conditions these two processes can be carried out simultaneously in one single granular sludge reactor, exploiting the coexistence of anoxic and aerobic zones within biofilm granules.

In order to comprehensively model the treatment process and the evolution of granules occurring in these reactors, all the main biological processes, microbial species and soluble substrates have been considered. Specifically, processes of nitritation, anammox, denitrification, organic carbon and nitrite oxidation are supposed to occur in biofilm granules and in the bulk liquid, induced by the metabolic activities of the planktonic and sessile biomass. Hence, the following active microbial species have been considered both in sessile and planktonic form: aerobic ammonia-oxidizing bacteria $AOB$, anaerobic ammonia-oxidizing (anammox) bacteria $AMX$, aerobic nitrite-oxidizing bacteria $NOB$ and facultative heterotrophic bacteria $HB$. All sessile species are supposed to decay and produce sessile inactive biomass $I$ which accumulate within the biofilm granules. Conversely, although it is assumed that planktonic biomasses also decay, the inactive biomass in planktonic form has not been included in the model because it is likely to play a negligible role in the development of the successive processes. Moreover, in order to describe the metabolic activity of the active microbial species, the following soluble compounds have been modelled: ammonium $NH_4$, nitrite $NO_2$, nitrate $NO_3$, soluble organic carbon $OC$ and oxygen $O_2$.

During the nitritation process, ammonium-oxidizing bacteria $AOB$ convert ammonium $NH_4$ and oxygen $O_2$ to nitrite $NO_2$ under aerobic conditions, according to the following reaction:

Under aerobic conditions, nitrite $NO_2$ is oxidized with $O_2$ to form nitrate $NO_3$ by nitrite-oxidizing bacteria $NOB$:

During anammox processes, under anoxic conditions, anammox bacteria $AMX$ convert ammonium $NH_4$ and nitrite $NO_2$ in nitrogen gas and little amounts of nitrate $NO_3$, as follows:

The model considers metabolic processes of facultative heterotrophic bacteria $HB$ as well. Specifically, $HB$ are supposed to growth both under aerobic and anoxic conditions. Under aerobic condition, they oxidate organic carbon $OC$ by using free oxygen $O_2$, while they carry out two denitrification reactions under anoxic conditions: first they oxidate $OC$ by using the nitrate bound oxygen and $NO_3$ reduces to nitrite $NO_2$; then they oxidize $OC$ by using nitrite bound oxygen and $NO_2$ reduces to molecular nitrogen.

In summary, the list of all model variables is reported below:

As mentioned in Section 1, a constant and appropriate oxygen level is maintained in PN/A granular reactors, in order to guarantee distinct zones for ideal growth of both anaerobic and aerobic species. To model this, oxygen concentration is assumed to be variable only within the biofilm $S_{O_2}(r, t)$, where it varies due to microbial consumption and diffusion phenomena. Instead, oxygen concentration in the bulk liquid does not represent a model variable and is fixed at a constant value $S^*_{O_2}(t) = {\overline S_{O_2}}$.

Biological pathways described above have been included in the model through the mathematical formulation of the reaction terms. According to D'Acunto et al. ^[46], specific growth rates due to sessile species $r_{M, i}$ in Eqs. (2.3) and (2.5), with $i \in \{AOB, AMX, NOB, HB\}$, are modelled as Monod-type kinetics, while the formation rate of inactive biomass $r_{M, i}$ in Eqs. (2.3) and (2.5), with $i \in \{I\}$, is modelled as first order kinetic and given by the sum of decay rates of each active species. The rates $r_i$ in Eqs. (2.3) and (2.5), with $i \in \{AOB, AMX, NOB, HB\}$, are modelled as Monod-type kinetics as well.

Moreover, the rates $r_{\psi, i}$ in Eq. (2.11), with $i \in \{AOB, AMX, NOB,$ $HB\}$, and $r_{S, j}$ in Eq.(2.13), with $j \in \{NH_4, NO_2, NO_3, OC,$ $O_2\}$, are formulated from the corresponding microbial growth rates through the stoichiometric coefficients and specific microbial yield $Y_{\psi, i}$ and $Y_i$, respectively.

Similarly, the rates $r^*_{\psi, i}$ in Eq. (2.15), with $i \in \{AOB,$ $AMX, NOB, HB\}$, and $r^*_{S, j}$ in Eq.(2.17), with $j \in \{NH_4,$ $NO_2, NO_3, OC\}$, are expressed as Monod-type kinetics and are assumed to be proportional with each other through the stoichiometric coefficients and specific microbial yield $Y_i$.

All kinetic rate equations have been reported in Supplementary material, while the values used for all stoichiometric and kinetic parameters have been summarized in Supplementary Table S3.

4.   Numerical studies and results

4.1.   Influent characteristics, reactor configuration and simulation parameters

Numerical simulations have been carried out to test the model behaviour, simulate the evolution and ecology of anammox granular biofilms and investigate the treatment process occurring in PN/A granular sludge reactors, with a focus on the start-up phase.

The modelled influent wastewater represents a typical high ammonium wastewater treated in PN/A granular sludge reactors. It is characterized by $300 \ gN \ m^{-3}$ of ammonium and $50 \ gCOD \ m^{-3}$ of soluble organic carbon, while nitrite and nitrate amounts are supposed to be negligible. Specifically, $S^{in}_{NO_2}$ and $S^{in}_{NO_3}$ are set to $0.0001 \ gN \ m^{-3}$ in order to avoid numerical errors arising from zero concentrations in the kinetic expressions. The constant oxygen level within the reactor is fixed at $0.75 \ gO_2 \ m^{-3}$. Microbial biomasses are assumed to be not present in the influent ($\psi^{in}_i = 0$).

The strategy of using two separate inocula is studied: at $t = 0$ the bioreactor is inoculated with the activated sludge coming from a conventional nitrification-denitrification reactor, where anammox biomass is not present. Once the process has started, an anammox inoculum is added at a fixed time instant $\bar{t}$. Then, the $AMX$ planktonic cells are supposed to colonize the granules through invasion phenomena and grow in sessile form in the innermost part, where anoxic conditions are guaranteed. The parameter $\bar{t}$ has been varied in the simulations to investigate the effect of granules dimension on anammox growth. The activated sludge inoculum is modelled by setting the initial concentration of planktonic biomasses in the bulk liquid: $\psi^*_{AOB, 0} = \psi^*_{NOB, 0} = \psi^*_{HB, 0} = 300 \ gCOD \ m^{-3}$, $\psi^*_{AMX, 0} = 0$. Meanwhile, in order to consider the addition of the anammox inoculum at $\bar{t}$, Eq. (2.15) for $\psi^*_{AMX}$ has been replaced by the following impulsive ordinary differential equation (IDE):

where $\psi^*_{AMX, \bar{t}}$ is the concentration of anammox planktonic biomass added in the bulk liquid at $\bar{t}$ and is related to the anammox inoculum size. $\psi^*_{AMX}(\bar{t}^+)$ and $\psi^*_{AMX}(\bar{t}^-)$ are the right and left limits of $\psi^*_{AMX}$ at time $\bar{t}$. Since the parameters $\psi^*_{AMX, \bar{t}}$ and $\bar{t}$ are varied in numerical studies, their values are provided below, case to case.

Reactor volume $V$ is assumed equal to $400 \ m^3$ ^[25,26,28] and fed with a constant flow rate $Q$ of $2000 \ m^3 \ d^{-1}$ (hydraulic retention time $HRT = 0.2 \ d$). The number of granules $N_G$ has been selected through an iterative procedure which involved the detachment coefficient $\lambda$ ^[31], with the aim to guarantee a $25$% filling ratio ^[25,26,28] by considering $1 \ mm$ as steady-state particle radius (an average size representative of the anammox granules ^[25,26,28]). In accordance with ^[39], diffusivity of soluble substrates in biofilm is assumed to be $80\%$ of diffusivity in water. The values reported above for all operating parameters are characteristic of PN/A granular reactors ^[47]. All model parameters have been summarized in Supplementary Table S3 and Table 1.

Four numerical studies have been performed with the aim of investigating how the invasion phenomena of anammox bacteria affect the formation process of anammox granules and the distribution of microbial species involved in the biological processes:

● The first study (S1) presents a reference case to test the model behaviour and explores the de novo granulation of anammox granules and the global treatment process of PN/A granular bioreactors;

● The second study (S2) examines the effect of the anammox addition time $\bar{t}$ on the start-up process and granules evolution;

● The third study (S3) investigates the influence of anammmox inoculum size $\psi^*_{AMX, \bar{t}}$ on the start-up process and granules evolution;

● The last study (S4) analyzes the combined effects of both the addition time and the size of the anammox inoculum on the start-up process.

4.2.   S1 - Partial nitritation/anammox process in granular-based reactors

In the first numerical study (S1) the de novo anammox granulation and the dynamics of solutes and planktonic biomasses within the PN/A granular system are investigated. As mentioned in Section 1, the granular sludge reactor is initially inoculated with activated sludge while an anammox inoculum is added later. Such study concerns a reference case ($RUN1$) where the addition time of anammox inoculum $\bar{t}$ is set to $10 \ d$ and the anammox inoculum size $\psi^*_{AMX, \bar{t}}$ is set to $500 \ gCOD \ m^{-3}$. Numerical results are summarized in Figures 1–4.

Figure 1 shows the time evolution of soluble substrates and planktonic biomasses within the bulk liquid during the initial 3 days. Initially, biofilm granules have small size and the dynamics of substrates are governed by planktonic biomass. In particular, under aerobic conditions nitritation and $NO_2$ and $OC$ oxidation occur in the reactor due to the metabolic activities of planktonic $AOB$, $NOB$ and $HB$. Ammonium $S^*_{NH_4}$ (blue in Figure 1-top) is consumed and converted into nitrite $S^*_{NO_2}$ (red in Figure 1-top) by planktonic $AOB$ $\psi^*_{AOB}$ (nitritation process), and subsequently nitrite is converted into nitrate $S^*_{NO_3}$ (yellow in Figure 1-top) by planktonic $NOB$ $\psi^*_{NOB}$ (nitrite oxidation). Meantime, the consumption of organic carbon $S^*_{OC}$ (cyan in Figure 1-top) indicates the activity of planktonic $HB$ $\psi^*_{HB}$. This initial trend is followed by a turnover phase in which the planktonic biomass rapidly decreases (Figure 1-bottom), due to the granulation process and the low hydraulic retention time (HRT), and is replaced by the sessile biomass of growing biofilm granules. Nevertheless, when the planktonic species concentrations approach to zero the amount of grown sessile biomass is still low to compensate the lost contribution of substrates conversion by planktonic biomass, and sudden increases in ammonium and organic carbon and reductions in nitrite and nitrate are observed (according to the composition of the influent wastewater). In real granular-based plants, the decrease of planktonic biomass is sometimes slowed down by considering variable HRTs or loading rates ^[10,13]. Although these procedures have been not included in the model, this does not compromise its reliability in describing the successive biological processes and substrates dynamics.

Figure 2 shows the time evolution of soluble substrates and planktonic $AMX$ concentrations, and sessile masses within the system, until a steady-state configuration is reached. After the initial days, the variation of the substrates concentration due to the sessile metabolism begins to be visible. In particular, since granules are still small, aerobic conditions are found in almost all the biofilm domain and dynamics of substrates are governed by aerobic sessile species: $AOB$ convert ammonium (blue in Figure 2-top) into nitrite (red in Figure 2-top) and $HB$ oxidate organic compounds (cyan in Figure 2-top). On the contrary, the conversion of nitrite to nitrate (yellow in Figure 2-top) by $NOB$ is not visible. This happens because $NOB$ have a high $O_2$ affinity constant and are less competitive than $AOB$ and $HB$ at low oxygen levels. Consequently, high masses of sessile $AOB$ and $HB$ and negligible amounts of $NOB$ are observed in the first 30 days (Figure 2-bottom). Furthermore, no anammox biomass is detected in the reactor until $\bar{t} = 10 \ d$, when the anammox inoculum is added and a discontinuity is generated in the graph of planktonic $AMX$ concentration (Figure 2-top). Then, planktonic $AMX$ invade the innermost layers of granules, where anoxic conditions optimal for their anaerobic metabolism are found and begin to grow in sessile form. The nitritation process by $AOB$ lead to the partial removal of ammonium, which reaches a temporary equilibrium at about 50% of the influent concentration, while all the organic matter is oxidized aerobically by $HB$. At this moment, a very long transition phase begins, in which ammonium and nitrite concentrations remain almost constant. Granules are fully developed and present anoxic conditions and shortage of organic carbon in the internal layers, which inhibit the $AOB$ and $HB$ growth.

At the same time, such anoxic conditions and the simultaneous presence of ammonium and nitrite in the reactor promote the $AMX$ growth. However, as the $AMX$ biomass is characterized by very low growth rates, the further ammonium and nitrite consumption induced by the anammox process begins to be relevant after 100 days. Specifically, a considerable growth of anammox biomass is observed between 100 and 200 days and leads to the consumption of ammonium and nitrite and small production of nitrate. The steady-state configuration shows a residual ammonium concentration lower than $40-50 \ gCOD \ m^{-3}$, and very low concentrations of nitrite and nitrate. In conclusion, the system presents an ammonium removal efficiency of about $90\%$, which is achieved through a two-stage treatment process: the first stage is governed by $AOB$ which halve the ammonium concentration and produce nitrite necessary for the successive stage; the second stage is governed by $AMX$ which consume a further relevant amount of ammonium by using nitrite.

In Figure 3, the distribution of sessile species within the granule is shown at 15, 50, 120, 150 and 300 days. At $T = 15 \ d$, the granule is constituted mostly by $HB$ (cyan), which have high growth rates, and $AOB$ (blue), favored by the availability of ammonium and oxygen. At $T = 50 \ d$, the granule is fully developed and is characterized by internal anoxic zones.

Therefore, aerobic $AOB$ and $HB$ accumulate in the outermost layers of the granule, while $AMX$ (red) grow in the centre. However, due to the low growth rate, a small $AMX$ core only begins to be visible at $T = 120 \ d$. Their growth intensifies strongly in the successive phases, up to a steady-state configuration where the granule is dominated by $AMX$ while $AOB$ and $HB$ are limited to the thin outermost layer. The amount of $NOB$ (yellow) present in the granule is negligible throughout the process. Such result is in agreement with Vlaeminck et al. ^[48], which shows that low oxygen levels limit the $NOB$ metabolic pathways. Indeed, as explained above, $NOB$ are less competitive than $AOB$ and $HB$ in the presence of low oxygen concentrations. This is beneficial for the ammonium removal process because $NOB$ metabolic activities would include the consumption of nitrite necessary for the $AMX$ growth.

The microbial distribution and the relative abundance here reported are in agreement with experimental observations in literature ^[48,49,50]. However, they can be influenced by various factors, such as the granule size, the oxygen level in the reactor, the influent composition. First, the $AMX/AOB$ ratio within the granules varies due to the granule size and the oxygen level. Specifically, large granules (described in this study) have more extended anoxic zones and consequently are characterized by an high $AMX/AOB$ ratio, while smaller granules have reduced anoxic cores and lower $AMX/AOB$ ratios ^[48,50]. Then, the extension of the anoxic zone and therefore the $AMX/AOB$ ratio increases as the oxygen level set in the reactor decreases ^[25,28]. Furthermore, the influent composition can affect the evolution of granules and the relative abundance of individual species. For example, as shown in Mozumder et al. ^[34], the presence of organic substance within the influent promotes the heterotrophic growth in the granules.

Lastly, the evolution of the granule radius over time $R(t)$ is reported in Figure 4. $R(t)$ starts from a vanishing initial value and it increases rapidly during the first days, when the granulation process is more intense. It reaches a temporary constant value after 30 days, and increases slightly again due to the anammox growth. The final steady-state value is approximately $1 \ mm$.

4.3.   S2 - Effects of the anammox addition time on the process

Results presented in Section 4.2 describe the evolution and dynamics of anammox granules in a granular-based reactor inoculated initially with a nitrifying/denitrifying activated sludge and later with an anammox sludge. The addition time of the anammox inoculum appears to be a key element in the process, since the development of the anammox biomass strongly depends on the conditions inside the reactor at the addition time.

Specifically, different addition times can lead to different scenarios. In such context, this study (S2) focuses on the effects that the addition time of anammox sludge $\bar{t}$ has on the growth of anammox granules and on the duration of the start-up period. For this purpose, six simulations ($RUN2 - RUN7$) have been carried out by setting the addition time $\bar{t}$ equal to 0, 3, 5, 10, 20 and 40 days, respectively. The inoculum size $\psi^*_{AMX, \bar{t}}$ has been fixed equal to $500 \ gCOD \ m^{-3}$ for all simulations. Results are reported in Figures 5–9.

Figure 5 shows the oxygen concentration within the granule at $\bar{t}$ for simulations $RUN3 - RUN7$. As can be seen, the extension of the anoxic zone increases as $\bar{t}$ increases. Indeed, when the anammox inoculum is introduced in the initial phase, granules are small and oxygen fully penetrated. On the other hand, when the inoculum is added later, granules are fully developed and almost completely anoxic, except for the most external layers. The extension of the anoxic zone at $\bar{t}$ affects the invasion and the growth of the anammox biomass.

This is visible in Figure 6, where the distribution of anammox sessile biomass at different times is reported for simulations $RUN2$, $RUN4$, $RUN5$, $RUN7$. When the $AMX$ addition occurs at the beginning of the process ($RUN2$), the invasion process is inhibited by the presence of oxygen throughout the biofilm. Hence the anammox growth is very slow. Conversely, when $\bar{t}$ is high ($RUN4$, $RUN5$, $RUN7$) anammox cells colonize the anoxic granule core and grow more rapidly in sessile form. Anyway, although the addition of the anammox inoculum can be delayed to promote more intense invasion phenomena, this also leads to a delay in the initiation of the anammox process. This means that the effects of $\bar{t}$ on the anammox growth are not unique. Indeed, the complete development of the anammox biomass for $\bar{t} = 10 \ d$ ($RUN5$) appears earlier than for $\bar{t} = 40 \ d$ ($RUN7$).

Since $\bar{t}$ influences the anammox growth, it affects also the rate of the biological processes within the reactor. The concentrations of soluble substrates and $AMX$ planktonic biomass within the reactor are reported in Figure 7. As can be seen, the first stage of the process (first 100 days) does not depend on $\bar{t}$, since the amount of anammox sessile biomass is limited and plays a negligible role. The successive trends of ammonium $S^*_{NH_4}$ (blue) and nitrite $S^*_{NO_2}$ (red) concentrations are influenced by the anammox biomass and, therefore, by $\bar{t}$, while not significant variations of nitrate $S^*_{NO_3}$ (yellow) and organic matter $S^*_{OC}$ (cyan) occur. When the anammox growth is faster, ammonium and nitrite are more rapidly consumed and the time necessary to reach the steady-state configuration decreases. Specifically, the minimum duration is achieved for $\bar{t} = 20 \ d$, while the start-up in the case of $\bar{t} = 0$ is not yet completed after 260 days. Furthermore, $\bar{t}$ does not affect the planktonic $AMX$ wash out after the addition. No variation is observed for planktonic $AOB$, $NOB$ and $HB$ concentrations (data not shown).

Figure 8 displays the relative abundances of active microbial within the granule at different times. It confirms the results previously described: the growth process of anammox biomass is deeply affected by $\bar{t}$ and is faster in the case of $\bar{t} = 20 \ d$ ($RUN 6$). However, in all cases except $\bar{t} = 0$ ($RUN 2$), the steady-state microbial distribution is reached within the observation period of 260 days.

Lastly, the evolution of the granule radius $R(t)$ over time is shown in Figure 9. $\bar{t}$ affects the granule evolution in the second stage of the process (Figure 9-right), since the further radius increase around 100-180 days is associated to the $AMX$ growth. However, the steady-state granule dimension achieved is not dependent on $\bar{t}$.

From all these results, it is clear that the addition time of the anammox inoculum $\bar{t}$ influences the rate of biological processes related to the $AMX$ biomass and therefore the duration of the start-up period. Anyway, it does not affect the steady-state configuration and the removal efficiency of the process.

4.4.   S3 - Effects of the anammox inoculum size on the process

As explained in Section 1, the most interesting strategy to carry out the start-up of a PN/A granular reactors is based on the use of anammox sludge inocula.

However, although this strategy allows to significantly reduce the start-up duration, it has high costs related to the limited availability of anammox biomass around the world. Therefore, in order to optimize timing and costs, an exhaustive knowledge of effects of anammox inoculum and its size on the process start-up is needed. To this aim, the present study (S3) analyzes numerically the role that the anammox inoculum size $\psi^*_{AMX, \bar{t}}$ plays in the process and investigates how this size affects the duration of the start-up period. Six simulations ($RUN8 - RUN13$) have been carried out with different values of $\psi^*_{AMX, \bar{t}}$ ($100, 300, 500, 700, 1000, 1500 \ gCOD \ m^{-3}$), while the addition time of anammox sludge $\bar{t}$ is fixed to $10 \ d$. Results are summarized in Figures 10-13.

Figure 10 reports the distribution of anammox species within the granule for $\psi^*_{AMX, \bar{t}}$ equal to $100$, $500$, $700$, $1000 \ gCOD \ m^{-3}$, at different times. Since the invasion process is proportional to the concentration of $AMX$ planktonic biomass, it becomes more intense as the $AMX$ inoculum size increases $\psi^*_{AMX, \bar{t}}$. Then, by increasing the $AMX$ inoculum size, the anammox sessile biomass grows faster and the steady-state microbial distribution is reached earlier. However, $\psi^*_{AMX, \bar{t}}$ has no effects on the steady-state configuration.

Trends of soluble substrates and $AMX$ planktonic biomass within the reactor are reported in Figure 11. The $AMX$ inoculum size influences the rate of anammox processes and then, the trend of substrates involved. The first stage seems to be independent of $\psi^*_{AMX, \bar{t}}$, since anammox processes are negligible and the grown $AMX$ biomass is not still sufficient to significantly affect the substrates concentration. In the second stage, anammox processes intensifies and the consumption rates of ammonium $S^*_{NH_4}$ (blue) and nitrite $S^*_{NO_2}$ (red) varies with the variation of $\psi^*_{AMX, \bar{t}}$. As the addition time of the $AMX$ inoculum $\bar{t}$, the $AMX$ inoculum size does not affect the steady-state concentrations of substrates and $AOB$, $NOB$ and $HB$ planktonic biomasses (data not shown).

Figure 12 displays the relative abundances of active microbial species within the granule at different times. It confirms that as $\psi^*_{AMX, \bar{t}}$ increases, the anammox growth is faster and the steady-state distribution is achieved earlier.

Figure 13 reports the evolution of the granule radius $R(t)$ over time for different $\psi^*_{AMX, \bar{t}}$. As noted above, an increase in $\psi^*_{AMX, \bar{t}}$ leads to a faster anammox growth and thus a faster increase in $R(t)$.

From the results shown it is possible to draw a general conclusion: the $AMX$ growth rate and the duration of the process start-up are directly proportional to the $AMX$ inoculum size. Anyway, it has no effect on the steady-state configuration occurring inside the reactor.

4.5.   S4 - Combined effects of anammox addition time and anammox inoculum size

The last numerical study (S4) analyzes the combined effect of the addition time $\bar{t}$ and the anammox inoculum size $\psi^*_{AMX, \bar{t}}$ on the biological process, in order to show how the choice of these operating parameters affects the process start-up of PN/A granular bioreactors. For this purpose, 25 simulations are performed by varying $\bar{t}$ and $\psi^*_{AMX, \bar{t}}$ (values reported in Table 2). The duration of the process start-up is assumed equal to the time necessary to achieve the steady-state ammonium concentration $T^ *$.

Figure 14 shows the anammox sessile mass $m_{AMX}$ for different $\bar{t}$ and $\psi^*_{AMX, \bar{t}}$. These results are shown for $T = 170 \ d$, which is the minimum time, among all the simulations carried out, necessary to reach the steady-state anammox sessile mass ($RUN25$: $\bar{t} = 20 \ d$; $\psi^*_{AMX, \bar{t}} = 1500 \ gCOD \ m^{-3}$). As the $\psi^*_{AMX, \bar{t}}$ increases, $m_{AMX}$ at $T = 170 \ d$ increases. However, $m_{AMX}$ has a less than linear behaviour with increasing $\psi^*_{AMX, \bar{t}}$. On the other hand, $m_{AMX}$ does not present a unique trend as $\bar{t}$ varies. Specifically, it increases as $\bar{t}$ increases up to $20 \ d$, while it decreases again for $\bar{t} > 20 \ d$. As explained in Section 4.3, this happens because the further delay in the addition of the $AMX$ inoculum prevails over the further acceleration in the process of invasion and anammox growth related to the increase of the anoxic zone within the granule.

The duration of the start-up $T^*$ for different $\bar{t}$ and $\psi^*_{AMX, \bar{t}}$ is reported in Figure 15. As expected, $T^*$ presents a behavior inversely proportional to the anammox mass $m_{AMX}$ at $T = 170 \ d$. When the invasion and growth of anammox biomass are more intense, $m_{AMX}$ grows faster. As a result, ammonium consumption is faster and the steady state value is achieved earlier.

It seems evident that the duration of the start-up can be reduced by increasing the anammox inoculum size and by adding it when most part of granules is under anoxic conditions. However, some considerations need to be done. First, higher is the anammox inoculum size, higher is the total cost. Furthermore, as can be seen in Figure 15, the time saved for the start-up is reduced as the inoculum size increases. For example, considering a 50% increase of $\psi^*_{AMX, \bar{t}}$ from $1000$ to $1500 \ gCOD \ m^{-3}$ (dash-dot line and thick line, respectively), the start-up period reduces only by a few days. Finally, the anammox inoculum should be added into the reactor at the time that maximizes the rates of anammox growth and ammonium removal.

5.   Conclusions

In the present work, we propose a mathematical model on the combined PN/A granular system. Such model describes the de novo anammox granulation through a spherical free boundary domain with radial symmetry, and simulates the start-up process of this wastewater system by considering both the biofilm and bulk liquid dynamics. In particular, growth and decay of sessile and planktonic biomasses, diffusion and conversion of substrates, invasion phenomena of planktonic biomass, attachment and detachment phenomena are modelled. A planktonic anammox inoculum is supposed to be added in the system at a specific time. Then, microbial invasion is accounted to initiate the anammox sessile growth.

Numerical results show that the evolution and composition of the biofilm granule mainly depend on oxygen trends. When the granule is small, high oxygen concentration throughout the granule promote the growth of aerobic species ($AOB$ and $HB$ bacteria). On the contrary, when the granule reaches larger dimensions, anoxic conditions arising in the granule center favor the anaerobic metabolic activity of $AMX$ bacteria. In the final steady-state configuration, granules are dominated by $AMX$ bacteria and small amounts of $AOB$ and $HB$ bacteria are detected in the outermost layers. Furthermore, $NOB$ bacteria are negligible throughout the process. The 90% of ammonium is removed through a two-stage treatment process: aerobic oxidation of ammonium and organic carbon by $AOB$ and $HB$; anaerobic ammonium consumption by $AMX$ bacteria.

Moreover, numerical studies have been carried out to investigate the effects of two significant factors on the start-up process: the anammox inoculum size and the anammox addition time. Such parameters affect the anammox invasion and growth, and then, the duration of the start-up period. The combined effect of the addition time and the anammox inoculum size on the start-up process is finally investigated, with the aim of providing a numerical support in the choice of the start-up strategy. However, it should be noted that such model results are qualitative, and a calibration procedure should be carried out in order to provide engineering relevance to the model output.

Acknowledgments

Fabiana Russo's research activity has been conducted in the context of D.D. n. 155 on 17 May 2018 additional PhD fellowships for 2018/2019 academic year, course XXXIV within the framework of POR Campania FSE 2014-2020 ASSE III - Specific objective 14 Action 10.5.2 - Public notice "Innovative PhD with industrial characterization".

Authors would like to acknowledge the project "Employing circular economy approach for OFMSW management within the Mediterranean countries - CEOMED" number A_B.4.2_0058, funded under the ENI CBC MED 2014-2020 programme and the project "METAGRO - bioMETanazione dei sottoprodotti della filiera AGROindustriale campana" CUP: B18H19005240009.

This paper has been performed under the auspices of the G.N.F.M. of I.N.d.A.M.

Conflict of interest

The authors declare there is no conflict of interest.

Supplementary material

According to D'Acunto et al. ^[46], the specific growth rates due to sessile species $r_{M, i}$ are modelled as Monod-type kinetics:

where $\mu_{\max, i}$ is the maximum growth rate for biomass $i$, $K_i$ is the affinity constant of the consumed substrate for biomass $i$, and $k_{d, i}$ is the decay constant for biomass $i$, $\beta_1$ and $\beta_2$ are the reduction factors for denitrification.

The formation rate of inactive biomass is given by the sum of decay rates of each active species, modelled as first order kinetics:

The specific growth rates due to planktonic species $r_i$, with $i \in \{AOB, AMX,$ $NOB, HB\}$, are defined as:

where $k_{col, i}$ is the maximum colonization rate of motile species $i$.

The conversion rates of planktonic species $r_{\psi, i}$ are expressed by:

where $Y_{\psi, i}$ denotes the yield of non-motile species $i$ on the corresponding motile species. While, the conversion rates for soluble substrates within the biofilm $r_{S, j}$, with $j \in \{NH_4, NO_2, NO_3, OC, O_2\}$, are listed below ^[46]:

where $Y_{AOB}$ is the yield of $AOB$ on $NH_4$, $Y_{AMX}$ is the yield of $AMX$ on $NH_4$, $Y_{NOB}$ is the yield of $NOB$ on $NO_2$ and $Y_{HB}$ is the yield of $HB$ on $OC$.

Moreover, the conversion rates of planktonic biomasses $r^*_{\psi, i}$ within the bulk liquid are defined as:

while, the conversion rates of soluble substrates $r^*_{S, j}$, with $j \in \{NH_4, NO_2, NO_3, OC\}$, within the bulk liquid are listed below:

The values used for all stoichiometric and kinetic parameters are reported in Supplementary Table S3.