A membrane fouling model based on pore adsorption

1 Introduction

Single mechanism fouling models are unable to fit complex experimental fouling data. First, the four classical fouling models (Complete, Standard, Intermediate, Cake) are unable to fit data where a non-linear dependence of the blocking constant on the foulant concentration is observed (Lee et al., 2013). All four models assume the blocking constant increases linearly with foulant concentration (Bolton G. et al., 2006).

Second, the four classical fouling models are unable to fit data with a high or low fouling index. The fouling index n relates the second and first inverse derivatives of volume versus time: ∂ 2 t ∂ V 2 ∝ ∂ t ∂ V n. It varies from zero to two for the four classical fouling models (Bolton G. R. et al., 2006). Fouling data exhibiting a negative fouling index have been observed experimentally (Lee et al., 2013; Ho and Zydney, 2000; Duclos-Orsello et al., 2006; Iritani et al., 2010; Costa et al., 2006; Bowen et al., 1995; Iritani, 2013; Cogan and Chellam, 2009) and has typically been fitted by assuming fouling occurs through two or three mechanisms sequentially using complex integral or approximate equations (Lee et al., 2013; Bolton G. et al., 2006; Bolton G. R. et al., 2006; Ho and Zydney, 2000; Duclos-Orsello et al., 2006).

To fit more complex data, a new model of membrane fouling was developed. This model assumes that foulant adsorption occurs because of a reaction between foulants and membrane pore surface area, reducing pore radius and increasing membrane resistance. The model utilizes the Hagen-Poiseuille law for flux, as has been done for previous models (Iritani, 2013; Bolton et al., 2005; Koonani and Amirinejad, 2019; Hwang et al., 2009; Du et al., 2022).

The rate of increase of foulant volume inside pores (V _*F*) was assumed to be governed by the reaction rate law d V F / d t = K C x S z, where K is the rate constant, C the foulant concentration, S the pore area, x the reaction order with respect to foulant concentration, and z the reaction order with respect to pore area. The dependence of pore radius on time was inserted into the Hagen-Poiseuille law for flux to obtain explicit equations to predict pore size, flux, resistance, pressure and volume versus time. The explicit equations allow the dependence on foulant concentration to be increased, decreased or eliminated by assuming the reaction order with respect to concentration x is above one, below one, or zero. The model was extended to treat membranes composed of multiple pore radii allowing flow in parallel.

Surprisingly, the four classical membrane fouling equations could be derived from the new model using specific reaction orders with respect to pore area (z = 1, 3, 5, 9, respectively), but with the improvement of allowing for a non-linear dependence on foulant concentration. The new model is therefore useful for determining if the optimal data fit will be provided by one of the classical fouling models by using these specific reaction orders, or by using other reaction orders. However, values of z above ∼3 are not physically realistic and are therefore considered fitting parameters.

The new model can fit a wider range of data than the four classical fouling models and is the first single mechanism fouling model capable of fitting data with a negative fouling index. The model provided excellent fits of flux versus time data using hydrophilic and hydrophobic PVDF membranes over a range of bovine serum albumin (BSA) concentrations, and flux versus time using hydrophilic PVDF membrane over a range of polyethylene glycol (PEG) concentrations. These data sets exhibited either a non-linear dependence on foulant concentration or a negative fouling index. These data sets had been fit previously with integral or approximate equations using four or more parameters and assuming fouling occurs through two or three mechanisms sequentially. The new model demonstrated broad utility by fitting these data using simple explicit equations, a single fouling mechanism, and only two to four fitted parameters.

This work summarizes the derivation of the new model, models with specific reaction orders with respect to membrane pore area, and the fit of these models to experimental data.

2 Modeling

2.1 Zeroth-order adsorption model

The adsorptive fouling model was previously derived by assuming that foulants react with the walls of membrane pores and are deposited, thereby constricting the pores and reducing flow. The reaction between foulants and pore walls was assumed to have a zeroth-order dependence on pore area (Bolton G. R. et al., 2006). The rate of deposition of foulant volume V _*F* on pore walls was described by:

d V F d t = K ′ C x S 0 (1)

Here K ′ is the adsorption rate constant; C the concentration of foulant; x is the reaction order with respect to foulant concentration; and S the membrane pore surface area. This is the same as Equation 6 in (Bolton G. R. et al., 2006), where the superscript 0 was erroneously omitted in the original publication. The initial pore radius r _*0* will decline to a reduced value r as foulants react and adsorb on pore walls. At any time, the pore area S = 2 π r L, where L is the pore length. The deposited foulant volume is given by V F = π r 0 2 − r 2 L, and d V F = − 2 π L d r. Assuming K A = K ′ / 2 π L r 0 and solving Equation 1:

r = r 0 1 − K A C x t (2)

Based on the Hagen-Poiseuille law, the flow of an incompressible uniform viscous liquid through a cylindrical tube with a constant circular cross-section will have a fourth-order dependence on pore radius. Equation 2 can be substituted into the Hagen-Poiseuille law to yield flux as a function time, where J = ∂ V / ∂ t:

J J 0 = r r 0 4 = 1 − K A C x t 4 (3)

Equation 3 can be integrated to determine volume as a function of time:

V = J 0 5 K A C x 1 − 1 − K A C x t 5 (4)

Equations 3, 4 are only applicable until the pore radius reaches zero, i.e., when K A C x t = 1. Beyond this point the flux will be zero and the volume will be equal to J 0 / 5 K A C x. Equation 4 can be inserted into Equation 3 to determine flux as a function of volume using Equation 5:

J J 0 = 1 − 5 K A C x V J 0 4 5 (5)

2.2 Derivation of the new fouling model

To improve upon previous work, a new adsorptive membrane fouling model has been developed, allowing for different reaction orders in terms of foulant concentration (x) and pore area (z). The rate of increase of foulant volume V _*F* on pore walls is governed by the reaction rate law:

d V F d t = K ′ C x S z (6)

Assuming V F = π r 0 2 − r 2 L, and inserting the pore area S = 2 π r L and d V F = − 2 π L d r into Equation 6 yields Equations 7, 8:

− 2 π L d r d t = K ′ C x 2 π r L z (7)

d r r z = − K ′ C x 2 π L z − 1 d t (8)

Integrating Equation 8:

1 r z − 1 = 1 r 0 z − 1 + z − 1 K ′ 2 π L z − 1 C x t (9)

r r 0 = 1 + z − 1 r 0 z − 1 K ′ 2 π L z − 1 C x t − 1 z − 1 (10)

Simplifying Equation 10 by introducing a new constant K = r 0 z − 1 K ′ 2 π L z − 1 yields an equation for pore radius versus time:

r r 0 = 1 + z − 1 K C x t − 1 z − 1 (11)

Using z = 1 in Equation 11 would result in a division by zero. In this case, Equation 8 can be solved to yield r = r 0 ⁡ exp ⁡ − K C x t. For z ≥ 1, the pore radius will decline but never reach zero, but for z < 1 the pore radius will reach zero at 1 − z − 1 K C x t = 0 and flow will stop.

Equation 11 can be substituted into the Hagen-Poiseuille law to yield flux versus time, J J 0, when operating and constant pressure, or the inverse of pressure versus time, P 0 P, when operating at constant flow rate:

P 0 P or J J 0 = r r 0 4 = 1 + z − 1 K C x t − 4 z − 1 (12)

Substituting J = d V d t and integrating to derive an equation for volume versus time yields Equations 13, 14:

V = ∫ J 0 1 + z − 1 K C x t − 4 z − 1 d t (13)

V = J 0 z − 5 K C x 1 + z − 1 K C x t z − 5 z − 1 − 1 (14)

Inserting Equation 14 into Equation 12 yields an equation for flux in terms of volume:

J J 0 = 1 + z − 5 K C x V J 0 − 4 z − 5 (15)

Using Equation 14 with z = 1 or z = 5 would result in division by zero. Explicit equations for these cases were derived from Equation 8–12 and are included in Table 1.

TABLE 1. New adsorption fouling models describing volume versus time, or flux or pressure versus time or volume, for different values of the reaction order with respect to pore area z. The fouling index is also provided as a function of z.

Typically the fouling rate constant K is assumed to increase or decrease linearly with foulant concentration (Bolton G. et al., 2006). Here, through the term C x, the dependence of fouling on concentration can be increased by assuming the reaction order with respect to concentration x is above 1, reduced by assuming x is between 0 and 1, or eliminated by assuming x is zero. A non-linear dependence of fouling on concentration has been observed experimentally (Duclos-Orsello et al., 2006), and was fit using the new models as described below.

Using different integer values of z in Equations 12, 14 allows generation of a range of mathematically explicit fouling models including the four classical fouling models. For z = 1, the model reduces to the equation for Complete Blocking (Bolton G. et al., 2006; Hermia, 1982). For z = 3, it reduces to Standard Blocking. For z = 5, it reduces to Intermediate Blocking. For z = 9, it reduces to Cake Fouling, though this high value of z is not physically realistic. This is the first time the four classical models have been derived assuming a reaction between foulants and pore walls. The new models, however, improve upon the four classical fouling models by allowing a non-linear dependence on foulant concentration through the parameter x.

2.3 Fouling model for multiple pore sizes

The new model was extended to treat membranes containing multiple pore radii (r _*a* , r _*b,* etc.). For flow through different pore sizes in parallel:

J J 0 = n a r a 4 + n b r b 4 + … n a r a 0 4 + n b r b 0 4 + … = ∑ n i r i 4 ∑ n i r i 0 4 (16)

Where n _*i* is the number of pores of radius i. Equation 11 is applicable to each pore radius r _*i*:

r i = r i 0 1 + z − 1 K i C x t − 1 z − 1 (17)

The reaction order in terms of membrane pore area z will be applicable to all pore sizes given that all pores have the same surface properties. The blocking constant K will be unique to each pore size as it depends on the specific initial pore radius r _*i* (K i = r i 0 z − 1 K ′ 2 π L z − 1.

The term f _*i* can be defined as the fraction of initial flow going through pores of radius r _*i* using Equation 18:

f i = n i r i 0 4 n a r a 0 4 + n b r b 0 4 + … (18)

Equation 16 can be updated by substituting Equation 17 for pore radius r _*i*:

P 0 P or J J 0 = ∑ i f i 1 + z − 1 K i C x t − 4 z − 1 (19)

Substituting J = d V d t into Equation 19 and integrating generates an equation for volume versus time:

V = ∑ i J 0 f i z − 5 K i C x 1 + z − 1 K i C x t z − 5 z − 1 − 1 (20)

For the case of two pore radii (r _*a* and r _*b*), f a + f b = 1 and Equations 19, 20 simplify to:

P 0 P or J J 0 = f a 1 + z − 1 K a C x t − 4 z − 1 + 1 − f a 1 + z − 1 K b C x t − 4 z − 1 (21)

V = J 0 z − 5 [ f a K a C x 1 + z − 1 K a C x t z − 5 z − 1 − 1 + 1 − f a K b C x 1 + z − 1 K b C x t z − 5 z − 1 − 1 ] (22)

Equations 21, 22 use five fitted parameters (z, x, f _*a* , K _*a* , K _*b*) versus three parameters (x, z, K) for the single pore size model (Equation 12).

The new model, and models for other specific values of z are summarized in Table 1.

3 Results and discussion

3.1 Volume versus time for different reaction orders with respect to pore area z

The performance of the new model was compared to those of the four classical fouling models and the previous adsorptive model (Bolton G. R. et al., 2006). For all models, the initial flux J _*0* was assumed to be equal to one and the value of the blocking constant K was varied until a volume of one was achieved at a filtration time of 10 seconds. Different values of the reaction order with respect to pore area z were used in Equation 14. The data are summarized in Figure 1.

Equation 14 will provide fits equal to those of the classical models by assuming specific values of z (1, 3, 5, 9) but can also fit a wider range of data depicted in the shaded area in Figure 1.

3.2 Fit of model to experimental flux versus time data

To evaluate the utility of the model, experimental flux versus time data from literature corresponding to a diversity of membranes, foulant, and foulant concentrations were fit. Equation 12 was fit to data from Figure 3B from literature (Lee et al., 2013) which depicts flux versus time for the fouling of hydrophilic PVDF membrane by different concentrations of polyethylene glycol (PEG). In addition, Equation 12 was fit to data from Figures 3A, 5A from literature (Duclos-Orsello et al., 2006), which depict flux versus time using hydrophilic and hydrophobic PVDF membranes over a range of bovine serum albumin (BSA) concentrations. Confocal microscopy has demonstrated that BSA internally fouls similar microporous membranes{Ferrando, 2005 #2131}{Kanani, 2008 #2132}. Other data from these publications were fit well but not included here either because the data were similar or the model systems were not representative of internal adsorptive fouling. The data from the three graphs were digitized using plotdigitizer.com with a precision of under 0.1%. The best fit was determined by minimizing the sum of squared residuals (SSR), where the residual was equal to the difference between a data point and the model prediction. Integer values of z and x were chosen when the best fit values were near integer values. The SSR exhibited a unique minimum over the parameter space in each case, indicating a unique best fit was obtained. In addition, the Root Mean Square Error (RMSE) and Degrees of Freedom Error (DFE) were evaluated.

3.2.1 Polyethylene glycol filtered through PVDF membranes

The data for the filtration of 20 kDa polyethylene glycol through 0.22 μm hydrophilic (based on the provided contact angle) PVDF membranes at concentrations from 50 to 1,000 mg/L (Lee et al., 2013) was fit using Equation 12. In this data set, the flux declines to 50% of the initial value in only about a third of the time for the high PEG concentration condition as compared to the low concentration, despite the 20-fold difference in concentrations (1,000 vs. 50 mg/L). This low dependence on foulant concentration may be due to the low interaction between the hydrophilic PEG and the hydrophilic PVDF membrane. It makes the data set an interesting test case for the non-linear dependence of foulant concentration enabled in Equation 12 using different values of x.

The new model provided a good fit to the four data curves using only 3 fitted parameters: z = 1.46 ± 0.10, x = 0.32 ± 6.7 × 10⁻³, and K = 1.46 × 10⁻³ ± 5.2 × 10⁻⁵, as the R ² value was 0.99, as shown in Figure 2. The SSR and RMSE were 0.167 and 0.0287. The DFE was 199, indicating the model was not overfitting the data. The new model fit the low dependence on foulant concentration using a value of x below one (0.32).

The previous work fitted the data to a model that assumed two different fouling mechanisms occurred sequentially and used an equation that approximated an integral solution. Seven parameters were used to fit each of the four curves, for a total of twenty eight parameters (Lee et al., 2013). The new model fit the data using three parameters and a simple explicit equation describing a single fouling mechanism.

3.2.2 BSA filtered through hydrophobic PVDF

The data from the filtration of 1, 2, 4, or 8 g/L 0.1 μm BSA through hydrophobic 0.22 μm GVHP membrane (Figure 3A from Duclos-Orsello et al. (2006)) was fit using Equation 12. The new model is a simple explicit equation that provided a good fit of the four data curves effectively using two fitted parameters: z = 4, and K = 1.21 × 10⁻³ ± 1.1 × 10⁻⁵, as the R ² value was 0.99, as shown in Figure 3. The SSR and RMSE were 0.0286 and 0.0177. The DFE was 91, indicating the model was not overfitting the data. The previous work used at least four parameters to fit the four curves using an integral equation that described three fouling mechanisms occurring sequentially (Duclos-Orsello et al., 2006). Using a value of z of 4 provides fits that are similar but slightly improved over those of the Standard (z = 3, SSR = 0.039) or Intermediate Blocking (z = 5, SSR = 0.0723) models. In this way, the new model provides a rapid way to evaluate if the best fit is provided by any of the four classical fouling models or other models with higher or lower reaction orders z. A value of x (reaction order with respect to foulant concentration) of one was used, as using higher or lower values did not improve the fit. This aligns with the linear dependence of the blocking constant K on foulant concentration assumed in the four classical fouling models.

3.2.3 BSA filtered through hydrophilic PVDF

The data from the filtration of 1, 2, 4, or 8 g/L BSA through hydrophilic 0.22 μm GVWP membrane (Figure 5A from Duclos-Orsello et al. (2006)) was fit using Equation 12. This data set is an interesting test case because it exhibited a negative fouling index. The new model is a simple explicit equation that provided a good fit of the data effectively using two fitted parameters, z = 11 and K = 9.72 × 10⁻⁴ ± 7.0 × 10⁻⁶, as the R ² value was 1.00, as shown in Figure 4. The SSR and RMSE were 0.00747 and 0.00927. The DFE was 87, indicating the model was not overfitting the data. The high reaction order with respect to pore area (z = 11) is not physically realistic and is therefore considered a fitting parameter. The previous work used four parameters to fit the four curves using an integral equation that described three fouling mechanisms occurring sequentially (Duclos-Orsello et al., 2006). The fit was similar but slightly worse when setting the parameters to be equivalent to Cake Fouling (z = 9, K = 8.55 × 10⁻⁴ ± 9.3 × 10⁻⁶, SSR = 0.0192, R ² = 0.99). In this way, the new model provides a rapid way to evaluate if the best fit is provided by any of the four classical fouling models or other models with higher or lower reaction orders z. Again, a value of x (reaction order with respect to foulant concentration) of one was used, as using higher or lower values did not improve the fit. This aligns with the linear dependence of the blocking constant K on foulant concentration assumed in the four classical fouling models.

It is surprising that the hydrophilic PVDF, which had a lower flux decline, would exhibit a higher reaction order with respect to pore area (z = 11 in Figure 4) than the hydrophobic PVDF (z = 4 in Figure 3). For both cases the membrane pore size was 0.22 μm, the pH was 7.0 where the BSA will have a negative charge, and the best fit reaction order with respect to foulant concentration x was one. It is unlikely that the observed difference in the value of z was due to caking above the membranes, as the best-fit value of z for the case of prefiltered BSA on hydrophobic PVDF (z = 3 for Figure 4A in Duclos-Orsello et al. (2006), data not shown) was similar to the value without pre-filtration (z = 4 in Figure 3).

It was also surprising that fouling of hydrophilic 0.22 μm PVDF with BSA exhibited a higher reaction order with respect to pore area (z = 11 in Figure 4) than the fouling of hydrophilic 0.22 μm PVDF by PEG (z = 1.46 in Figure 2). The membranes are similar and, in both cases, the foulants are hydrophilic and should have limited binding to the membranes.

Given the unrealistically high value of z used to fit the data in Figure 4, it is possible that the membrane is composed of a mix of pores that foul rapidly and pores that foul more slowly. Therefore, the flux data were fit with the two-pore model provided in Equation 21, as shown in Figure 5. This model provided a good fit using four fitted parameters: z = 3, K _*a* = 1.82 × 10⁻³ ± 9.4 × 10⁻⁵, K _*b* = 9.44 × 10⁻⁵ ± 1.4 × 10⁻⁵, f _*a* = 0.525 ± 0.021. The R ² value was 1.00. The DFE was 85, indicating the model was not overfitting the data. The SSR and RMSE (0.00766 and 0.00950) were slightly worse than the fits in Figure 4 (0.00747 and 0.00927), despite the use of two additional fitted parameters. However, the parameters used in the dual pore model were more physically realistic (z = 3) than the parameters used in Figure 4 (z = 11).

The fouling data may be affected by factors such as membrane pore size, surface chemistry, and the heterogeneity of foulants. The initial flux data may be influenced by membrane wetting, membrane swelling or compaction. A more detailed physical characterization of the membranes, foulants and their interactions would be informative for elucidating the causes of these observed differences in flux decline curves and to validate the model assumptions.

4 Conclusion

Traditionally, complex membrane fouling data have been fitted using models that assume up to three mechanisms occur sequentially using at least four fitted parameters. In this work a model of membrane fouling was derived that assumes foulant adsorption occurs due to a reaction between foulant and membrane pore surface area, reducing pore radius and increasing membrane resistance. Explicit equations were developed to predict flux, resistance or volume versus time for different values of reaction orders. The model was extended to treat membranes composed of multiple pore radii allowing flow in parallel. Surprisingly, the model reduces to the four classical fouling models using specific reaction orders with respect to membrane area, but with a new capability for fitting a non-linear dependence on foulant concentration. The new model is therefore useful for determining if the optimal data fit will be provided by one of the classical fouling models, or if a better fit will be provided using other reaction orders.

Using data for the fouling of PVDF membranes with PEG or BSA, the new model demonstrated the ability to fit a wider range of data than the four classical models, fit data with a negative fouling index, and fit data with a non-linear dependence on foulant concentration using only two to four fitted parameters. Data with a negative fouling index was fit either with a single pore model using an unrealistically high value of z, or with a dual pore size model using more physically realistic constants. A more detailed characterization of the relationship between the fitted parameter values and the physical characteristics of the membrane (pore size, pore distribution, surface properties), foulant (size, shape, and surface properties) and their interactions would be useful for further validating the modeling assumptions.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Nataraj Sanna Kotrappanavar, Jain University, India

Reviewed by: Gyorgy Szekely, King Abdullah University of Science and Technology, Saudi Arabia

Yun Teng, Singapore University of Social Sciences, Singapore

Footnotes

C, foulant concentration (kg/m³); f, fraction of initial flow going through pores of a specific radius; J, flux (m/s); K, adsorption blocking constant; L, membrane thickness (m); n, number of pores of a specific radius; P, trans-membrane pressure drop (kg/ms²); r, pore radius (m); S, pore surface area (m²); t, time (s); V, volume filtered (m³/m²); V _*F*, volume of foulant deposited by reacting with pore walls (m³); x, reaction orders in terms of foulant concentration C; z, reaction order in terms of membrane pore area S.