### 1. Introduction

### 1.1. Mathematical Model

### 1.2. Model assumptions

### 1.3. Model Equations

*D*,

*Y*,

*X*

*,*

_{v}*δ*,

*μ*

_{max},

*K*

*,*

_{S}*K*

*,*

_{I}*α*,

*m*and

*x*symbolize the diffusion coefficient in the biofilm, the biomass yield coefficient, the dry cell density of the biofilm, the biofilm thickness, the maximum specific growth rate, the half-saturation constant, the inhibition constant for alpha-pinene in the presence of MeOH, the coefficient for the effect of MeOH on alpha-pinene biodegradation, the air-biofilm partition coefficient and the distance in the biofilm, respectively. While

*S*,

*S*

*and*

_{i}*C*denote the concentration in the biofilm, the concentration at the air-biofilm interface and the concentration in the air stream, respectively. The subscripts

*m*and

*P*represent MeOH and alpha-pinene, respectively.

*Biofilm phase model equations:*

##### (1)

$${D}_{m}\frac{{d}^{2}{S}_{m}}{d{x}^{2}}=\frac{{X}_{v}}{{Y}_{m}}\frac{{\mu}_{\text{max}(m)}{S}_{m}}{{K}_{S(m)}+{S}_{m}}$$##### (2)

$${D}_{p}\frac{{d}^{2}{S}_{p}}{d{x}^{2}}=\alpha \frac{{X}_{v}}{{Y}_{p}}\frac{{\mu}_{\text{max}(p)}{S}_{p}}{{K}_{S(p)}+{S}_{p}}$$##### (3)

$${S}_{m}=\frac{{C}_{m}}{{m}_{m}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}and\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{S}_{p}=\frac{{C}_{p}}{{m}_{p}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x=0$$##### (4)

$$\frac{d{S}_{m}}{dx}=\frac{d{S}_{p}}{dx}=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x=\delta $$*U*

*,*

_{g}*A*

*,*

_{s}*h*and

*C*

*stand for the air superficial velocity, the biofilm surface area per unit volume of the BF, the position along the BF height and the concentration in the inlet air stream, respectively.*

_{i}*Gas phase model equations:*

##### (7)

$${C}_{m}={C}_{i(m)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}h=0$$##### (8)

$${C}_{p}={C}_{i(p)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}h=0$$*Dimensionless biofilm phase model equations:*

##### (9)

$$\frac{{d}^{2}{S}_{m}^{*}}{d{\sigma}^{2}}={\phi}_{1}\frac{{S}_{m}^{*}}{1+{\beta}_{1}{C}_{m}^{*}{S}_{m}^{*}}$$##### (10)

$$\frac{{d}^{2}{S}_{p}^{*}}{d{\sigma}^{2}}=\alpha {\phi}_{2}\frac{{S}_{p}^{*}}{1+{\beta}_{2}{C}_{p}^{*}{S}_{p}^{*}}$$##### (11)

$${S}_{m}^{*}={S}_{p}^{*}=1\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\sigma =0$$##### (12)

$$\frac{d{S}_{m}^{*}}{d\sigma}=\frac{d{S}_{p}^{*}}{d\sigma}=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}at\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\sigma =1$$*Dimensionless gas phase model equations:*

##### (13)

$$\frac{d{C}_{m}^{*}}{d\xi}={\kappa}_{1}{C}_{m}^{*}{\left(\frac{d{S}_{m}^{*}}{d\sigma}\right)}_{\sigma =0}$$##### (14)

$$\frac{d{C}_{p}^{*}}{d\xi}={\kappa}_{2}{C}_{p}^{*}{\left(\frac{d{S}_{p}^{*}}{d\sigma}\right)}_{\sigma =0}$$### 2. Presentation of the Iterative Scheme

*et al*. [28] to find approximated solutions for linear and nonlinear third order ODEs that arise in physical applications. The advantages of using this iterative scheme include fast convergence to the solution; and it gives highly accurate results without the need for spatial discretization or any assumptions that might change the physical representation of the model.

### 2.1. The Iterative Scheme for the Solution of First Order ODEs

*L*is a linear differential operator representing the first order derivative and

*f*(

*x*) is the inhomogeneous term.

*y*

*is the homogeneous solution which satisfies*

_{h}*L*[

*y*] = 0 and the initial condition is given by (18). While

*y*

*is the particular solution which satisfies*

_{o}*L*[

*y*] =

*f*(

*x*) and the following initial condition

*B*[

*y*] = 0.

*G*(

*x*,

*s*) is defined as:

*x*≠

*s*, the Green’s function is a homogeneous solution to (17). However, when

*x*=

*s*the Green’s function has singular behavior. Therefore, it comprises of two identical solutions for each side of the Dirac delta function given as:

##### (23)

$$G(x,s)=\{\begin{array}{ll}{e}_{1}\hspace{0.17em}\text{exp\hspace{0.17em}}(-\underset{s}{\overset{x}{\int}}p\hspace{0.17em}(x)\hspace{0.17em}dx)\hfill & a<x<s\hfill \\ {e}_{2}\hspace{0.17em}\text{exp\hspace{0.17em}}(-\underset{s}{\overset{x}{\int}}p\hspace{0.17em}(x)\hspace{0.17em}dx)\hfill & x>s\hfill \end{array}$$*e*

_{1}and

*e*

_{2}are constants; these constants are to be found using the following properties:

*G*(

*x*,

*s*) satisfies the homogeneous initial condition:

*G*(

*x*,

*s*) at

*x*=

*s*:

##### (25)

$${e}_{2}\hspace{0.17em}\text{exp}\left(-\underset{s}{\overset{x}{\int}}p\hspace{0.17em}(x)\hspace{0.17em}dx\right)-{e}_{1}\hspace{0.17em}\text{exp}\left(-\underset{s}{\overset{x}{\int}}p\hspace{0.17em}(x)\hspace{0.17em}dx\right)=1$$##### (26)

$$K[y]={y}_{h}+\underset{a}{\overset{\infty}{\int}}G(x,s)[{y}^{\prime}(s)+p(s)y(s)]\hspace{0.17em}ds$$##### (27)

$$K[y]={y}_{h}+\underset{a}{\overset{\infty}{\int}}G(x,s)[{y}^{\prime}(s)+p(s)y(s)-f(s)]\hspace{0.17em}ds+\underset{a}{\overset{\infty}{\int}}G(x,s)f(s)\hspace{0.17em}ds$$##### (28)

$$K[y]={y}_{h}+\underset{a}{\overset{\infty}{\int}}G(x,s)\hspace{0.17em}[{y}^{\prime}(s)+p(s)y(s)-f(s)]\hspace{0.17em}ds+y-{y}_{h}$$##### (29)

$$K[y]=y+\underset{a}{\overset{\infty}{\int}}G(x,s)\hspace{0.17em}[{y}^{\prime}(s)+p(s)y(s)-f(s)]\hspace{0.17em}ds$$##### (30)

$${y}_{n+1}=(1-{\alpha}_{n}){y}_{n}+{\alpha}_{n}K[y]\mathrm{\hspace{0.17em}\u200a\u200a}\mathrm{\hspace{0.17em}\u200a\u200a}\mathrm{\hspace{0.17em}\u200a\u200a}for\mathrm{\hspace{0.17em}\u200a\u200a}\mathrm{\hspace{0.17em}\u200a\u200a}\mathrm{\hspace{0.17em}\u200a\u200a}n\ge 0$$##### (31)

$${y}_{n+1}={y}_{n}+{\alpha}_{n}\underset{a}{\overset{\infty}{\int}}G(x,s)\hspace{0.17em}[{{y}^{\prime}}_{n}(s)+p(s){y}_{n}(s)-f(s)]\hspace{0.17em}ds$$*n*= 0, the starting function

*y*

_{0}for the iterations is selected to be the solution to the homogeneous equation

*L*[

*y*] = 0 subject to (18). Due to the lack of proper justification on the choice of the sequence

*α*

*, it is approximated by trying arbitrary values close to one and compare the results to find the value that meets a certain specified tolerance.*

_{n}### 2.2. The Iterative Scheme for the Solution of Second Order ODEs

*L*is a linear differential operator representing the second order derivative and

*f*(

*x*) is the inhomogeneous term.

*y*

*is the homogeneous solution which satisfies*

_{h}*L*[

*y*] = 0 and the boundary conditions given by (33) and (34). While

*y*

*is the particular solution which satisfies*

_{o}*L*[

*y*] =

*f*(

*x*) with the following boundaries

*B*

*[*

_{a}*y*] =

*B*

*[*

_{b}*y*] = 0.

*G*(

*x*,

*s*) is defined as:

*y*

_{1}and

*y*

_{2}form the fundamental set of solutions for

*L*[

*y*] = 0, which is the homogeneous solution. For

*x*≠

*s*, the Green’s function is a homogeneous solution to (32); however, when

*x*=

*s*the Green’s function has singular behavior. Therefore, the Green’s function comprises of two identical solutions for each side of the Dirac delta function given as:

##### (39)

$$G(x,s)=\{\begin{array}{ll}{f}_{1}{y}_{1}+{f}_{2}{y}_{2}\hfill & a\le x<s\hfill \\ {g}_{1}{y}_{1}+{g}_{2}{y}_{2}\hfill & s<x\le b\hfill \end{array}$$*f*

_{1},

*f*

_{2},

*g*

_{1}and

*g*

_{2}are constants; these constants are to be found using the following properties:

*G*(

*x*,

*s*) satisfies the homogeneous boundary conditions:

*G*(

*x*,

*s*) at

*x*=

*s*:

*G*(

*x*,

*s*) at

*x*=

*s*:

##### (42)

$${g}_{1}{{y}_{1}}^{\prime}(s)+{g}_{2}{{y}_{2}}^{\prime}(s)-{f}_{1}{{y}_{1}}^{\prime}(s)-{f}_{2}{{y}_{2}}^{\prime}(s)=1$$##### (43)

$$K[y]={y}_{h}+\underset{a}{\overset{b}{\int}}G(x,s)\hspace{0.17em}[{y}^{\u2033}(s)+q(s){y}^{\prime}(s)+r(s)y(s)]\hspace{0.17em}ds$$##### (44)

$$K[y]={y}_{h}+\underset{a}{\overset{b}{\int}}G(x,s)\hspace{0.17em}[{y}^{\u2033}(s)+q(s){y}^{\prime}(s)+r(s)y(s)-f(s)]\hspace{0.17em}ds\underset{a}{\overset{b}{\int}}G(x,s)f(s)\hspace{0.17em}ds$$##### (45)

$$K[y]={y}_{h}+\underset{a}{\overset{b}{\int}}G(x,s)\hspace{0.17em}[{y}^{\u2033}(s)+q(s){y}^{\prime}(s)+r(s)y(s)-f(s)]\hspace{0.17em}ds+y-{y}_{h}$$##### (46)

$$K[y]=y+\underset{a}{\overset{b}{\int}}G(x,s)\hspace{0.17em}[{y}^{\u2033}(s)+q(s){y}^{\prime}(s)+r(s)y(s)-f(s)]\hspace{0.17em}ds$$##### (47)

$${y}_{n+1}=(1-{\alpha}_{n}){y}_{n}+{\alpha}_{n}K[y]\mathrm{\hspace{0.17em}\u200a\u200a}\mathrm{\hspace{0.17em}\u200a\u200a}\mathrm{\hspace{0.17em}\u200a\u200a}for\mathrm{\hspace{0.17em}\u200a\u200a}\mathrm{\hspace{0.17em}\u200a\u200a}\mathrm{\hspace{0.17em}\u200a\u200a}n\ge 0$$### 2.3. Construction of the Iterative Scheme for the Biofilter Model Equations

*The iterative schemes for the biofilm phase model equations are:*

##### (49)

$$\begin{array}{c}{S}_{m(n+1)}^{*}={S}_{m(n)}^{*}+{\alpha}_{n}[{\int}_{0}^{\sigma}(-s)\hspace{0.17em}(\frac{{d}^{2}{S}_{m(n)}^{*}}{d{\sigma}^{2}}-\frac{{\phi}_{1}{S}_{m(n)}^{*}}{1+{\beta}_{1}{C}_{m(n)}^{*}{S}_{m(n)}^{*}})ds+\\ +{\int}_{\sigma}^{1}(-\sigma )\hspace{0.17em}(\frac{{d}^{2}{S}_{m(n)}^{*}}{d{\sigma}^{2}}-\frac{{\phi}_{1}{S}_{m(n)}^{*}}{1+{\beta}_{1}{C}_{m(n)}^{*}{S}_{m(n)}^{*}})ds]\end{array}$$##### (50)

$$\begin{array}{c}{S}_{p(n+1)}^{*}={S}_{p(n)}^{*}+{\alpha}_{n}[{\int}_{0}^{\sigma}(-s)\hspace{0.17em}(\frac{{d}^{2}{S}_{p(n)}^{*}}{d{\sigma}^{2}}-\frac{\alpha {\phi}_{2}{S}_{p(n)}^{*}}{1+{\beta}_{2}{C}_{p(n)}^{*}{S}_{p(n)}^{*}})ds+\\ {\int}_{\sigma}^{1}(-\sigma )\hspace{0.17em}(\frac{{d}^{2}{S}_{p(n)}^{*}}{d{\sigma}^{2}}-\frac{\alpha {\phi}_{2}{S}_{p(n)}^{*}}{1+{\beta}_{2}{C}_{p(n)}^{*}{S}_{p(n)}^{*}})ds]\end{array}$$*The iterative schemes for the gas phase model equations are:*

##### (51)

$${C}_{m(n+1)}^{*}={C}_{m(n)}^{*}+{\alpha}_{n}[{\int}_{0}^{\xi}(-1)\hspace{0.17em}(\frac{d{C}_{m(n)}^{*}}{d\xi}-{\kappa}_{1}{C}_{m(n)}^{*}{\left(\frac{d{S}_{m}^{*}}{d\sigma}\right)}_{\sigma =0})\hspace{0.17em}ds]$$### 3. Simulation Results and Discussion

### 3.1. Biodegradation of Methanol and Alpha-pinene

### 3.2. Biofiltration of Methanol

^{−4}is imposed on the residual error as a stopping criterion for the iterations. Several simulations runs have been done with a different number of iterations and various sequence

*α*

*values. The corresponding maximum residual errors were estimated along the BF height to determine the number of iterations and the proper sequence setting that yields to the specified tolerance. The maximum residual errors for 15, 20 and 25 iterations at*

_{n}*α*

*equal to 0.6, 0.4 and 0.2 are reported in Tables S1 to S3. The results showed that at a given value of*

_{n}*α*

*, the maximum residual error decreases with increasing number of iterations. To plot the concentration profile of MeOH, 15 iterations and a sequence of 0.4 were used as they meet the specified stopping criteria as mentioned above. Fig. 2 shows that the proposed technique yields a highly accurate solution and it can describe the biofiltration of MeOH when compared against the original model solution and the experimental data from the literature.*

_{n}### 3.3. Biofiltration of Alpha-pinene

*α*

*were determined by carrying out several simulations runs to meet a tolerance of 10*

_{n}^{−4}. These simulations were also done to study the impact of the number of iterations and the sequence value on the residual error. The residual error results are shown for 50, 55 and 60 iterations at

*α*

*equal to 0.2, 0.18 and 0.1 in Tables S4 to S6. The results showed that the residual error improves when more iteration are carried out.*

_{n}### 3.4. Sensitivity Analysis

### 3.5. Effect of Empty Bed Residence Time of the RE

^{−3}. RE increased rapidly when the EBRT was varied from 5 s to around 30 s. For short EBRT of around 10 s and less, the BF exhibited low RE of MeOH (≤ 60%). This is attributed to the fact that, the contact between the biomass and MeOH was too quick at which the microorganisms had insufficient time to perform the required degradation on the available amount of MeOH. In contrast, for long EBRT the RE increased until it reached a complete elimination of MeOH at an EBRT of 50 s. The increase in the RE can be explained by the fact that the time spent by MeOH in the BF is sufficiently large for the process to effectively remove the pollutant.

### 3.6. Effect of Biofilm Surface Area on the RE

^{−3}and an EBRT of 50 s in Fig. 5. The results showed that the RE is an increasing function of the biofilm surface area. The analysis demonstrated that a complete removal of MeOH can be accomplished at a biofilm surface area of 63 m

^{2}.m

^{−3}. This is intuitively expected because for a given biofilm thickness increased surface area increases the reaction volume and area for mass transfer.

### 3.7. Effect of Air-biofilm Partition Coefficient on the RE

^{−3}and an EBRT of 50 s. Although the air-biofilm partition coefficient is constant for a given compound, this parameter could be altered by factors such as temperature and surfactant addition. Therefore, it is essential to examine the effect of this parameter on the RE. At low values of air-biofilm partition coefficient of up to 0.01, the BF showed 100% removal of MeOH. This phenomenon is due to the fact that at low values of the air-biofilm partition coefficient, the solubility of MeOH in the biofilm increases resulting in rapid biodegradation in the biofilm. However, the increase in the air-biofilm partition coefficient to 0.1 substantially reduced the RE to 40.87%. This analysis confirms that the RE of the biological system diminishes as the air-biofilm partition coefficient increases since this hinders the diffusion of pollutants into the biofilm where the elimination process occurs.

### 3.8. Effect of Diffusion Coefficient on the RE

^{−6}m

^{2}.h

^{−1}to 6.00 × 10

^{−6}m

^{2}.h

^{−1}at an inlet MeOH concentration of 1 g.m

^{−3}and EBRT of 50 s. The simulation results show that the RE increases as the diffusion coefficient is increased. This behavior is expected as at high values of the diffusion coefficient, MeOH diffuses faster into the biofilm where the biodegradation takes place. A maximum RE of 98.06% was noticed at a diffusion coefficient of 6.00 × 10

^{−6}m

^{2}.h

^{−1}. While at low diffusivity of 1.00 × 10

^{−6}m

^{2}.h

^{−1}the BF removed about 82.83% of the MeOH.