The Modified Trickling Filter: A Mathematical Model

Abstract

The primary objective of this study is to develop a modified trickling filter (MTF) mathematical model that addresses limitations of conventional biofilm models, including their inability to capture microbial stratification, simultaneous nitrification–denitrification, and diffusion effects in porous sponge media. Under fixed-bed circumstances, the MTF uses a polyurethane sponge medium to promote biomass adhesion and proliferation. The two main zones of the MTF are aerobic and anoxic. Heterotrophic organic carbon oxidation and nitrification are both carried out concurrently in the aerobic zone. The two main zones of the MTF are aerobic and anoxic: heterotrophic organic carbon oxidation and nitrification occur in the aerobic zone, whereas denitrification proceeds in the anoxic zone. The model integrates three components: a compartmental reactor flow model (axial plug-flow representation of completely mixed biofilm reactors), a biofilm kinetics module (substrate conversion by heterotrophs, nitrifiers, and denitrifiers), and a porous-media diffusion module (effective diffusion via the Millington–Quirk relationship). Within the biofilm, both bilayer distributions (heterotrophs outer, nitrifiers inner) and homogeneous distributions were evaluated. Model predictions were validated against experimental data (influent COD = 119–161 mg L⁻¹, NH₃-N = 29–35 mg L⁻¹, hydraulic residence times = 2–4 h, with/without recycle) using parity charts and statistical indices. The bilayer biofilm with pore diffusion provided the most accurate results, with Willmott’s index of agreement (d) = 0.95 for COD, 0.95 for Ammoniacal Nitrogen (NH3-N), and 0.94 for Nitrate Nitrogen (NO₃-N) and relative errors (RE) of 0.05–0.07, compared with d ≤ 0.47 and RE ≥ 0.24 for the homogeneous model. These findings confirm that incorporating stratification and pore diffusion significantly improves model performance. The framework not only clarifies COD, NH₃, and NO₃ removal mechanisms in sponge-based MTFs but also offers design guidance on recycle ratios and porosity thresholds for enhanced operation.

Graphical Abstract

Keywords

diffusion-reaction kinetics mathematical modeling multispecies biofilm porous media trickling filter reactor

Introduction

Wastewater treatment has been a global concern for decades. The first known biological wastewater treatment systems utilized beds filled with rocks in 1901; however, the mechanism involved was unknown for many decades of the 20th century (Peters and Alleman, 1982). The concept of trickling filters that commonly utilize pebbles and Raschig rings originated from this ancient way of water filtration. The simple mechanism of biomass immobilization on a solid support, involved in the trickling filter, makes it even more easily applicable. With the invention of various microscopic techniques, researchers became able to visualize the slime layer coated on the rock surface, known as biofilm (Bakke and Olsson, 1986; Chatterjee et al., 2014; Cooke, 1959; Costerton et al., 1978; Palmer and Sternberg, 1999; Seymour et al., 2004). Many scientists worked on improving the efficiency of trickling filters; however, the mechanism of pollutant removal was still unknown. Two researchers named Germain and Schulze established empirical correlations for the treatment efficiencies of trickling filters (Germain, 1966; Schulze, 1960). These correlations helped determine the experimental efficiencies of typical trickling filters for limited conditions, but still could not explain the mechanism of nutrient removal.

The mathematical model can offer a far more comprehensive approach to system analysis. Rather than doing tests under particular settings, researchers may examine all treatment-related elements and processes using a mathematical model as well as how changing one parameter might impact the entire system. The 1970s saw the development of fundamental analytical techniques by several authors who used mathematical modeling in wastewater treatment. These techniques assume that the biomass in the biofilm is equally distributed across the biofilm matrix and that the biofilm is a continuum. For suspended culture systems, the simple Monod growth kinetics can be the most effective (Alsaeed et al., 2023; Nuhoglu et al., 2005), but it would not work for biofilm systems. To simulate a biofilm reactor, one of the main goals is to examine biofilm kinetics because the transport is diffusion limited. The mass movement of nutrients across the biofilm matrix via microbial breakdown and diffusion was taken into account by the fundamental biofilm model (Costerton et al., 1978; Rittmann and McCarty, 1980). Diffusion and microbial reactions with Monod kinetics are included in the biofilm model. The biggest obstacle is that the Monod growth kinetics is a variable order, and there is no analytical solution to that equation (Rittmann and McCarty, 1978). These analytical limitations of classical Monod-based biofilm formulations underscore the need for advanced modeling frameworks that explicitly account for diffusion limitations and biofilm structure.

The analytical solution of the model exists only under two circumstances. Either it can be presumed to be the first order when the concentration of the substrate is less than the Monod half-saturation constant, or Monod kinetics simply decreases to zeroth order when the substrate concentration is significantly higher than the half-saturation constant (Haug and Mccarty, 1972; Mulcahy and La Motta, 1978). Since substrate concentrations in biofilm kinetics are not always in high or low ranges and both situations are severe, two methods might be used: half-order kinetics, which involves averaging the analytical solutions of first-order and zeroth-order kinetics (Harris and Hansford, 1976; Rittmann and McCarty, 2001), which is an erroneous approximation. The second method, on the other hand, is numerical and uses an appropriate software program to solve the model with a minimal margin of error (Kissel et al., 1984; Versteeg and Malalasekera, 2007). The solution of the biofilm model is necessary for the biofilm reactor modeling. To determine the substrate mass flow into the biofilm surface, the biofilm model solution is utilized. The model of the reactor is developed using the mass flux. The biofilm model was seen as a uniform layer of microorganisms in the 1970s and early 1980s. These models may be used to investigate a number of phenomena, such as mass transport, molecular diffusion, and substrate utilization kinetics. With a single species and a single limiting substrate, these models presupposed a steady-state methodology (Jennings et al., 1976; Williamson and McCarty, 1976). The idea of steady-state biofilm thickness is presented by Rittmann and McCarty. In their study of substrate usage to biofilm growth, they found that biofilms maintain a consistent thickness, have no net growth or decay, and have identical rates of biomass growth and decay (Williamson and McCarty, 1976).

The biofilm modeling method used in previous studies worked well for assessment systems that have a single dominating process (such as COD removal or nitrification). These techniques were not very effective in predicting the precise distributions of various microorganisms, the conversion of many substrates, the overall buildup of biofilms, or the intricate structure of biofilms. Using a continuum method, several authors provided a more explicit depiction of a multispecies dynamic biofilm model. The model forecasts the impact of variations in biofilm thickness (Wanner and Gujer, 1986). The basic model was expanded into a more complex representation around 10 years later, representing biomass as particulate components and an extracellular matrix as a liquid phase. The model also incorporated impacts on the liquid fraction of the biofilm created by the displacement of particulates within the matrix. The biomass was considered to be fixed in the biofilm in the previous model. Particulate flow was also taken into account by the model (Wanner and Reichert, 1996). A concise timeline summarizing the evolution from conventional trickling filters to sponge-based modified trickling filters (MTFs) is provided in Supplementary Figure S1. The biofilm reactor may be modeled with the help of the recently established biofilm kinetics model, which is based on a dimensionless quantity called “biofilm number B” (Sharma et al., 2022).

Sponge carrier media are widely reported to enhance biomass retention and improve nitrification and denitrification compared with conventional rock or plastic media, primarily due to their high porosity and large specific surface area (Deguchi and Kashiwaya, 1994; Mattei et al., 2018; Morgenroth, 2020; So et al., 2014). The extent of improvement varies with operating conditions such as loading, recycle ratio, and sponge porosity, whereas excessive pore clogging can reduce permeability and slow reaction rates (So et al., 2014). Recent studies further confirm the relevance of sponge-based carriers: Jiang et al. (2021) demonstrated improved pollutant removal using modified polyurethane (PU) sponge fillers in moving bed biofilm reactors (Jiang et al., 2021); Wang et al. (2022) showed that PU sponge carriers enhance substrate diffusion and removal efficiency in biofilm electrode reactors (Wang et al., 2022); and Chalipa et al. (2024) reported effective pharmaceutical removal in sponge-based MBBRs (Chalipa et al., 2024). These findings highlight that sponge morphology, porosity, and maintenance practices (e.g., backwashing to prevent clogging) remain critical design variables in current applications.

Building on these insights, this study develops and validates a mechanistic multispecies model of the MTF that, unlike prior biofilm models, explicitly incorporates microbial stratification (homogeneous vs. bilayer structures) and pore-scale diffusion in PU sponge carriers. Adapted from Sharma et al. (2022), the model integrates heterotrophic and nitrifying kinetics and was tested against nine experimental datasets (influent COD = 119–161 mg L⁻¹, NH₃-N = 29–35 mg L⁻¹, hydraulic residence time [HRT] = 2–4 h, with and without recycle). Validation showed that the bilayer model with pore diffusion achieved Willmott’s index of agreement (d) ≈ 0.95 and relative error (RE) ≤ 0.07 for COD, NH₃-N, and NO₃-N removal compared with d ≤ 0.47 and RE ≥ 0.24 for the homogeneous model.

Materials and Methods

The MTF’ description

The two zones that make up the MTF are the anoxic zone at the bottom and the aerobic zone at the top (Fig. 1). After entering the aerobic zone, the influent wastewater descends into the anoxic zone. As a carrier material for biomass immobilization in both zones, the PU sponge is used. Anoxic conditions are maintained in the bottom zone, whereas holes are positioned on the aerobic zone’s edge for natural aeration. In the aerobic zone, heterotrophs and nitrifiers use oxygen to oxidize substrates of organic carbon and ammonia, respectively. For heterotrophic denitrification, the oxidized products are transferred to an anoxic zone.

FIG. 1.

Schematic diagram of a modified trickling filter.

Biofilm processes

The aerobic zone consists of a biofilm containing two species: (1) nitrifiers and 2) heterotrophs. To flourish, the heterotrophs need oxygen and organic carbon. Heterotrophic metabolism is represented by the following equation:

C_{6} H_{12} O_{6} + 6 O_{2} \to 6 C O_{2} + 6 H_{2} O

This oxidation process provides energy for the growth of bacteria. In conjunction with heterotrophs, the nitrifiers also utilize oxygen for their growth. In two steps, the two nitrifier species convert ammonia substrate into nitrate: (1) Ammonia is converted to nitrite by nitrosomonas, nitrosospira, and nitrosococcus species (Burrell et al., 1999), and (2) nitrite is oxidized into nitrate by nitrobacter; the reaction is as follows:

N H_{3} + 1.5 O_{2} \to N O_{2}^{-} + H^{+} + H_{2} O (Nitrosomonas)

N O_{2}^{-} + 0.5 O_{2} \to N O_{3}^{-} (Nitrobacter)

Both species, as autotrophic nitrifiers in wastewater treatment systems, can be categorized as a single, wide species (Horan, 2004; Wanner et al., 2006). Biofilm containing facultative heterotrophs makes up the anoxic zone, which receives the wastewater from the aerobic zone. The anoxic zone lacks oxygen; therefore, the heterotrophs change their metabolism to use nitrates as terminal electron acceptors. This results in heterotrophic denitrification. The reaction is shown as follows:

{5 C}_{6} H_{12} O_{6} (COD) + 24 N O_{3}^{-} + 24 H^{+} \to 30 C O_{2} + 42 H_{2} O + 12 N_{2}

Denitrification reduces oxidized nitrogenous compounds to nitrogen gas (Payne, 1981).

Model development

To develop the mathematical model for the MTF, the following elements are needed: (1) Biofilm model, (2) reactor model, and (3) reactor flow model and are covered in this section.

The following presumptions were used for modeling: 1.

The biomass within the biofilm is treated as a continuum (Sharma et al., 2022).

The substrate has negligible mass-transfer resistance in the bulk liquid phase.

Substrate movement inside the biofilm is described by Fick’s law of diffusion (Wanner et al., 2006).

The diffusivity within the biofilm is assumed to be 80% of the diffusivity of that substrate in water (Morgenroth, 2020; Rittmann and McCarty, 2001; Wanner et al., 2006).

The reaction within the biofilm follows Monod growth kinetics (Wanner et al., 2006).

Steady-state conditions exist everywhere. The basis of steady-state biofilm is taken from the literature (Sharma et al., 2022).

Only diffusive transport and microbial response inside the biofilm are taken into account (Wood and Whitaker, 1998).

All components are assumed to be transported normally to the substratum.

The heterotrophic and nitrifying metabolism are substrate limited, and oxygen is assumed to be abundant.

10.

The biofilm has a uniform biomass density.

Biofilm model

Two species make up the biofilm in the aerobic zone: nitrifiers and heterotrophs, as discussed in the “Biofilm Processes” section. Because both groups depend on oxygen, there is potential competition for oxygen consumption. For simplicity, oxygen was assumed to be abundant in the biofilm. Given their distinct growth rates (heterotrophs as fast growers, nitrifiers as slow growers), heterotrophs are typically found near the outer layer and nitrifiers deeper within the biofilm. In this study, bilayer and homogeneous arrangements were prescribed as alternative configurations rather than emergent predictions, to explicitly test their impact on model outcomes. The bilayer structure represents stratification supported by earlier experimental and modeling studies (Rittmann and Manem, 1992; Wanner and Gujer, 1985; Wanner and Reichert, 1996; Watanabe et al., 1995), whereas the homogeneous case was included as a control scenario where heterotrophs and nitrifiers are evenly distributed across the biofilm volume. Click or tap here to enter text. For the anoxic zone, it was assumed that only facultative heterotrophs were present, following facultative aerobic metabolism. This zone was therefore modeled as a single-species biofilm.

The biofilm model’s derivation

To determine the substrate consumption kinetics inside the differential section of a biofilm adhering to the carrier, the biofilm model is derived using the microscopic mass balance. This section covers the derivation of the two types of aerobic biofilm models (1) bilayer biofilm model [Fig. 2(a)], (2) homogeneous biofilm [Fig. 2(c)], and (3) denitrifying biofilm.

FIG. 2.

(a) Biofilm bilayer model (green = nitrifiers, blue = heterotrophs). (b) Heterotroph and nitrifier boundary conditions for the bilayer model. (c) A homogeneous biofilm model, where green represents nitrifiers and blue represents heterotrophs. (d) Boundary conditions for the bilayer model for heterotrophs and nitrifiers.

Bilayer biofilm model

General steady-state component mass balance within a biofilm can be written as follows:

(Rate o f component i n) - (Rate o f component out) = (Accumulation) + (Consumption) - (Production)

(2.1)

At the steady state, accumulation = 0. The mass balance equation’s component degradation and production terms describe the kinetics of the corresponding component-degrading/producing microbial group.

Consider a layer of biofilm thickness Δx within the biofilm, applying mass balance for a species whose mass flux J_F moving across the Δx element is written as:

A J_{F} |_{x} - A J_{F} |_{x + Δ x} - Δ x A r_{F} = 0

(2.2)

Applying Taylor’s series, for $J_{F} A |_{x + △ x}$ can be written as follows:

J_{F} {A |}_{x + △ x} = J_{F} {A |}_{x} + △ x \frac{d J_{F}}{d x} A

(2.3)

Substituting (2.3) in (2.2),

J_{F} {A |}_{x} - J_{F} {A |}_{x} - △ x \frac{d J_{F}}{d x} A + r_{F} A Δ x = 0

(2.4)

- \frac{d J_{F}}{d x} + r_{F} = 0

(2.5)

Rearrangement of the equation, as r_F is a degrading term,

According to Monod kinetics $r_{F} = - \frac{μ_{\max}}{Y} (\frac{S}{S + K_{S}}) X_{F}$ and Fick’s law $J_{F} = - D_{F} \frac{d S}{d x}$

D_{F} \frac{d^{2} S}{d x^{2}} - \frac{μ_{\max}}{Y} (\frac{S}{S + K_{S}}) X_{F} = 0

(2.6)

For bilayer biofilm Equation (2.6) is as follows:

Heterotrophs

D_{F, H} \frac{d^{2} S_{H}}{d x^{2}} - \frac{μ_{\max, H}}{Y_{H}} (\frac{S_{H}}{S_{H} + K_{S, H}}) X_{F} = 0

(2.7)

Rearranging the equation into dimensionless form by substituting $S_{H}^{*} = \frac{S_{H}}{K_{S, H}}$ and multiplying and dividing by K_{S, H} on both sides, we get

\frac{d^{2} S_{H}^{*}}{d {x^{*}}^{2}} - B_{H}^{2} (\frac{S_{H}^{*}}{S_{H}^{*} + 1}) = 0

(2.8)

Where $B_{H}^{2} = \frac{μ_{\max, H}}{Y_{H} D_{F, H} K_{S, H}} L_{H}$ , and taking the square root of $B_{H}^{2}$ , we get the biofilm number of heterotrophs within the biofilm

B_{H} = L_{H} \sqrt{\frac{μ_{\max, H} X_{f, H}}{Y D_{H} K_{S, H}}}

(2.9)

The boundary conditions [refer Fig. 2(b)] for the heterotrophic layer are

BC.1 : A t x^{*} = 0 (Biofilm surface), S_{H} = S_{H, surf}^{*}, J_{F}^{*} = J_{F, H surf}^{*}

(2.10)

BC 2 : A t x^{*} = \frac{x_{1}}{L_{H}} (Heterotroph nitrifier interface), S_{H}^{*} = S_{H, wall}^{*}, J_{F}^{*} = 0 (N o flux conditions)

(2.11)

Note: The inner layer, which primarily consists of the nitrifiers, is thought to be completely penetrated by the COD substrate, leaving the heterotroph–nitrifier interface flux free.

Similarly, for nitrifiers:

\frac{d^{2} S_{N}^{*}}{d {x^{*}}^{2}} - B_{N}^{2} (\frac{S_{N}^{*}}{S_{N}^{*} + 1}) = 0

(2.12)

Like Equation (2.8), the biofilm number for the nitrifiers within the biofilm can be written as

B_{N} = L_{N} \sqrt{\frac{μ_{\max, N} X_{f, H}}{Y D_{N} K_{S, N}}}

(2.13)

The boundary conditions [refer Fig. 2(b)] for the nitrifiers are

BC.1 : A t x^{*} = \frac{x_{1}}{L_{H}} (Heterotroph nitrifier interface), S_{N}^{*} = S_{N, surf}^{*}, J_{F}^{*} = J_{N, surf}^{*}

(2.14)

Note: Similarly, the heterotrophic layer has been entirely penetrated for ammonia substrate (SN*), which leads to substrate consumption and provides no hindrance to mass transfer.

BC 2 : A t x^{*} = \frac{x_{2} - x_{1}}{L_{N}} (Substratum), S^{*} = S_{wall}^{*}, J_{F, N}^{*} = 0 (N o flux c onditions)

(2.15)

From the previous work (Sharma et al., 2022), the derived expression for the dimensionless flux for the single rate-limiting substrate and single-species model is given by

J_{F}^{*} = B^{2} \int_{0}^{1} (\frac{S^{*}}{S^{*} + 1}) d x^{*} = B^{2} Ψ (S^{*})

(2.16)

The expression for dimensionless flux is derived from the Equations (2.5) and (2.15) as follows:

For heterotrophs

J_{F, H}^{*} = B_{H}^{2} \int_{0}^{\frac{x_{1}}{L_{H}}} (\frac{S_{H}^{*}}{S_{H}^{*} + 1}) d x^{*}

(2.17)

J_{F, H}^{*} = B_{H}^{2} Ψ (S_{H}^{*})

(2.18)

Where $Ψ (S_{H}^{*}) = \int_{0}^{\frac{x_{1}}{L_{H}}} (\frac{S_{H}^{*}}{S_{H}^{*} + 1}) d x^{*}$

Similarly, for nitrifiers

J_{F, N}^{*} = B_{N}^{2} \int_{\frac{x_{1}}{L_{H}}}^{\frac{x_{2} - x_{1}}{L_{N}}} (\frac{S_{N}^{*}}{S_{N}^{*} + 1}) d x^{*}

(2.19)

J_{F, N}^{*} = B_{N}^{2} Ψ (S_{N}^{*})

(2.20)

Where $Ψ (S_{H}^{*}) = \int_{0}^{\frac{x_{1}}{L_{H}}} (\frac{S_{H}^{*}}{S_{H}^{*} + 1}) d x^{*}$ , $Ψ (S_{N}^{*}) = \int_{\frac{x_{1}}{L_{H}}}^{\frac{x_{2 - x_{1}}}{L_{N}}} (\frac{S_{H}^{*}}{S_{H}^{*} + 1}) d x^{*} .$ These integrals can be quickly evaluated from the expressions derived in the work by Sharma et al. (2022). The integral limits can be reduced for both layers to 0 and 1 by replacing $x_{1} = L_{H}$ so that $\frac{x_{1}}{L_{H}} = 1$ . For nitrifiers, if we take the reference limit 0, then $x_{2} - x_{1} = L_{N}$ so that $\frac{x_{2} - x_{1}}{L_{N}} = 1$ .

Homogeneous biofilm model

Given that the nitrifiers and heterotrophs have volume fractions f_H and f_N [refer Fig. 2(c)], respectively, the biomass density of the two groups inside the biofilm is determined by

X_{H} = X_{F} f_{H}

(2.21)

X_{N} = X_{F} f_{N}

(2.22)

Heterotroph model equations can be derived from Equation (2.8) as follows:

\frac{d^{2} S_{H}^{*}}{d {x^{*}}^{2}} = B_{H}^{2} (\frac{S_{H}^{*}}{S_{H}^{*} + 1})

(2.23)

In the same way, for nitrifiers

\frac{d^{2} S_{N}^{*}}{d {x^{*}}^{2}} = B_{N}^{2} (\frac{S_{N}^{*}}{S_{N}^{*} + 1})

(2.24)

The biofilm numbers for heterotrophs and nitrifiers can be formulated for the homogeneous model by substituting Equations (2.21) and (2.22) as follows:

B_{H} = L_{H} \sqrt{\frac{μ_{\max, H} X_{F} f_{H}}{Y_{H} D_{H} K_{S, H}}}

(2.25)

B_{N} = L_{N} \sqrt{\frac{μ_{\max, N} X_{F} f_{N}}{Y D_{N} K_{S, N}}}

(2.26)

The expression for the dimensionless mass flux for heterotrophs and nitrifiers from Equation (2.16) can be written as

J_{H}^{*} = B_{H}^{2} Ψ (S_{H}^{*})

(2.27)

J_{N}^{*} = B_{N}^{2} Ψ (S_{N}^{*})

(2.28)

Heterotrophs [refer to Fig. 2(d)]:

BC.1 : A t x^{*} = 0 (Biofilm surface), S_{H}^{*} = S_{, surf}^{*}, J_{F H}^{*} = J_{H, surf}^{*}

(2.29)

BC 2 : A t x^{*} = 1, S_{H}^{*} = S_{wall}^{*}, J_{F, H}^{*} = 0 (N o flux conditions a t substratum)

(2.30)

Nitrifiers [refer to Fig. 2(d)]:

BC 1 : A t x^{*} = 0 (Biofilm surface), S_{N}^{*} = S_{, surf}^{*}, J_{F, N}^{*} = J_{N, surf}^{*}

(2.31)

BC 2 : A t x^{*} = 1, S_{H}^{*} = S_{wall}^{*}, J_{F, N}^{*} = 0 (N o flux conditions a t substratum)

(2.32)

Biofilm model for denitrification

From Equation (2.8), the model equation for denitrifiers, using stoichiometric parameters in Table 1, can be written for nitrate substrate as follows:

\frac{d^{2} S_{N O_{3}}^{*}}{d {x^{*}}^{2}} = B_{N O_{3}}^{2} (\frac{S_{N O_{3}}^{*}}{S_{N O_{3}}^{*} + 1})

(2.33)

B_{{N O}_{3}} = L_{f} \sqrt{\frac{{(1 - Y_{H}) μ}_{\max, H} X_{H}}{2.86 Y_{H} D_{N O_{3}} K_{N O_{3}}} η_{H}}

(2.34)

Table 1.

Stoichiometric Table for the Microbial Processes

Processes	S_S	S_N	$S_{N O_{3}}$	Process rate
Denitrification metabolism	$- \frac{1}{Y_{H}}$		$- \frac{(1 - Y_{H})}{2.86 Y_{H}}$	$μ_{\max, H} η_{H} (\frac{S_{F, H}}{S_{F, H} + K_{S, H}}) X_{F, H}$
Autotrophs metabolism		$\frac{1}{Y_{N}}$		$μ_{\max, N} \frac{S_{F, N}}{S_{F, N} + K_{S, N}} X_{F, N}$
Heterotrophic substrate removal	$- \frac{1}{Y_{H}}$			$μ_{\max, H} \frac{S_{F, H}}{S_{F, H} + K_{S, H}} X_{F, H}$

Heterotrophs (H), nitrifiers (N), and anoxic organic carbon oxidation efficiency (η).

Stoichiometrically, 1 gram of nitrate oxidizes 2.86 grams of COD substrate

J_{N O_{3}}^{*} = B_{N O_{3}}^{2} Ψ (S_{N O_{3}}^{*})

(2.35)

Model for porous media

The PU sponge media used in the MTF is different from regular Raschig rings and other nonporous media. Biofilm in Raschig rings or rocks may be consistent in thickness; however, in the pore network of a PU sponge, biofilm may vary in architecture depending on the microscopic pore structure. The hydrodynamic conditions and pore structure have a major impact on pore diffusion and biomass morphology (Rittmann, 1993). The attachment mechanism in the porous medium of biomass differs from that of the smooth surface. The effective diffusion coefficient will take the place of the standard diffusion coefficient if pore diffusion is taken into account in the system. The definition of the effective diffusion coefficient is

D_{eff, i} = \frac{D_{f} E}{τ}

(2.36)

where D _f is the molecular diffusion coefficient of species i in the bulk fluid phase, E is the porosity of the medium, and τ is the tortuosity of the pore channel inside the porous media. The Millington and Quirk Model explains the connection between porosity and tortuosity (Chou et al., 2012).

τ = E^{- 0.333}

(2.37)

While the Millington–Quirk formulation is widely applied for soils and granular porous media, its direct application to sponge carriers should be interpreted with caution. Sponge matrices differ structurally, exhibiting interconnected pores, higher tortuosity, and variable compression under flow, which may affect effective diffusivity estimates (Jiang et al., 2021; Mattei et al., 2018; So et al., 2014). Here, it is used as a first-order approximation, consistent with prior biofilm modeling studies.

This representation assumes constant porosity and tortuosity for the sponge medium. Although pore clogging and occlusion are known to reduce effective diffusivity and permeability in practice (So et al., 2014), these dynamic effects were not explicitly modeled here. Instead, the present model uses fixed E and τ values derived from literature and experimental measurements. Observed trends from prior studies suggest that excessive pore clogging can lower effective diffusivity and thereby reduce reaction rates, but in this framework, such effects are only acknowledged qualitatively and not mechanistically simulated.

Reactor model

The global mass balance on a differential axial portion of the reactor column is used to create the reactor model. (1) The compartment model and (2) the flow model comprise the reactor model. Each differential section or compartment of a trickling filter can be considered a completely mixed biofilm reactor (CMBR) (Rittmann and McCarty, 2001; Wik, 2003). CMBRs are assumed to be connected in series in this flow model, resulting into a plug-flow reactor (PFR) design.

Compartment model: CMBR

The CMBR could potentially be considered as a typical biofilm reactor. Bulk flows through a CMBR after the biofilm has been completely mixed in one compartment. The whole resultant would behave like a chemostat.

Flow model: CMBRs in series

The trickling filter modeling assumes that the CMBRs are connected in series. CMBRs connected in succession would produce a PFR overall. The uniform distribution of the influent in a deceiving filter eliminates the assumption of radial flux fluctuation. The rate of reaction or flux will fluctuate in the reactor’s axial direction due to the plug-flow nature (Rittmann and McCarty, 2001; Wik, 2003).

Derivation of the reactor model

The reactor model is derived by applying the steady-state mass balance axially to the differential reactor volume element ΔV (Fig. 3a). This element assumes that the substrate’s conversion is ΔX, as shown below:

(Rate o f mass o f substrate entering i n the volume element Δ V) - (Rate o f mass o f substrate leaving the volume element Δ V) + (Net Rate o f generation o f mass o f substrate within the volume element Δ V) = 0

Q {S |}_{V} - Q {S |}_{V + Δ V} + R_{i} Δ V = 0

(2.38)

FIG. 3.

(a) Mass balance over differential volume element ΔV over the reactor height, (b) modified Levenspiel plot for recycle system.

Where a is the specific surface area (cm⁻¹) of the media on which biofilm is attached, R_i = -J_ia, J_i = mass flux (mg cm⁻² day⁻¹) in the biofilm in the ith compartment of the CMBR, and Q is the volumetric flow rate (cm³/day).

Dividing S|_V in the equation on both sides, we get

Q (\frac{{S |}_{V} - {S |}_{V + Δ V}}{{S |}_{V}}) - \frac{R_{i} Δ V}{{S |}_{V}} = 0

(2.39)

Substituting $Δ X = (\frac{{S |}_{V} - {S |}_{V + Δ V}}{{S |}_{V}})$ and $R_{i} = J_{i} a$

Q Δ X - \frac{J_{i} a Δ V}{{S |}_{V}} = 0

(2.40)

The Levenspiel equation can be obtained as follows:

\frac{Δ V}{Δ X} - Q \frac{{S |}_{V}}{J_{i} a} = 0

(2.41)

Taking limits $Δ V \to 0, and Δ X \to 0$ ,

\frac{d V}{d X} - Q \frac{{S |}_{V}}{J_{i} a} = 0

(2.42)

If the reactor has void volume fraction Φ, then the equation becomes

\frac{d V (1 - Φ)}{d X} - Q \frac{{S |}_{V}}{J_{i} a} = 0

(2.43)

Rearranging and separating the variables and integrating under limits X₁ = 0 at V₁ = 0, X₂ = X to V₂ = V, we get

\int d V = Q S_{i, o} α \int_{0}^{X} \frac{d X}{J_{i}},

(2.44)

The equation is known as the Levenspiel equation for computing the volume and area under the curve $\frac{{Q S}_{o} α}{J}$ and dX to be calculated. This graph is known as the Levenspiel plot (Levenspiel, 1998; Scott Fogler, 2020; Skogestad, 2009).

Where group constant,

α = \frac{1}{(1 - Φ) a}

Dividing the equation by Q into both sides converts the Levenspiel plot in HRT, θ, as follows. This can be considered as modified Levenspiel equation and is given by

\int d θ = S_{i, o} α \int_{0}^{X} \frac{d X}{J_{i}}

(2.45)

The curve between $\frac{S_{o} α}{J}$ and dX could potentially be plotted to create the equation’s modified Levenspiel graphic. The area covered by the modified Levenspiel map represents the reactor’s HRT. Here, the numerical integration technique can be used to assess the area under the modified Levenspiel plot; it applies Simpson’s 1/3^rd integration rule (Atkinson, 1989).

Recycle reactor model

As previously stated, the trickling filter can be seen as a PFR. For PFR with recycle, the recycling ratio (RR) is

R R = \frac{Volume o f fluid returned t o the reactor entrance}{Volume Leaving the system}

(2.46)

The range of the RR is zero to infinity. The reactor exhibits pure plug-flow behavior when RR = 0; at RR = ∞, the system acts more like a mixed-flow reactor. In a PFR, recycling results in reverse mixing.

Rearranging the reactor’s equation with recycling and using the mass balance is

\int d θ = S_{i, o} α \int_{X_{1}}^{X_{2}} \frac{d X}{J_{i}}

(2.47)

In this equation,

X_{1} = (\frac{R R}{R R + 1}) X_{2}

(2.48)

The recycle in the trickling filter results in back mixing. The HRT, θ calculated using the modified Levenspiel plot is as follows:

For the nonrecycled system, the integral dθ can be evaluated under the limits ${X = X}_{1} = 0$ to ${X = X}_{2}$

For recycle system, the integral dV is evaluated under the limits $X_{1} = (\frac{R}{R + 1}) X_{2}$ to ${X = X}_{2}$

The modified Levenspiel plot (Fig. 3b) demonstrates that the area under the curve reaches its maximum when the initial restriction is zero. When recycling is used, the initial limit shifts to the right, resulting in a smaller area under the curve with an increasing recycle ratio. As a result, less HRT is required to achieve the same conversion as if the system had not been recycled.

Solution procedure

The rate law expression equations obtain the solution of the biofilm model for heterotrophs (Eqs. 2.18, 2.24), nitrifiers (Eqs. 2.20, 2.25), and denitrifiers (Eqs. 2.34) for two types of biofilms (bilayer and homogeneous). The values of biofilm kinetic function Ψ(S_H*), Ψ(S_N*) $Ψ (S_{{NO}_{3}}^{*})$ are calculated from the correlations given in the literature. The B_H, B_N, and $B_{N O_{3}}$ and Ψ(S_H*), Ψ(S_N*), $Ψ (S_{{NO}_{3}}^{*})$ together give J*_H, J*_N, and $J_{N O_{3}}^{*}$ for each discretized CMBR compartment (Rittmann and McCarty, 2001; Wik, 2003). The reactor was discretized into 100 CMBR compartments. The model was solved using numerical integration using MATLAB 2018 package. For each aerobic and anoxic zone for Equation (2.43), J*_H, J*_N, and $J_{N O_{3}}^{*}$ can be converted into dimensional form as follows:

J_{F i} = \frac{J_{F, i}^{*} D_{i} K_{i}}{L_{i}}

(2.49)

Where i = H, N, NO₃

The algorithm for the modeling process is demonstrated in Fig. 4.

FIG. 4.

Algorithm for the modeling of MTF. MTF, modified trickling filter.

Results and Discussion

Theoretical performance calculations

The mass fluxes of COD and ammonia substrates are calculated using the rate law expressions described by Equations (2.18) and (2.20) for bilayer biofilm and by Equations (2.26) and (2.27) for homogeneous biofilm, respectively. Activated Sludge Model No. 1 contains the parameters needed to calculate the biofilm numbers for each species (Henze et al., 2015) and from the literature (Morgenroth, 2020). The experimental results of the biomass dispersion in the PU sponge utilized in a MTF that are used to calculate the biofilm thickness (Shukla, 2021) are mentioned in the table, and the biofilm numbers calculated are shown in Table 2.

Table 2.

Reactor Configuration and Biofilm Parameters Are Used to Determine Biofilm Numbers for Heterotrophs, Nitrifiers, and Denitrifiers*

Biofilm parameters	Heterotrophs	Nitrifiers	Denitrifiers (organic carbon)	Denitrifiers (nitrate)
L, Biofilm thickness (cm)	0.01495	0.01495	0.01495	0.01495
μ_max, Specific growth rate (day⁻¹)	6.0	0.8	6.0	6.0
K, Half-saturation constant (mg cm⁻³)	0.020	0.001	0.020	0.0005
Xf, Biomass density (mg cm⁻³)	35	35	35	35
D, Diffusivity in water (cm² day⁻¹)	0.60	1.70	0.60	1.60
Df, Diffusivity in biofilm (cm² day⁻¹)	0.48	1.36	0.48	1.28
Y, Yield (mg biomass mg⁻¹ substrate)^a	0.50	0.24	0.67	$2.86 (\frac{1 - Y_{H}}{Y_{H}})$
ε, Porosity (polyurethane foam)	0.80	0.80	0.80	0.80
Φ, Reactor void volume fraction	0.67	0.67	0.67	0.67
a, Specific surface area (cm⁻¹)	27.27	27.27	27.27	27.27
Granule diameter (cm)	0.0299	0.0299	0.0299	0.0299
Pore diameter (cm)	0.0550	0.0550	0.0550	0.0550
Volume fraction of one pore occupied by granule	0.1622	0.1622	0.1622	0.1622
Fraction of pores filled with granules	0.50	0.50	0.50	0.50
Biomass distribution within one granule (%)	70	20	70	70

*Henze et al., 1987; Shukla, 2021.

The yield coefficient (Y) represents the true biomass yield expressed as mg biomass (as VSS) formed per mg substrate (as COD) removed.

The heterotrophic yield is Y_h = 0.50. For nitrate, the yield value has been corrected following $Y_{N O_{3}} = 2.86 (\frac{1 - Y_{H}}{Y_{H}}),$ ensuring stoichiometric consistency between substrate oxidation and biomass growth.

The biofilm kinetic functions are Ψ(S_H*), Ψ(S_N*), and $Ψ (S_{{NO}_{3}}^{*})$ . Both models are evaluated from the correlations given in Table 3.

Table 3.

Biofilm Numbers for Heterotrophs and Nitrifiers

Biofilm model	Biofilm numbers
Aerobic zone	Heterotrophs		Nitrifiers
Bilayer biofilm	0.8962	1.0400^a	0.9012	1.0459^a
Homogeneous biofilm	2.2676	2.6313^a	1.9647	2.2799^a
Anoxic zone	Denitrifiers
Bilayer biofilm	1.1791	1.3627^a
Homogeneous biofilm	2.9833	3.0963^a

Biofilm numbers were calculated when microscopic pore diffusion effects were considered.

Using the technique described in the “Model for Porous Media” section, the theoretical removal efficiencies of the MTF at different HRTs were computed and compared with the experimental removal efficiencies of COD, NH₃-N, and NO₃-N at the identical HRT shown in Table 4.

Table 4.

Comparison of Experimental and Theoretical Removal Efficiencies in the Modified Trickling Filter Under Varying Hydraulic Residence Times and Inlet Concentrations

HRT (h)	Inlet COD (mg L⁻¹)	Exp COD (%)	Bilayer COD (%)	Bilayer COD + diff (%)	Homog. COD (%)	Homog. COD + diff (%)
4	119.2	76	82	100	100	100
4	140.1	74	77	100	97	100
4	161.3	71	72	99	90	100
2.5	119.2	47	56	79	68	100
2.5	140.1	46	53	73	63	100
2.5	161.3	45	49	67	58	100
2.5 + 0.5RR	119.2	73	73	100	93	100
2.5 + 0.5RR	140.1	68	70	97	89	100
2.5 + 0.5RR	161.3	61	66	75	81	100
2 + 0.5RR	119.2	60	62	87	77	100
2 + 0.5RR	140.1	56	59	82	73	94
2 + 0.5RR	161.3	52	56	64	67	88

COD removal efficiencies at different HRTs and inlet COD concentrations.

Theoretical predictions are shown for bilayer and homogeneous biofilm models, with and without pore diffusion. “+0.5RR” denotes operation at recycle ratio = 0.5.

HRT, hydraulic residence time; RR, relative error.

The theoretical removal efficiencies are estimated using the formulation of the following conditions: Case I:

Microscopic pore diffusion effects are not taken into account while assuming a bilayer biofilm (Fig. 5a).

FIG. 5.

Parity plots comparing theoretical and experimental removal efficiencies (X-axis: experimental, %; Y-axis: theoretical, %). (a) Bilayer model without pore diffusion effects; (b) bilayer model* with pore-scale diffusion; (c) homogeneous model without pore diffusion effects; (d) homogeneous model* with pore-scale diffusion. Blue diamonds indicate heterotrophs, green triangles nitrification, and red crosses denitrification. The solid line represents 1:1 parity. The asterisk (*) denotes inclusion of effective diffusivity.

Case II:

Microscopic pore diffusion effects are taken into account while assuming a bilayer biofilm (Fig. 5b).

Case III:

Microscopic pore diffusion effects are not taken into account while assuming a homogeneous biofilm (Fig. 5c).

Case IV:

Microscopic pore diffusion effects are taken into account while assuming a homogeneous biofilm (Fig. 5d).

Following the reactor configurations shown in Table 2, the mass fluxes for each column compartment were obtained. The % removal of each substrate (COD, ammonia, and nitrates) at different HRTs and inlet concentrations was computed using Equations (2.44) and (2.46). The void volume % of the column filled with the PU sponge is determined by the experimental results from the SEM pictures of the biomass immobilized in the sponge network (Shukla, 2021), and the group parameter for the reactor design is determined by α and comprises a specified surface area shown in Table 2. The theoretical and experimental removal efficiencies of each substrate at different HRTs and inlet concentrations are summarized in Table 4. The MTF was tested under a range of COD and ammonia loadings and at multiple HRTs to evaluate treatment effectiveness. For clarity, the results are divided into two subtables: Table 4 presents COD removal efficiencies, whereas Table 5 shows ammonia and nitrate removal efficiencies. In both subtables, experimental data are compared with predictions from the bilayer and homogeneous biofilm models, evaluated with and without pore diffusion. To avoid masking trends by high percent removals, the full datasets, including effluent concentrations and, for nitrogen species, complete mass balances (influent, effluent, and N removed), are reported in Supplementary Tables S1 and S2. These supplementary data allow direct evaluation of compliance with discharge standards and provide additional insight into denitrification performance.

Table 5.

Comparison of Experimental and Theoretical Removal Efficiencies in the Modified Trickling Filter Under Varying Hydraulic Residence Times and Inlet Concentrations

HRT(h)	Inlet NH₃-N (mg L⁻¹)	Inlet NO₃-N (mg L⁻¹)	Exp NH₃-N (%)	Exp NO₃-N (%)	Bilayer NH₃-N (%)	Bilayer NH₃-N + diff (%)	Bilayer NO₃-N (%)	Bilayer NO₃-N + diff (%)	Homog. NH₃-N (%)	Homog. NH₃-N + diff (%)	Homog. NO₃-N (%)	Homog. NO₃-N + diff (%)
4	28.8	57.6	100	55	100	100	59	100	100	100	100	100
4	31.9	62.3	98	40	100	100	55	100	90	100	100	100
4	35	65.9	94	25	94	100	52	100	88	100	100	100
2.5	28.8	29.2	62	55	71	95	72	100	100	100	100	100
2.5	31.9	29.8	54	45	65	87	69	100	100	100	100	100
2.5	35	31.2	48	34	59	80	66	100	100	100	100	100
2.5 + 0.5RR	28.8	35.6	91	84	100	100	60	100	95	100	95	100
2.5 + 0.5RR	31.9	37.5	96	80	94	100	57	100	92	100	92	100
2.5 + 0.5RR	35	39.4	82	72	86	100	55	100	88	100	88	100
2 + 0.5RR	28.8	35	87	89	84	100	48	100	80	100	80	100
2 + 0.5RR	31.9	35.6	81	83	76	97	44	100	76	100	76	100
2 + 0.5RR	35	36.8	76	76	71	90	45	97	74	100	74	100

Nitrogen species (NH₃-N and NO₃-N) removal efficiencies at different HRTs and inlet concentrations.

Theoretical predictions are shown for bilayer and homogeneous biofilm models, with and without pore diffusion. “+0.5RR” denotes operation at recycle ratio = 0.5.

HRT, hydraulic residence time; RR, relative error.

Experimental validation

The model was validated against published experimental data obtained from a 170 L pilot-scale MTF (Shukla, 2021). The available dataset included COD, NH₃-N, and NO₃-N removal efficiencies measured under different inlet concentrations and HRTs. These reported values were used directly for model validation, as replicate statistics (SD/SE) were not available in the source publications.

The theoretical and experimental removal efficiencies are summarized in Table 4 and illustrated in the form of parity charts in Fig. 5. Across all tested conditions, the bilayer biofilm model provided closer agreement with reported experimental removal efficiencies than the homogeneous biofilm model. When pore diffusion was incorporated, the bilayer model achieved up to ∼25% higher removal efficiency compared with the no-diffusion case. It should be noted that the composition of PU sponge media limits true microscopic pore diffusion; however, representing this effect in the model improved predictive agreement.

Model validation was quantified using Willmott’s (d) (Willmott, 1982) and RE, defined as the normalized RMSE with respect to the mean of the experimental data. These metrics provide a robust evaluation of model–experiment agreement in the absence of replicate-based variability measures (Stöckle et al., 2004; Willmott, 1982).

The computed values of Willmott’s Index and RE are shown in Table 6. Without taking into account the effects of pore diffusion, it has been found that the bilayer model has values that show a good model prediction with the experimental value. The bilayer biofilm is also justified when the difference in the kinetics of the two species is significantly higher than the outer layer that is occupied by rapid growers (for heterotrophs, μ_max = 6 d⁻¹), as they tend to use substrate more rapidly than slow growers (nitrifiers, μ_max = 0.8 d⁻¹) (Stöckle et al., 2004). The growth dynamics of both species make the homogeneous model impossible. The oxidation of ammonia into nitrates is carried out by two different groups of organisms. In this model, they were considered collectively as one species. This assumption is based on the kinetics of both species. The difference in the kinetics of these species is extremely insignificant; thus, they can exist evenly dispersed in the nitrifying layer beneath the heterotrophic layer (Blackburne et al., 2007). Figure 5 shows that the experimental values for denitrification are higher than the theoretical values.

Table 6.

Index of Agreement of Goodness and Relative Error for Model Predictions and Experimental Data

	COD				NH₃-N				NO₃-N
Parameters of agreement of goodness fit	Bilayer biofilm		Homogeneous biofilm		Bilayer biofilm		Homogeneous biofilm		Bilayer biofilm		Homogeneous biofilm
Willmott’s index, d	0.954	0.453^a	0.104	0.147^a	0.946	0.469^a	0.94	0.021^a	0.335	−0.771^a	−0.34	−0.008^a
Relative error, RE	0.0456	0.278^a	0.370	0.414^a	0.070	0.241^a	0.07	0.248^a	0.4125	0.308^a	0.393	0.437^a

The values are calculated while considering microscopic pore diffusion effects.

In the aerobic zone, ammonia is oxidized to generate nitrates (Shukla, 2021); according to the results of the research, the MTF contains a variety of denitrifying species, including around 10% of the anammox species found in the biomass that were separated from the anoxic zone. Anammox are slow-paced growers, with an extremely high doubling time of 29 days (van de Graaf, 1997). Adding explicit representations of biomass growth, decay, and wash-off would make the model better in a more direct way than adding anammox kinetics. These processes influence solid retention time (SRT), biofilm stability, and the competition between fast-growing heterotrophs and slower nitrifiers that dominate the MTF under aerobic and intermittently anoxic conditions. In the present study, the model was developed to describe heterotrophic oxidation of organic matter with simultaneous nitrification in the aerobic zone and heterotrophic denitrification in the anoxic zone, which represent the principal mechanisms of pollutant removal in the system. To make the analysis easier and focus on these main pathways, we assumed that the biofilm thickness and biomass density stayed the same. Anammox bacteria, on the other hand, grow very slowly (doubling time ≈ 29 days) and need long SRTs and stable anoxic environments, which are not what the MTF is currently doing. Consequently, forthcoming model enhancements will concentrate on integrating biomass growth, decay, and detachment dynamics before expanding the framework to encompass anammox under suitable environmental conditions. The fauna of the biofilm may become more diverse as a result of the PU sponge medium. In the model, only heterotrophic denitrification is taken into account. Pure heterotrophic denitrification is the basis for calculating the theoretical values. The observation suggests that autotrophic denitrification or other microbial activities inside the biofilm may be contributing to the denitrification process.

Classical multispecies biofilm and trickling filter models (Rittmann and McCarty, 2001; Wanner and Gujer, 1986) have long provided frameworks for stratification and substrate diffusion in nonporous media, but they typically omit detailed treatments of porous sponge media, tortuosity, or bilayer versus homogeneous structures. More recent studies have begun to address these gaps. For example, Bottomley and Eberl, (2024) developed a coupled PDE–ODE model for biofilm growth in porous media, showing that neglecting longitudinal diffusion can underestimate substrate depletion and required reactor length. Building on these insights, the present model advances beyond existing approaches by (1) incorporating an effective diffusion coefficient explicitly linked to sponge porosity and tortuosity, (2) comparing bilayer and homogeneous biofilm structures under identical conditions, (3) validating model predictions with experimental data using rigorous metrics such as Willmott’s (d) and RE, and (4) identifying design-relevant operating windows (recycle ratio ≈ 0.5, biofilm thickness 0.015–0.020 cm, porosity > 0.7) that can be directly applied to sponge-based MTF systems. One limitation is that replicate statistics were not available in the source datasets; future experimental studies with replicated measurements would enable confidence intervals to be established and further strengthen the robustness of model validation. In addition, future extensions of the model could incorporate slow-growing functional groups such as anammox by applying hybrid kinetic approaches, thereby capturing their long doubling times and potential contribution to nitrogen removal under anoxic microzones.

Summary and Conclusions

Two categories of biofilm structures are covered by the rate law expressions obtained from the research of Sharma et al. (2022): (1) homogeneous biofilm and (2) bilayer biofilm. When the kinetic speeds of the two species differ significantly, a multilayer structure is feasible due to competition for space; heterotrophs inhabit the outer layer, and nitrifiers occupy the interior layer. A homogeneous biofilm can form when the kinetics of the two species are similar or when there is no interaction between the two species. The axial part of a trickling filter functions as a CMBR. In a trickling filter, these CMBRs are linked together in series. Overall, the trickling filter functions as a PFR. The experimental findings verified the model as the top zone is composed of two species: heterotrophs and nitrifiers, whereas the bottom zone is composed solely of facultative heterotrophs. The bilayer model hypothesis in the upper zone provides promising results, with theoretical values estimated using the model and experimental data of acceptable correctness. With RE values of 0.04564 and 0.0703 for heterotrophs and nitrifiers, respectively, Willmott’s (d) is 0.954 and 0.946, indicating an ideal prediction.

The experimental nitrate removal efficiency for denitrification is greater than the theoretical values. Ammonia undergoes oxidation in the aerobic zone, resulting in nitrates. According to the experimental results of Shukla (2021), the MTF contains a variety of denitrifying species, including around 10% of the anammox species found in the biomass that was isolated from the anoxic zone. The biofilm’s flora may diversify as a result of the PU sponge medium. In the model, only heterotrophic denitrification is taken into account. Pure heterotrophic denitrification is the basis for calculating the theoretical values. It is likely that autotrophic denitrification or other microbial activities inside the biofilm are contributing to the observed denitrification. It has been noted that denitrification does not proceed by a single pathway, often heterotrophic denitrification; additional denitrifying processes in wastewater treatment are a promising area of research. Theoretically, the quantity of biofilms increases due to pore diffusion. Microscopic pore diffusion is not supported by the PU sponge media utilized in MTF, according to observations (Fig. 5b, d). The available specific area (27.27 cm²/cm³) for biomass growth and attachment is being macroscopically enhanced by the PU sponge. Under the conditions tested, predicted removal efficiencies with microscopic pore diffusion are ∼25% higher than those without microscopic pore diffusion.

From a design perspective, the model highlights several actionable levers for optimizing MTF performance:

Recycle ratio: The modified Levenspiel analysis (Section Recycle Reactor Model) shows that moderate recycle (RR ≈ 0.5) can reduce HRT requirements by 20–25% while maintaining comparable COD and NH₃-N removal. Excessive recycle, however, approaches mixed-flow behavior and reduces plug-flow efficiency.

Biofilm stratification: Maintaining bilayer structures (heterotrophs outside, nitrifiers inside) yields more accurate predictions and higher removal efficiency compared with assuming homogeneous distributions. This suggests that operating conditions that favor layered microbial growth (e.g., controlled loading rates, oxygen penetration depth) improve system stability.

Porosity thresholds: High porosity (ε ≈ 0.8 in the present sponge media) supports enhanced flux, but sensitivity analysis indicates that a drop of 15–20% in effective porosity significantly increases the biofilm number BBB, reducing substrate penetration. Periodic backwashing or media replacement may therefore be necessary to avoid pore occlusion.

Design range for thickness: The model indicates that biofilm thickness in the range of 0.015–0.020 cm provides a balance between substrate penetration and biomass retention; thicker biofilms risk diffusion limitation without proportional gains in removal.

In summary, the model not only reproduces COD, NH₃-N, and NO₃-N removal trends but also provides quantitative operating windows (recycle ratio ≈ 0.5, porosity > 0.7, biofilm thickness 0.015–0.020 cm) that can guide the design and operation of sponge-based MTFs toward reliable and energy-efficient treatment.

Authors’ Contributions

D.S.: Conceptualization, methodology, software, formal analysis, investigation, writing—original draft, and numerical modeling. S.G.: Data curation, validation, and writing—review and editing. R.C.: Resources and writing—review and editing. A.G.: Software and data visualization. Z.A.S.: Supervision, project administration, investigation, and resources. T.R.S.: Supervision, project administration, and investigation.

Footnotes

Author Disclosure Statement

The authors declare that they have no competing financial interests.

Funding Information

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Supplemental Material

References

Alsaeed

, Nelson

, Edwards

, et al. A mathematical model for the activated sludge process with a sludge disintegration unit. Chem Prod Process Model, 2023; 18(2):325–353; doi: 10.1515/cppm-2021-0064

Atkinson

. An Introduction to Numerical Analysis, Second Ed. Wiley; 1989.

Bakke

, Olsson

. Biofilm thickness measurements by light microscopy. J Microbiol Methods, 1986; 5(2):93–98; doi: 10.1016/0167-7012(86)90005-9

Blackburne

, Vadivelu

, Yuan

, et al. Determination of growth rate and yield of nitrifying bacteria by measuring carbon dioxide uptake rate. Water Environ Res, 2007; 79(12):2437–2445; doi: 10.2175/106143007X212139

Bottomley

, Eberl

. A PDE–ODE coupled model for biofilm growth in porous media that accounts for longitudinal diffusion and its effect on substrate degradation. Mathematics and Computation in Applications, 2024; 29(6):116; doi: 10.3390/mca29060116

Burrell

, Keller

, Blackall

. Characterisation of the bacterial consortium involved in nitrite oxidation in activated sludge. Water Sci Technol, 1999; 39(6):45–52; doi: 10.2166/wst.1999.0257

Chalipa

, Hosseinzadeh

, Nikoo

. Performance evaluation of a new sponge-based moving bed biofilm reactor for the removal of pharmaceutical pollutants from real wastewater. Sci Rep, 2024; 14(1):14240; doi: 10.1038/S41598-024-64442-5

Chatterjee

, Biswas

, Datta

, et al. Atomic force microscopy in biofilm study. Microscopy (Oxf), 2014; 63(4):269–278; doi: 10.1093/JMICRO/DFU013

Chou

, Wu

, Zeng

, et al. Evaluation of solute diffusion tortuosity factor models for variously saturated soils. Water Resour Res, 2012; 48(10); doi: 10.1029/2011WR011653

10.

Cooke

. Trickling filter ecology. Ecology, 1959; 40(2):273–291; doi: 10.2307/1930039

11.

Costerton

, Geesey

, Cheng

. How bacteria stick. Sci Am, 1978; 238(1):86–95; doi: 10.1038/scientificamerican0178-86

12.

Deguchi

, Kashiwaya

. Study on nitrified liquor recycling process operations using polyurethane foam sponge cubes as a biomass support medium. Water Sci Technol, 1994; 30(6):143–149; doi: 10.2166/wst.1994.0261

13.

Germain

. Economical treatment of domestic waste by plastic-medium trickling filters. J Water Pollut Control Fed, 1966; 38(2):192–203.

14.

Harris

, Hansford

. A study of substrate removal in a microbial film reactor. Water Res, 1976; 10(11):935–943; doi: 10.1016/0043-1354(76)90070-1

15.

Haug

, Mccarty

. Nitrification with submerged filters. Water Pollut Control Fed, 1972; 44(11):2086–2102.

16.

Henze

, Grady

CPL

, Gujer

, et al. A general model for single-sludge wastewater treatment systems. Water Res, 1987; 21(5):505–515; doi: 10.1016/0043-1354(87)90058-3

17.

Henze

, Gujer

, Mino

, et al. Activated sludge models ASM1, ASM2, ASM2d and ASM3. Water Intell Online, 2015; 5(0):9781780402369; doi: 10.2166/9781780402369

18.

Horan

. Handbook of Water and Wastewater Microbiology. Academic Press; 2004; doi: 10.5860/choice.41-4037

19.

Jennings

, Snoeyink

, Chian

ESK

. Theoretical model for a submerged biological filter. Biotechnol Bioeng, 1976; 18(9):1249–1273; doi: 10.1002/bit.260180906

20.

Jiang

, Zhang

, Lin

, et al. Modification of polyurethane sponge filler using medical stones and application in a moving bed biofilm reactor for ex situ remediation of polluted rivers. J Water Process Eng, 2021; 42:102189; doi: 10.1016/J.JWPE.2021.102189

21.

Kissel

, McCarty

, Street

. Numerical simulation of mixed‐culture biofilm. J Environ Eng, 1984; 110(2):393–411; doi: 10.1061/(ASCE)0733-9372(1984)110:2(393)

22.

Levenspiel

. Chemical Reaction Engineering, 3rd Edition. Wiley: New York; 1998; p. 704.

23.

Mattei

, Frunzo

, D’Acunto

, et al. Continuum and discrete approach in modeling biofilm development and structure: A review. J Math Biol, 2018; 76(4):945–1003; doi: 10.1007/s00285-017-1165-y

24.

Morgenroth

. Modelling biofilms. In: Biological Wastewater Treatment: Principles, Modeling and Design. ( Chen

; Ekama

; van Loosdrecht

MCM

; Brdjanovic

, eds.) IWA Publishing; 2020; pp. 757–812; doi: 10.2166/9781789060362_0757

25.

Mulcahy

, La Motta

. Mathematical model of the fluidized bed biofilm reactor. Enzyme Microb Technol, 1978; 4(4):269–275.

26.

Nuhoglu

, Keskinler

, Yildiz

. Mathematical modelling of the activated sludge process—the Erzincan case. Process Biochem, 2005; 40(7):2467–2473; doi: 10.1016/J.PROCBIO.2004.09.011

27.

Palmer

, Sternberg

. Modern microscopy in biofilm research: Confocal microscopy and other approaches. Curr Opin Biotechnol, 1999; 10(3):263–268; doi: 10.1016/S0958-1669(99)80046-9

28.

Payne

. Denitrification. Wiley: New York; 1981.

29.

Peters

, Alleman

. History of Fixed-Film Wastewater Treatment Systems. Iowa State University; 1982; pp. 60–88.

30.

Rittmann

. The significance of biofilms in porous media. Water Resour Res, 1993; 29(7):2195–2202; doi: 10.1029/93WR00611

31.

Rittmann

, Manem

. Development and experimental evaluation of a steady‐state, multispecies biofilm model. Biotechnol Bioeng, 1992; 39(9):914–922; doi: 10.1002/bit.260390906

32.

Rittmann

, McCarty

. Variable-order model of bacterial-film kinetics. J Envir Engrg Div, 1978; 104(5):889–900; doi: 10.1061/JEEGAV.0000825

33.

Rittmann

, McCarty

. Model of steady-state-biofilm kinetics. Biotechnol Bioeng, 1980; 22(11):2343–2357; doi: 10.1002/BIT.260221110

34.

Rittmann

, McCarty

. Environmental Biotechnology: Principles and Applications. McGraw-Hill; 2001.

35.

Schulze

BKL

. Load and efficiency of trickling filters. J Water Pollut Control Fed, 1960; 32(3):245–261.

36.

Scott Fogler

. Elements of Chemical Reaction Engineering. Sixth Ed. Pearson; 2020.

37.

Seymour

, Codd

, Gjersing

, et al. Magnetic resonance microscopy of biofilm structure and impact on transport in a capillary bioreactor. J Magn Reson, 2004; 167(2):322–327; doi: 10.1016/J.JMR.2004.01.009

38.

Sharma

, Sreekrishnan

, Ahammad

. Modeling biofilm and development of rate law expressions for biofilm kinetics. Chem Eng J Adv, 2022; 12:100419; doi: 10.1016/j.ceja.2022.100419

39.

Shukla

. Efficiency comparison of conventional wastewater treatment systems in removing emerging pollutants. Indian Institute of Technology: Delhi; 2021.

40.

Skogestad

. Chemical and Energy Process Engineering. Routledge; 2009; p. 436.

41.

, Naka

, Goel

, et al. Modelling clogging and biofilm detachment in sponge carrier media. Water Sci Technol, 2014; 69(6):1298–1303; doi: 10.2166/wst.2014.025

42.

Stöckle

, Kjelgaard

, Bellocchi

. Evaluation of estimated weather data for calculating penman-monteith reference crop evapotranspiration. Irrig Sci, 2004; 23(1):39–46; doi: 10.1007/s00271-004-0091-0

43.

van de Graaf

. Biological Anaerobic Ammonium Oxidation. Delft University Press; 1997.

44.

Versteeg

, Malalasekera

. An Introduction to Computational Fluid Dynamics : The Finite Volume Method. Pearson Education Ltd; 2007.

45.

Wang

, Chen

, Xu

, et al. Three-dimensional biofilm electrode reactors with polyurethane sponge carrier for highly efficient treatment of pharmaceuticals wastewater containing tetrahydrofuran. Water (Basel), 2022; 14(22):3792; doi: 10.3390/W14223792

46.

Wanner

, Eberl

, Morgenroth

, et al. Mathematical Modeling of Biofilms. Scientific and technical report no. 18. IWA Publishing; 2006.

47.

Wanner

, Gujer

. Competition in biofilms. Water Sci Technol, 1985; 17(2–3):27–44; doi: 10.2166/WST.1985.0117

48.

Wanner

, Gujer

. A multispecies biofilm model. Biotechnol Bioeng, 1986; 28(3):314–328; doi: 10.1002/bit.260280304

49.

Wanner

, Reichert

. Mathematical modeling of mixed-culture biofilms. Biotechnol Bioeng, 1996; 49(2):172–184; doi: 10.1002/(SICI)1097-0290(19960120)49:2<172::AID-BIT6>3.0.CO;2-N

50.

Watanabe

, Okabe

, Hirata

, et al. Simultaneous removal of organic materials and nitrogen by micro-aerobic biofilms. Water Sci Technol, 1995; 31(1):195–203; doi: 10.2166/wst.1995.0044

51.

Wik

. Trickling filters and biofilm reactor modelling. Rev Environ Sci Biotechnol, 2003; 2(2–4):193–212; doi: 10.1023/B:RESB.0000040470.48460.bb

52.

Williamson

, McCarty

. A model of substrate utilization by bacterial films. J Water Pollut Control Fed, 1976; 48(1):9–24.

53.

Willmott

. Some comments on the evaluation of model performance. Bull Amer Meteor Soc, 1982; 63(11):1309–1313; doi: 10.1175/1520-0477(1982)063<1309:SCOTEO>2.0.CO;2

54.

Wood

, Whitaker

. Diffusion and reaction in biofilms. Chem Eng Sci, 1998; 53(3):397–425; doi: 10.1016/S0009-2509(97)00319-9

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.03 MB

0.02 MB