Approximate Models of Microbiological Processes in a Biofilm Formed on Fine Spherical Particles

Abstract

This paper concerns the dynamical modeling of the microbiological processes that occur in the biofilms that are formed on fine inert particles. Such biofilm forms e.g. in fluidized-bed bio-reactors, expanded bed biofilm reactors and biofilm air-lift suspension reactors. An approximate model that is based on the Laplace–Carson transform and a family of approximate models that are based on the concept of the pseudo-stationary substrate concentration profile in the biofilm were proposed. The applicability of the models to the microbiological processes was evaluated following Monod or Haldane kinetics in the conditions of dynamical biofilm growth. The use of approximate models significantly simplifies the computations compared to the exact one. Moreover, the stiffness that was present in the exact model, which was solved numerically by the method of lines, was eliminated. Good accuracy was obtained even for large internal mass transfer resistances in the biofilm. It was shown that significantly higher accuracy was obtained using one of the proposed models than that which was obtained using the previously published approximate model that was derived using the homotopy analysis method.

1. Introduction

Numerous problems in the field of chemical and process engineering, for example those concerning heat transfer or simultaneous diffusion and reaction, are described by nonlinear differential equations [1,2,3]. Typically, such differential problems are solved using numerical algorithms, which may cause problems with convergence and require time-consuming computations [4,5,6]. For that reason, there is an interest in using ap-proximation methods for solving nonlinear differential equations. The application of these methods in dynamic simulations of two-scale objects is particularly significant [7]. Fluidized-bed catalytic reactors, reactors with a stationary catalytic bed and all microbiological reactors with an immobilized biofilm belong to this group [8]. In such devices, the microbiological process occurs in the liquid phase where the diffusional resistances are negligible and in the biofilm where the internal mass transfer resistances cannot be neglected. As a result, the mathematical description of the process is based on the application of different types of equations for both of the phases. According to Nicolella et al. [9], the significance of particle-based biofilm reactors is large as they provide the potential to develop compact and high-rate processes.

Studies concerning approximation models of catalytic processes in porous catalysts can be found in the literature. An example is the study by Szukiewicz [10], who proposed a method based on the Laplace–Carson transform for the determination of the dynamic changes of the average reactant concentration in a catalytic pellet. Another example is the study by Sun et al. [11], who used the Adomian decomposition method (ADM) for the determination of the stationary profiles of the reactant concentration in a catalytic pellet and its effectiveness factor. Abbasbandy [12] used the homotopy analysis method (HAM) for the modelling of a catalytic pellet at stationary conditions. In our previous study [13], the homotopy analysis method was applied for the simulation of a biofilm at a steady state.

Several studies can be found in which a simplified approach was used in the mod-eling of bioreactors with a biofilm. For example, Russo et al. [14] applied the concept of a so-called zero-dimensional biofilm model for a dynamic simulation of a continuous biofilm bioreactor. The quoted authors introduced a lumped-state variable expressing the overall attached biomass concentration in the bioreactor. This is a far-reaching approach, because the external and internal mass transfer resistances of the growth-limiting substrate are neglected. A similar approach was applied in the study by Dokianakis et al. [15] and in other recent studies [16,17,18,19] concerning the modeling of biofilm bioreactors. It was repeatedly proven that mass transfer resistances cause the distribution of a substrate concentrations inside the biofilm [20,21]. Neglecting this phenomenon in the mathematical model leads to obtaining results that are inconsistent with experimental observations [22]. Despite numerous papers concerning the applica-tion of approximation methods in the modeling of catalytic pellets and several studies concerning the modeling of biofilms at stationary conditions [13,23], there are no studies which demonstrate the applicability of such an approach in the modeling of dynamical biofilm growth.

This study presents an approximate model for biofilm dynamics, which is based on the Laplace–Carson transform and a family of approximate models that are based on the concept of pseudo-stationary substrate concentration profiles in the biofilm. The com-putational efficiency and accuracy of the proposed models were evaluated. The solution of the exact model that was obtained by using the method of lines was chosen as a reference. Moreover, the accuracy of the approach presented was compared with a previously published model, which was derived using the homotopy analysis method.

2. Exact Mathematical Model

A microbiological process limited by a single substrate is considered. Mass balance of substrate A in a spherical biofilm is the following:

$$\frac{\partial c_{\mathrm{A}}^b}{\partial t}=D_{e\mathrm{A}}\frac{{\partial }^2{c}_{\mathrm{A}}^b}{\partial x^2}+D_{e\mathrm{A}}\frac{2}{x}\frac{\partial c_{\mathrm{A}}^b}{\partial x}-r_{\mathrm{A}}^b(c_{\mathrm{A}}^b)$$

(1)

Equation (1) is subject to the following boundary and initial conditions:

$$\frac{\partial c_{\mathrm{A}}^b(r_0,t)}{\partial x}=0,\hspace{1em}\hspace{1em}t\ge 0$$

(2)

$$c_{\mathrm{A}}^b(r_0+L_b,t)=c_{\mathrm{A}s}^b(t)$$

(3)

$$c_{\mathrm{A}}^b\left(x,0\right)=c_{0\mathrm{A}}^b(x)\hspace{1em}x\in [r_0,r_0+L_b]$$

(4)

where $c_{0\mathrm{A}}^b$ denotes the initial substrate A concentration in the biofilm, while $c_{\mathrm{A}s}^b$ de-notes the concentration at the biofilm’s surface.

The uptake rate of the carbonaceous substrate A in the biofilm can be formulated as follows:

$$r_{\mathrm{A}}^b(c_{\mathrm{A}}^b)=\frac{1}{w_{\mathrm{BA}}}\cdot f(c_{\mathrm{A}}^b){\rho }_b$$

(5)

In this study, two kinetic models of microbiological processes were taken into account, i.e., Monod and Haldane. The former describes the utilization of many non-toxic substrates [24], whereas the latter is commonly used in modelling the biodegradation of toxic compounds [25]. Function $f(c_{\mathrm{A}}^b)$ for Monod kinetics is defined as

$$f(c_{\mathrm{A}}^b)=\frac{k\cdot c_{\mathrm{A}}^b}{K_{\mathrm{A}}+c_{\mathrm{A}}^b}$$

(6)

For Haldane kinetics, one obtains

$$f(c_{\mathrm{A}}^b)=\frac{k\cdot c_{\mathrm{A}}^b}{K_{\mathrm{A}}+c_{\mathrm{A}}^b+\frac{{(c_{\mathrm{A}}^b)}^2}{K_{in}}}$$

(7)

The following dimensionless space coordinate in the biofilm was introduced to the mathematical model (1)–(2):

$$z=\frac{x-r_0}{L_b(t)}\in [0,1]$$

(8)

where L_b(t) denotes the biofilm thickness variable over time.

After that, the substrate concentration was described by a new function which depends on the dimensionless space coordinate and time:

$$c_{\mathrm{A}}^b(x,t)\to {\tilde{c}}_{\mathrm{A}}^b(z(x,t),t)$$

(9)

In accordance with differential calculus, the derivatives of the substrate concentration with respect to x and t were transformed according to the following relationships:

$$\frac{\partial c_{\mathrm{A}}^b}{\partial x}=\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial z}\frac{\partial z}{\partial x}=\frac{1}{L_b}\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial z}$$

(10)

$$\frac{\partial c_{\mathrm{A}}^b}{\partial t}=\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial z}\frac{\partial z}{\partial t}+\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial t}=-\frac{(x-r_0)L_{b,t}^{\prime }}{L_b^2}\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial z}+\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial t}=-\frac{z{L}_{b,t}^{\prime }}{L_b}\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial z}+\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial t}$$

(11)

The dynamic changes of a substrate concentration in the biofilm can be, therefore, described by the following partial differential equation (see Appendix A for the derivation):

$$\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial t}=\frac{D_{e\mathrm{A}}}{L_b^2}\left[\frac{{\partial }^2{\tilde{c}}_{\mathrm{A}}^b}{\partial z^2}+\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial z}\left(\frac{2L_b}{(r_0+z{L}_b)}+\frac{1}{D_{e\mathrm{A}}}(z{L}_b)L_{b,t}^{\prime }\right)-\frac{L_b^2}{D_{e\mathrm{A}}}{r}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}}^b)\right]$$

(12)

The equation of the biofilm thickness dynamics was based on the biomass balance from the paper [26]:

$$\frac{d{L}_b}{dt}=\frac{1}{3{\rho }_b}\frac{r_b^3-r_0^3}{r_b^2}{r}_{\mathrm{B},av}^b-\frac{k_{det}{L}_b}{3}\frac{r_b^3-r_0^3}{r_b^2}$$

(13)

where:

$$r_{\mathrm{B},av}^b=\frac{3}{r_b^3-r_0^3}{\displaystyle \underset{r_0}{\overset{r_b}{\int }}{x}^2{r}_{\mathrm{B}}^b(c_{\mathrm{A}}^b})dx=\frac{3L_b}{r_b^3-r_0^3}{\displaystyle \underset{0}{\overset{1}{\int }}{(r_0+L_{b}z)}^2{r}_{\mathrm{B}}^b(c_{\mathrm{A}}^b})dz$$

(14)

with the initial condition:

$$L_b(0)=L_{b0}$$

(15)

3. Approximate Mathematical Model Based on Laplace-Carson Transform

Numerical simulations were performed with the use of the exact model of biofilm growth (Equation (12) with conditions (2)–(4) and Equation (13) with the condition (15)). The computations were realized using Matlab 2016a. The results proved that the term $\frac{1}{D_{e\mathrm{A}}}(z{L}_b)L_{b,t}^{\prime }$ is lower than $\frac{2L_b}{(r_0+z{L}_b)}$ by more than three orders of magnitude. The former term can be therefore neglected. Then, Equation (12) simplifies to:

$$\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial t}=\frac{D_{e\mathrm{A}}}{L_b^2}\left[\frac{{\partial }^2{\tilde{c}}_{\mathrm{A}}^b}{\partial z^2}+\frac{2L_b}{(r_0+z{L}_b)}\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial z}-\frac{L_b^2}{D_{e\mathrm{A}}}{r}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}}^b)\right]$$

(16)

In order to reduce the number of parameters in the above equation, the variable τ was introduced, which is related to time t according to the relationship:

$$d\tau =\frac{D_{e\mathrm{A}}}{L_b^2}dt$$

(17)

After introducing variable τ, Equation (16) will transform to:

$$\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial \tau }=\frac{{\partial }^2{\tilde{c}}_{\mathrm{A}}^b}{\partial z^2}+\frac{2L_b}{(r_0+L_{b}z)}\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial z}-\frac{L_b^2}{D_{e\mathrm{A}}}{r}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}}^b)$$

(18)

The following boundary and initial conditions are associated with Equation (18):

$$\frac{\partial {\tilde{c}}_{\mathrm{A}}^b(0,\tau )}{\partial z}=0$$

(19)

$${\tilde{c}}_{\mathrm{A}}^b(1,\tau )={\tilde{c}}_{\mathrm{A}s}^b$$

(20)

$${\tilde{c}}_{\mathrm{A}}^b(z,0)={\tilde{c}}_{\mathrm{A}0}^b$$

(21)

Equation (18) can be written as:

$$\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial \tau }=\frac{{\partial }^2{\tilde{c}}_{\mathrm{A}}^b}{\partial z^2}+\frac{2L_b}{(r_0+L_{b}z)}\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial z}-\Phi {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}}^b)$$

(22)

where: $\Phi =\frac{L_b^{2}k}{D_{e\mathrm{A}}}$

The forms of the kinetic expressions in Equation (22) were as follows:

• Monod model:

$${\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}}^b)=\frac{1}{w_{\mathrm{BA}}}\frac{{\tilde{c}}_{\mathrm{A}}^b}{K+{\tilde{c}}_{\mathrm{A}}^b}{\rho }_b$$

(23)

• Haldane model:

$${\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}}^b)=\frac{1}{w_{\mathrm{BA}}}\frac{{\tilde{c}}_{\mathrm{A}}^b}{K+{\tilde{c}}_{\mathrm{A}}^b+{({\tilde{c}}_{\mathrm{A}}^b)}^2/K_{in}}{\rho }_b$$

(24)

In order to derive the approximate model, the kinetic expression was linearized using Taylor series expansion:

$${\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}}^b)={\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)+\frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}({\tilde{c}}_{\mathrm{A}}^b-{\tilde{c}}_{\mathrm{A}s}^b)$$

(25)

After introducing Expression (25) into Equation (22), one obtains:

$$\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial \tau }=\frac{{\partial }^2{\tilde{c}}_{\mathrm{A}}^b}{\partial z^2}+\frac{2L_b}{r_0+L_{b}z}\frac{\partial {\tilde{c}}_{\mathrm{A}}^b}{\partial z}-\Phi {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)+\Phi \frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}({\tilde{c}}_{\mathrm{A}s}^b-{\tilde{c}}_{\mathrm{A}}^b)$$

(26)

The Laplace–Carson transform was used for the derivation of the approximate model. It is an integral transform defined by the following formula:

$$Q(p,x)=p\cdot {\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{-p\tau }(\tau ,x)d}\tau$$

(27)

It has, historically, been applied in solving various mathematical and physical problems. For example, Kumar and Qureshi [27] used it to obtain an analytical solution for the non-integer order initial value problems with the Caputo operator. Makarov [28] applied it in order to reduce a model of nonstationary flow of a viscoplastic medium. Aggarwal et al. [29] proved that the Laplace–Carson transform can be used for solving system of convolution-type Volterra integro-differential equations of the first kind. In turn, Szukiewicz published several papers describing the application of this transform to the problems of diffusion and reaction in a catalytic pellet [10,30,31]. His studies were the motivation for developing an approximate model of the process occurring in a biofilm.

Equation (26) is converted by the Laplace–Carson transform into:

$$pC-p{\tilde{c}}_{\mathrm{A}s}^b=\frac{{\partial }^{2}C}{\partial z^2}+\frac{2L_b}{r_0+L_{b}z}\frac{\partial C}{\partial z}-\Phi {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)+\Phi \frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}({\tilde{c}}_{\mathrm{A}s}^b-C)$$

(28)

with boundary conditions:

$$\frac{dC(0)}{dz}=0$$

(29)

$$C(1)={\tilde{c}}_{\mathrm{A}s}^b$$

(30)

where:

$$C=p{\displaystyle \underset{0}{\overset{\infty }{\int }}{\tilde{c}}_{\mathrm{A}}^b{e}^{-p\tau }d\tau }$$

(31)

The solution of Equation (28) with boundary conditions (29) was the following:

$$\begin{array}{l}C=\frac{1}{{\mu }^3}\left({\tilde{c}}_{\mathrm{A}0}^{b}p-\Phi {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)+\Phi \frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}{\tilde{c}}_{\mathrm{A}s}^b\right)+\\ \frac{\left(L_b+r_0\right)\left({\tilde{c}}_{\mathrm{A}s}^{b}p-{\tilde{c}}_{\mathrm{A}0}^{b}p+\Phi {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)\right)\left({\mathrm{e}}^{\mu \left(z+1\right)}\left(L_b+r_0\mu \right)-{\mathrm{e}}^{-\mu \left(z-1\right)}\left(L_b-r_0\mu \right)\right)}{{\mu }^2\left(r_0+L_{b}z\right)\left(r_0\mu -L_b+{\mathrm{e}}^{2\mu }(L_b+r_0\mu )\right)}\end{array}$$

(32)

where:

$$\mu =\sqrt{p+\Phi \frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}}$$

(33)

The average concentration ${\tilde{c}}_{\mathrm{A},av}^b$ was calculated as an integral average for a spherical bioparticle:

$${\tilde{c}}_{\mathrm{A},av}^b=\frac{3L_b}{{(r_0+L_b)}^3-r_0^3}{\displaystyle \underset{0}{\overset{1}{\int }}{\tilde{c}}_{\mathrm{A}}^b(z)}{(r_0+L_{b}z)}^{2}dz$$

(34)

It can be written in the domain of complex numbers as:

$$C_{av}=\frac{3L_b}{{(r_0+L_b)}^3-r_0^3}{\displaystyle \underset{0}{\overset{1}{\int }}C(p,z)}{(r_0+L_{b}z)}^{2}dz$$

(35)

Equation (32) was substituted to Equation (35). Afterwards, using Matlab’s symbolic toolbox, it was integrated and transformed into the following form:

$$p{C}_{av}-p{\tilde{c}}_{\mathrm{A}0}^b=\Psi (p)({\tilde{c}}_{\mathrm{A}s}^b-C_{av})-\Phi {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)+\Phi \frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}({\tilde{c}}_{\mathrm{A}s}^b-C_{av})$$

(36)

where:

$$\Psi (p)=\frac{3{\mu }^2\left(L_b+r_0\right)\left(L_b^2\mu \cosh \left(\mu \right)-L_b^2\sinh \left(\mu \right)+r_0^2{\mu }^2\sinh \left(\mu \right)+L_b{r}_0{\mu }^2\sinh \left(\mu \right)\right)}{\left(\begin{array}{l}3L_b^3\sinh \left(\mu \right)+3L_b^2{r}_0\sinh \left(\mu \right)-3L_b^3\mu \cosh \left(\mu \right)+L_b^3{\mu }^2\sinh \left(\mu \right)-3r_0^3{\mu }^2\sinh \left(\mu \right)+\\ 3r_0^3{\mu }^3\cosh \left(\mu \right)-3L_b^2{r}_0\mu \cosh \left(\mu \right)+3L_b{r}_0^2{\mu }^3\cosh \left(\mu \right)+L_b^2{r}_0{\mu }^3\cosh \left(\mu \right)-3L_b{r}_0^2{\mu }^2\sinh \left(\mu \right)\end{array}\right)}$$

(37)

The left-hand side of Equation (36) is a transform of the derivative of the average substrate concentration with respect to time, while the right-hand side includes the transform of the driving force $(c_{\mathrm{A}s}^b-C_{av})$. Therefore, Equation (21) is a Laplace–Carson transform of the linear driving force (LDF) equation [10].

In order to verify the obtained expression, r₀ = 0 was substituted. It was expected that Equation (37) would then simplify to the equation that was obtained by Szukiewicz [10] for a spherical catalytic pellet. As expected, the following result was obtained:

$$\Psi (p)=\frac{3{\mu }^2\left(\mu \cosh \left(\mu \right)-\sinh \left(\mu \right)\right)}{3\sinh \left(\mu \right)+{\mu }^2\sinh \left(\mu \right)-3\mu \cosh \left(\mu \right)}$$

(38)

where: $\mu =\sqrt{p+\Phi \frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}}$

For the determination of the inverse transform of Equation (36), the value p = 0 was substituted, which is equivalent to time τ tending to infinity. The form of the coefficient Ψ in the domain of real numbers depends on the derivative $\frac{\partial {\tilde{r}}_{\mathrm{A}}({\tilde{c}}_{\mathrm{A}}^b)}{\partial {\tilde{c}}_{\mathrm{A}s}^b}$ at the biofilm surface, which can accept positive or negative values, or zero:

• for $\frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}}^b)}{\partial {\tilde{c}}_{\mathrm{A}s}^b}>0$

$$\Psi =\frac{3{\phi }^2\left(L_b+r_0\right)\left(L_b^2\phi \cosh \left(\phi \right)-L_b^2\sinh \left(\phi \right)+r_0^2{\mu }^2\sinh \left(\phi \right)+L_b{r}_0{\mu }^2\sinh \left(\phi \right)\right)}{\left(\begin{array}{l}3L_b^3\sinh \left(\phi \right)+3L_b^2{r}_0\sinh \left(\phi \right)-3L_b^3\phi \cosh \left(\phi \right)+L_b^3{\mu }^2\sinh \left(\phi \right)-3r_0^3{\phi }^2\sinh \left(\phi \right)+\\ 3r_0^3{\phi }^3\cosh \left(\phi \right)-3L_b^2{r}_0\mu \cosh \left(\phi \right)+3L_b{r}_0^2{\mu }^3\cosh \left(\phi \right)+L_b^2{r}_0{\mu }^3\cosh \left(\phi \right)-3L_b{r}_0^2{\mu }^2\sinh \left(\phi \right)\end{array}\right)}$$

(39)

where:

$$\phi =\sqrt{\Phi \frac{\partial r_{\mathrm{A}}({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}}$$

• for $\frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}=0$

$$\Psi =\frac{15\left(L_b+r_0\right)\left(L_b^2+3L_b{r}_0+3r_0^2\right)}{L_b^3+6L_b^2{r}_0+15L_b{r}_0^2+15r_0^3}$$

(40)

• for $\frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}<0$

$$\Psi =\frac{-3{\phi }^2\left(L_b+r_0\right)\left(\phi L_b^2\cos \left(\phi \right)-L_b^2\sin \left(\phi \right)+{\phi }^2{r}_0^2\sin \left(\phi \right)+{\phi }^2{L}_b{r}_0\sin \left(\phi \right)\right)}{\left(\begin{array}{l}-3L_b^3\sin \left(\phi \right)-3L_b^2{r}_0\sin \left(\phi \right)+3\phi L_b^3\cos \left(\phi \right)-{\phi }^2{L}_b^3\sin \left(\phi \right)+3{\phi }^2{r}_0^3\sin \left(\phi \right)+\\ 3{\phi }^3{r}_0^3\cos \left(\phi \right)+3\phi L_b^2{r}_0\cos \left(\phi \right)+3{\phi }^3{L}_b{r}_0^2\cos \left(\phi \right)+{\phi }^3{L}_b^2{r}_0\cos \left(\phi \right)+3{\phi }^2{L}_b{r}_0^2\sin \left(\phi \right)\end{array}\right)}$$

(41)

where:

$$\phi =\sqrt{-\Phi \frac{\partial {\tilde{r}}_{\mathrm{A}}^b({\tilde{c}}_{\mathrm{A}s}^b)}{\partial {\tilde{c}}_{\mathrm{A}}^b}}$$

Applying the inverse Laplace–Carson transform to Equation (36) results in the LDF equation, in the domain of real numbers, as follows:

$$\frac{d{c}_{\mathrm{A},av}^b}{d\tau }=\Psi (c_{\mathrm{A}s}^b-c_{\mathrm{A},av}^b)-\Phi {\tilde{r}}_{\mathrm{A}}(c_{\mathrm{A},av}^b)+\Phi \frac{\partial {\tilde{r}}_{\mathrm{A}}^b(c_{\mathrm{A}s}^b)}{\partial c_{\mathrm{A}}^b}(c_{\mathrm{A}s}^b-c_{\mathrm{A},av}^b)$$

(42)

Or, by substituting time t:

$$\frac{d{c}_{\mathrm{A},av}^b}{dt}=\frac{D_{e\mathrm{A}}}{L_b^2}\left[\Psi (c_{\mathrm{A}s}^b-c_{\mathrm{A},av}^b)-\Phi {\tilde{r}}_{\mathrm{A}}(c_{\mathrm{A},av}^b)+\Phi \frac{\partial {\tilde{r}}_{\mathrm{A}}^b(c_{\mathrm{A}}^s)}{\partial c_{\mathrm{A}}^b}(c_{\mathrm{A}s}^b-c_{\mathrm{A},av}^b)\right]$$

(43)

with the initial condition:

$$c_{\mathrm{A},av}^b(0)=c_{\mathrm{A}0}^b$$

(44)

Based on Equation (25), it can be written as follows:

$${\tilde{r}}_{\mathrm{A}}^b(c_{\mathrm{A},av}^b)\cong {\tilde{r}}_{\mathrm{A}}^b(c_{\mathrm{A}s}^b)-\frac{\partial {\tilde{r}}_{\mathrm{A}}^b(c_{\mathrm{A}s}^b)}{\partial c_{\mathrm{A}}^b}(c_{\mathrm{A}s}^b-c_{\mathrm{A},av}^b)$$

(45)

After using Equation (45), the LDF Equation (43) transformed into the following form [10]:

$$\frac{d{c}_{\mathrm{A},av}^b}{dt}=\frac{D_{e\mathrm{A}}}{L_b^2}\left[\Psi (c_{\mathrm{A}s}^b-c_{\mathrm{A},av}^b)-\Phi \cdot {\tilde{r}}_{\mathrm{A}}^b(c_{\mathrm{A},av}^b)\right]$$

(46)

with the initial condition (44).

The model given by Equation (43) with initial condition (44) is linear, whereas the model given by Equation (46) with initial condition (44) is nonlinear.

It is convenient to use approximate models in a dimensionless form. To this end, a dimensionless substrate concentration was used:

$$\eta =\frac{{\tilde{c}}_{\mathrm{A}}^b}{c_{\mathrm{A}}^c}$$

(47)

The linear approximate model in a dimensionless form can be presented as follows:

$$\frac{d{\eta }_{av}}{dt}=-\frac{1}{c_{\mathrm{A}}^c}\frac{d{c}_{\mathrm{A}}^c}{dt}\eta +\frac{D_{e\mathrm{A}}}{L_b^2}\Psi ({\eta }_s-{\eta }_{av})-\frac{D_{e\mathrm{A}}}{L_b^2}\overline{\Phi }{\tilde{r}}_{\mathrm{A}}({\eta }_{av})+\frac{D_{e\mathrm{A}}}{L_b^2}\overline{\Phi }\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }({\eta }_s-{\eta }_{av})$$

(48)

where: $\overline{\Phi }=\frac{L_b^{2}k{\rho }_b}{D_{e\mathrm{A}}{w}_{\mathrm{BA}}{c}_{\mathrm{A}}^c}$, with the initial condition:

$${\eta }_{av}(0)={\eta }_0$$

(49)

While the nonlinear one is expressed as:

$$\frac{d{\eta }_{av}}{dt}=-\frac{1}{c_{\mathrm{A}}^c}\frac{d{c}_{\mathrm{A}}^c}{dt}{\eta }_{av}+\frac{D_{e\mathrm{A}}}{L_b^2}\Psi ({\eta }_s-{\eta }_{av})-\frac{D_{e\mathrm{A}}}{L_b^2}\overline{\Phi }\cdot {\tilde{r}}_{\mathrm{A}}({\eta }_{av})$$

(50)

with the initial condition (49).

The form of kinetic expressions in a dimensionless form is as follows:

• Monod model:

$${\tilde{r}}_{\mathrm{A}}(\eta )=\frac{\eta }{\overline{K}+\eta }$$

(51)

• Haldane model:

$${\tilde{r}}_{\mathrm{A}}(\eta )=\frac{\eta }{\overline{K}+\eta +{\eta }^2/{\overline{K}}_{in}}$$

(52)

where: $\overline{K}=\frac{K_{\mathrm{A}}}{c_{\mathrm{A}}^c}$, ${\overline{K}}_{in}=\frac{K_{in}}{c_{\mathrm{A}}^c}$.

It is assumed that the average biomass growth rate can be approximated by using the following expression:

$$r_{\mathrm{B},av}^b\cong r_{\mathrm{B}}^b(c_{\mathrm{A},av}^b)$$

(53)

In order to evaluate the accuracy of this approximation, dynamical simulations of the biofilm growth were carried out using the exact model and the average biomass growth rate was calculated using the exact equation—i.e., (53). Matlab 2016a was used for this purpose. The results of the comparison are illustrated in Figure 1 for two values of the bulk liquid substrate concentration for Monod kinetics (Figure 1a) and Haldane kinetics (Figure 1b). It can be seen that the accuracy decreases with the increase in the value of parameter $\overline{\Phi }$ and with the increase in $c_{\mathrm{A}}^c$. In each case, a very good approximation was obtained for $\overline{\Phi }<2$.

Applying the relationship described in Equation (53) to Equation (13) resulted in:

$$\frac{d{L}_b}{dt}=\frac{1}{3{\rho }_b}\frac{r_b^3-r_0^3}{r_b^2}{r}_{\mathrm{B}}^b(c_{\mathrm{A},av}^b)-\frac{k_{det}{L}_b}{3}\frac{r_b^3-r_0^3}{r_b^2}$$

(54)

with the initial condition:

$$L_b(0)=L_{b0}$$

(55)

The approximate model is a system of ordinary differential equations. The linear version consists of Equation (48) with boundary condition (49) and Equation (54) with boundary condition (55). The nonlinear version consists of Equation (50) with boundary condition (49) and Equation (54) with boundary condition (55). The model that was derived using the Laplace–Carson transform was used in the comparisons that are presented in the further parts of the study and will be abbreviated to AM1.

4. Approximate Model Based on Pseudo-Steady State Approximation

The processes of diffusion and reaction are significantly faster than the growth and detachment of a biofilm, thus the distribution of the limiting substrate concentration in the biofilm is close to a steady-state profile, even when the biofilm thickness is changing [3,8,32].

A stationary profile of the substrate concentration was obtained by solving its mass balance, i.e., Equation (26), after equating the time derivative to zero. It was proposed, in this study, to treat the biofilm as though it were flat. As a result of this, a simpler analytical expression was obtained that was related to the lower number of parameters. In order to evaluate the accuracy of such an approach, the exact biofilm model was solved assuming flat and spherical geometry. The average biomass growth rate, (Equation (14)) was calculated and compared for a range of inert particle radius and biofilm thickness, what is presented in Figure 2. It can be seen that the relative difference between obtained values does not exceed 1.5%.

The dimensionless form of the linearized substrate mass balance in the flat biofilm was expressed as:

$$0=\frac{{\partial }^2\eta }{\partial z^2}-\overline{\Phi }\cdot {\tilde{r}}_{\mathrm{A}}(\eta )$$

(56)

with boundary conditions:

$$\frac{\partial \eta (0)}{\partial z}=0,\eta (1)={\eta }_s$$

(57)

The analytical solution of the linearized form of Equation (56) was:

$$\eta (z)=\frac{\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }-{\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }}+\frac{\cosh \left(z\sqrt{\overline{\Phi }\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }}\right)\left({\tilde{r}}_{\mathrm{A}}({\eta }_s)-\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }+\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }{\eta }_s\right)}{\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }\cosh \left(\sqrt{\overline{\Phi }\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }}\right)}$$

(58)

Equation (59), describing the biofilm dynamics, was derived based on the relationship that connects the overall substrate uptake rate in the biofilm and the substrate mass flux into the biofilm:

$$V^b{r}_{\mathrm{A},av}^b=S^b{D}_{e\mathrm{A}}\frac{d{c}_{\mathrm{A}}^b(L_b)}{dx}$$

(59)

Equation (59) can be transformed to:

$$r_{\mathrm{A},av}^b=\frac{3r_b^2}{(r_b^2-r_0^3)}\frac{D_{e\mathrm{A}}{c}_{\mathrm{A}}^c}{L_b}\frac{d\eta (1)}{dz}$$

(60)

The derivative at the biofilm surface (i.e., for z = 1) was obtained from Equation (58):

$$\frac{d\eta (1)}{dz}=\frac{\tanh \left(\sqrt{\overline{\Phi }\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }}\right)\sqrt{\overline{\Phi }\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }}\left({\tilde{r}}_{\mathrm{A}}({\eta }_s)-\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }+\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }{\eta }_s\right)}{\frac{\partial {\tilde{r}}_{\mathrm{A}}({\eta }_s)}{\partial \eta }}$$

(61)

The average biomass growth rate can be calculated by its relationship with the substrate uptake rate, according to:

$$r_{\mathrm{B},av}^b=w_{\mathrm{BA}}\cdot r_{\mathrm{A},av}^b=w_{\mathrm{BA}}\cdot \frac{3r_b^2}{(r_b^2-r_0^3)}\frac{D_{e\mathrm{A}}{c}_{\mathrm{A}}^c}{L_b}\frac{d\eta (1)}{dz}$$

(62)

This approximate model, including Equations (54), (55), (61) and (62), will be abbreviated, in the further parts of the present study, to AM2-1.

The derivative $\frac{d\eta (1)}{dz}$ can be determined in different ways, which will influence the accuracy of the approximate model. In this study, it was proposed to find this derivative analytically through the finite difference method for one internal node. For the Monod kinetic model, the following algebraic equation was obtained:

$$\frac{{\eta }_0-2{\eta }_{0.5}+{\eta }_s}{{(\Delta z)}^2}-\overline{\Phi }\frac{{\eta }_{0.5}}{\overline{K}+{\eta }_{0.5}}=0$$

(63)

where ${\eta }_0=\frac{4\eta -{\eta }_s}{3}$ and Δz = 0.5.

Equation (63) had two solutions, one of which is positive. It can be expressed as:

$${\eta }_{0.5}=\frac{8-8\overline{K}-3\overline{\Phi }+\sqrt{9{\overline{\Phi }}^2-48\overline{\Phi }+64+48\overline{\Phi }\overline{K}+128\overline{K}+64{\overline{K}}^2}}{16}$$

(64)

The application of the differential quotient of the second order resulted in:

$$\frac{d\eta (1)}{dz}=\frac{8}{3}-\frac{8-8\overline{K}-3\overline{\Phi }+\sqrt{(9{\overline{\Phi }}^2-48\overline{\Phi }+64)+48c_{\mathrm{A}}^c\overline{\Phi }\overline{K}+128\overline{K}+64{\overline{K}}^2}}{6}$$

(65)

The above derivative was used to obtain the average biomass growth rate according to Equation (62). This approximate model, involving Equations (54), (55), (62) and (65), will be abbreviated, in the further parts of the present study, to AM2-2.

Another way to determine the derivative $\frac{d\eta (1)}{dz}$, which is the original aim of this study, is based on a two-step application of the finite difference method. Namely, in the first step, the η_0.5 is determined, as above. In the second step, the same method is used to determine η_0.75, based on the previously found value of η_0.5. Afterwards, the backward difference quotient of the second order can be used to find the derivative, obtaining:

$$\frac{d\eta (1)}{dz}=\frac{32+32\overline{K}-\sqrt{{\mathrm{A}}_2+\overline{K}(48{\mathrm{A}}_1+640\overline{K})+{\mathrm{A}}_3}+\overline{\Phi }}{8}$$

(66)

where:

$${\mathrm{A}}_1=\sqrt{9{\overline{\Phi }}^2-48\overline{\Phi }+64+48\overline{\Phi }\overline{K}+128\overline{K}+64{\overline{K}}^2}$$

$${\mathrm{A}}_2=640-240\overline{\Phi }+25{\overline{\Phi }}^2$$

$${\mathrm{A}}_3={\mathrm{A}}_1(48-8\overline{\Phi })+\overline{K}(1280-16\overline{\Phi })$$

This approximate model, involving Equations (54), (55), (62) and (66), will be abbreviated, in the further parts of the present study, to AM2-3.

A similar derivation can be performed for the Haldane kinetic model. If the procedure that was previously abbreviated to AM2-2 is used, then the following expression is obtained:

$$\frac{d\eta (1)}{dz}=\frac{8}{3}-\frac{8x}{3}$$

(67)

where:

$$x={\left\{q+{\left[q^2+{(r-p^2)}^3\right]}^{1/2}\right\}}^{1/3}+{\left\{q-{\left[q^2+{(r-p^2)}^3\right]}^{1/2}\right\}}^{1/3}+p$$

(68)

$$p=\frac{-b}{3a},q=p^3+\frac{\left(bc-3ad\right)}{6a^2},r=\frac{c}{3a}$$

and:

$$a=1,\hspace{1em}b={\overline{K}}_{in}-1,\hspace{1em}c=3/8\overline{\Phi }{\overline{K}}_{in}+\overline{K}\cdot {\overline{K}}_{in}-{\overline{K}}_{in},\hspace{1em}d=-\overline{K}\cdot {\overline{K}}_{in}.$$

The cubic root $\sqrt[3]{f}$ in Equation (68) denotes the real cube root of f, if f is real. If f is not real, it is then the cube root with the maximum real part.

If the procedure that has been previously referred to as AM2-3 is used, then one gets:

$$\frac{d\eta (1)}{dz}=2{\eta }_{0.5}-8x+6$$

(69)

where: x, p, q and r are defined by Equation (68) and:

$$\begin{gathered} a=1,b=-({\eta }_{0.5}-2\overline{K}+1)/2,c=(32\overline{K}{\overline{K}}_{in}-16{\overline{K}}_{in}-16{\eta }_{0.5}{\overline{K}}_{in}+\overline{\Phi }{\overline{K}}_{in})/32, \\ d=-({\eta }_{0.5}\overline{K}{\overline{K}}_{in})/2-(\overline{K}{\overline{K}}_{in})/2 \end{gathered}$$

The value η_0.5 in the Equation (69) was determined in the first step of the procedure. It was expressed as a real solution of a third-order polynomial, i.e.,

$${\eta }_{0.5}={\left\{q+{\left[q^2+{(r-p^2)}^3\right]}^{1/2}\right\}}^{1/3}+{\left\{q-{\left[q^2+{(r-p^2)}^3\right]}^{1/2}\right\}}^{1/3}+p$$

(70)

where: p, q and r are defined by Equation (68) and:

$$a=1,b={\overline{K}}_{in}-1,c=\frac{3}{8}\overline{\Phi }{\overline{K}}_{in}+\overline{K}\cdot {\overline{K}}_{in}-{\overline{K}}_{in},d=-\overline{K}\cdot {\overline{K}}_{in}.$$

5. Accuracy and Efficiency Evaluation of the Approximate Model Based on the Laplace-Carson Transform

Calculations were performed for two of the microbiological processes that are of large technological significance: the oxidation of nitrites by nitrifying bacteria and the phenol degradation that is performed by Pseudomonas putida. The former follows Monod kinetics (Equation (23)), whereas the latter follows Haldane kinetics (Equation (24)). The kinetic parameters are presented in

Table 1. Values of kinetic parameters for processes following Monod kinetic model [15].

Process	k [1/h]	KA [kg/m3]	wBA [kg B/kg A]
Nitrite oxidation	7.917 × 10−3	2.4 × 10−4	0.23

and

Table 2. Values of kinetic parameters for biodegradation of phenol following Haldane kinetic model [33].

k [1/h]	KA [kg/m3]	Kin [kg/m3]	wBA [kg B/kg A]
0.26	0.0254	0.173	0.616

, respectively. A stepwise change of the substrate concentration at the biofilm surface $c_{\mathrm{A}s}^b=c_{\mathrm{A}}^c=0.1{\mathrm{kg}/\mathrm{m}}^3$ and a constant biofilm thickness were used for the purpose of the evaluation. The method of lines was used for solving the exact model (see Appendix B).

Figure 3 presents the relationship, η_av(τ), that was obtained according to the exact model and the approximate models known as AM1 (i.e., the linear one—Equations (48) and (49)) and the nonlinear one—Equations (50) and (49). Figure 3a illustrates the results that were obtained for the oxidation of nitrites, while Figure 3b shows those that were obtained for the phenol degradation. The results that were obtained indicate that:

• The accuracy of the approximate models decreased with the increase in $\overline{\Phi }$;
• The nonlinear approximate model was more accurate than the linear one;
• The differences between the approximate models and the accurate model decreased with time.

The above observations are in agreement with those presented by Szukiewicz [10].

The proposed approximate mathematical models, known as AM1, were evaluated with regard to their computational efficiency. Matlab 2016a was used for this purpose. The execution times of the computer programs using the exact model and the approximate models were compared. The exact model was solved, as previously explained, by the method of lines (MOL). Two integration algorithms were used in all of the programs. The first one was the Dormand–Prince (4/5) algorithm, which is an explicit method that is used for solving non-stiff and moderately stiff differential equations [34]. The second one was Gear’s algorithm, which is commonly used for stiff problems [26,35,36]. The results of the comparison are shown in

Table 3. Execution times of computer programs solving the exact model and the approximate models.

Model	Exact Model	Linear Approximate Model	Nonlinear Approximate Model
Integration algorithm	Dormand-Prince (4/5) (MOL)	Dormand-Prince (4/5)
Normalized time of execution	1	0.0468	0.0327
Integration algorithm	Gear’s method (MOL)	Gear’s method
Normalized time of execution	0.0816	0.0561	0.0535

. The values of the execution times were normalized through the process of dividing by the value that was obtained for the exact model and the explicit integration algorithm. The application of the implicit algorithm, i.e., Gear’s algorithm, significantly increased the performance of the exact model. This suggests that the system of equations that were obtained by the method of lines was stiff. In order to verify this hypothesis, the stiffness ratio was calculated, according to the formula:

$$SR=\frac{\max \left[\mathrm{Re}\left|{\lambda }_i\right|\right]}{\min \left[\mathrm{Re}\left|{\lambda }_i\right|\right]}$$

(71)

In Equation (71), λ_i represents the eigenvalues of the Jacobi matrix of the right hand sides of the system of the equations that was considered. If the stiffness ratio (SR) is less than 100, then the system of differential equations is not regarded to be stiff [37]. It can be seen, in Figure 4, that the stiffness ratio was over 1000 for the whole range of dimensionless time.

It can be seen, from Table 3, that when the explicit algorithm was used for the integration, the computations that used the approximate models were about 20–30 times faster than those that used the exact one. This makes it advantageous for engineering purposes.

6. Application of Approximate Models to Conditions of the Dynamical Growth of Biofilm

Figure 5 presents the results of the dynamical simulations of the biofilm growth that were obtained with the use of the exact model and the approximate models AM1 for three different values of the substrate concentration in bulk liquid. Figure 5a,b illustrate the average substrate concentration in the biofilm and biofilm thickness that were obtained for the process of nitrite oxidation that was performed by the nitrifying bacteria. It can be seen that very low values of biofilm thicknesses were obtained, a result which arises from the low value of the maximum specific growth rate of this process. Due to low internal mass transfer resistances, excellent accuracy of the approximate model was obtained. Figure 5c,d present the results that were obtained from the glucose utilization by Pseudomonas aeruginosa bacteria (the kinetic parameters of which are shown in

Table 4. Values of kinetic parameters for glucose utilization following Monod kinetic model [38].

k [1/h]	KA [kg/m3]	wBA [kg B/kg A]
0.4	0.002	0.34

). This process follows Monod kinetics, so too does the nitrite oxidation process; however, the maximum specific growth rate was over 50 times greater. Significantly larger biofilm thicknesses were obtained for this than for the nitrite oxidation. The values of parameter $\overline{\Phi }$ at a steady-state condition (i.e., for a sufficiently long time) equal $\overline{\Phi }=70.1$ for $c_{\mathrm{A}}^c=0.01{\mathrm{kg}/\mathrm{m}}^3$, $\overline{\Phi }=56.5$ for $c_{\mathrm{A}}^c=0.02{\mathrm{kg}/\mathrm{m}}^3$ and $\overline{\Phi }=43.7$ for $c_{\mathrm{A}}^c=0.04{\mathrm{kg}/\mathrm{m}}^3$. Due to the very large internal mass transfer resistances, the accuracy of the approximate model AM1 was significantly lower than that of the nitrite oxidation. The values of $\overline{\Phi }$ considerably exceeded the ones for which the distributions of the dimensionless substrate concentration in biofilm were obtained with good accuracy (see Figure 3). Figure 5d,e show the results that were obtained for phenol degradation by Pseudomonas putida. In this process, good accuracy was obtained for the smaller values of the substrate concentration in the liquid. The values of parameter $\overline{\Phi }$ at a steady-state condition were significantly lower than those of the glucose utilization process, equaling $\overline{\Phi }=9.76$ for $c_{\mathrm{A}}^c=0.01{\mathrm{kg}/\mathrm{m}}^3$, $\overline{\Phi }=9.58$ for $c_{\mathrm{A}}^c=0.02{\mathrm{kg}/\mathrm{m}}^3$ and $\overline{\Phi }=8.40$ for $c_{\mathrm{A}}^c=0.04{\mathrm{kg}/\mathrm{m}}^3$. The lower values of this parameter were related to the biokinetic parameters. The phenol degradation was characterized by a smaller value of the specific growth rate and a larger value of the growth yield coefficient than glucose utilization.

Figure 6 presents the results that were obtained with the use of the AM2 models and the results that were obtained with the use of MOL, for comparison. The glucose utilization (Figure 6a–c) and phenol degradation (Figure 6d–f), with different values of substrate concentration in the liquid, were taken into consideration. It can be seen that the highest accuracy for both processes was obtained with the AM2-3 model, which makes it suitable for modeling microbiological processes following various kinetic models.

Abbas and Eberl [39] derived an analytical approximate model of a microbiological process in a biofilm by using an original method that was based on the homotopy analysis method. The obtained substrate profile in the biofilm can be expressed as:

$$c_{\mathrm{A}}^b(z)=c_{\mathrm{A}}^c+\frac{{\rho }_{b}k}{w_{\mathrm{BA}}{D}_{e\mathrm{A}}}\frac{L_b^2(z^2-1)}{2}\frac{c_{\mathrm{A}}^c}{K+c_{\mathrm{A}}^c}$$

(72)

The above model was used afterwards by the quoted authors in order to determine the diffusive flux into a biofilm and to obtain the equation of the biofilm dynamics. It can be presented in the following way:

$$\frac{d{L}_b}{dt}=k\cdot \frac{c_{\mathrm{A}}^c{L}_b}{K+c_{\mathrm{A}}^c}-k_{det}{L}_b^2$$

(73)

The disadvantage of the above model is that it can only be used for Monod kinetics and for flat biofilm geometry. The simulations of the biofilm growth in the present study were performed using the method of lines, model of Abbas and Eberl and AM2-3. The dynamical changes of the biofilm thickness are presented in Figure 7. The low accuracy of the model that was formulated by Abbas and Eberl was caused by the large internal mass transfer resistances that are present in the biofilm. It was shown in a previous study [13] that solutions that are obtained using the homotopy analysis method require a very large number of terms in the approximate series in order to achieve acceptable accuracy. The approximate series in the solution that was formulated by Abbas and Eberl takes into account only two terms. It can be seen that a significantly higher level of accuracy can be obtained by using the model AM2-3 that is proposed in this study.

7. Conclusions

This paper presents original approximate models for the dynamics of the microbiological processes that occur in a biofilm that is formed on fine spherical particles. The model that was based on the Laplace–Carson transform can be used for the modeling of processes following any single-substrate kinetics. Good accuracy of the model was obtained when it was used for low and internal mass transfer resistances. However, it could not be applied to thick biofilms that are characterized by large internal mass transfer resistances. It was shown that, for such conditions, significantly higher accuracy was obtained with the use of the proposed models, based on the concept of the pseudo-stationary substrate concentration profile in the biofilm.

The advantage of the dynamical biofilm growth model that is based on the pseudo-stationary concept is that it can be modified by using any approximate analytical solution of stationary substrate concentration profile. The approach that has been presented in this study was compared with a previously published dynamical biofilm model. It was shown that significantly higher accuracy was obtained by using the proposed model.

The application of the approximate models simplifies the computations significantly and eliminates the stiffness from the exact model, solved by the method of lines. These models can be used for the modeling and analysis of microbiological processes and for the design of automatic control systems.