Dynamic Modeling of Fouling in Reverse Osmosis Membranes

Abstract

During reverse osmosis (RO) membrane filtration, performance is dramatically affected by fouling, which concurrently decreases the permeate flux while increasing the energy required to operate the system. Comprehensive design and optimization of RO systems are best served by an understanding of the coupling between membrane shape, local flow field, and fouling; however, current studies focus exclusively on simplified steady-state models that ignore the dynamic coupling between fluid flow, solute transport, and foulant accumulation. We developed a customized solver (SUMs: Stanford University Membrane Solver) under the open source finite volume simulator OpenFOAM to solve transient Navier–Stokes, advection–diffusion, and adsorption–desorption equations for foulant accumulation. We implemented two permeate flux reduction models at the membrane boundary: the resistance-in-series (RIS) model and the effective-pressure-drop (EPD) model. The two models were validated against filtration experiments by comparing the equilibrium flux, pressure drop, and fouling pattern on the membrane. Both models not only predict macroscopic quantities (e.g., permeate flux and pressure drop) but also the fouling pattern developed on the membrane, with a good match with experimental results. Furthermore, the models capture the temporal evolution of foulant accumulation and its coupling with flux reduction.

1. Introduction

Reverse osmosis (RO) filtration systems are widely applied in seawater desalination [1,2,3,4,5], landfill leachate treatment [6,7], and wastewater reclamation [8,9,10,11,12,13,14]. Typically, RO performs one of the final stages of water treatment and is designed to remove ions or soluble substances.

Due to the extraction of the solvent (e.g., water) on the feed side close to the membrane surface, the solute concentration rises, which is known as concentration polarization (CP) [15,16,17,18,19]. Some solutes can precipitate or crystallize on the membrane surface, while other solutes adsorb to the membrane, hindering permeation of the solvent and reducing the efficiency of the membrane [20,21]. Such fouling processes cause reduction of clean water permeate flux. By increasing the applied pressure, one can increase the pressure gradient across the membrane to force a larger permeate flux, but the energy input per unit flux increases as a result. Fouling depends on the solute and membrane properties; for instance, biologically active foulants can produce thick, relatively low permeability biofilms [22]. RO membrane modules require spacers to separate membrane leaves and create flow channels; these spacers play an important role in fouling development. For example, the most commonly used net-like spacers create dead zones where foulant cake growth is accentuated [23].

Different mechanisms can be utilized at different scales to control fouling. These include: (i) changes in the morphology (shape) of the flow channel at the system scale (∼cm) [24,25,26,27,28]; (ii) modifications to the topology of the membrane surface at the micro-scale (∼mm–μm) [29,30,31,32,33,34,35]; (iii) chemical or surface treatment which changes the interaction force between foulant and membrane at the nano-scale (∼nm) [36,37,38].

It has been shown that morphological changes can provide in-situ fouling mitigation; a number of studies [35,39,40,41] have demonstrated that flow and solute transport at both the macro- and micro-scales can be controlled by modifying the membrane/spacer morphology. However, most analyses still optimize the system by trial and error since a general framework to study foulant deposition and in situ control is still not available. Due to experimental difficulties and cost, performing extensive studies on different configurations is challenging.

Computational models [5,25,42,43] represent an attractive alternative to more expensive experimentation as they allow one to virtually span the entire design space at a fraction of the cost; however, foulant dynamical behavior is elusive for most existing models. The major challenges associated with modeling dynamic fouling processes are (i) the temporal evolution of foulant deposition and (ii) the strong coupling between flow, bulk solute concentration, and foulant deposition. Complex spacer geometry complicates the matter even further, and, while inherently essential to RO system optimization, modeling the spatio-temporal evolution of fouling remains an open challenge.

Most of the models that account for the temporal evolution of the foulant layer do so without a full coupling between flow, transport, and foulant deposition. For example, Bucs et al. [22] model the thickness of foulant as an empirically postulated function of time with a constant growth rate, yet the velocity and bulk concentration fields are determined from steady-state equations. Xie et al. [44] model fouling accumulation using a temporal adsorption/desorption equation under the hypothesis that the adsorption rate depends on the local bulk concentration. The authors also introduce in their model the process of mechanical removal of foulant due to hydrodynamic shear by introducing a stress term into the growth equation for the foulant. Again, not only are the flow and concentration fields solved by steady-state equations, but the flow field is imposed as a background field without accounting for the feedback from fouling processes. Lyster and Cohen [45] propose a set of equations and boundary conditions that couple the velocity component orthogonal to the membrane surface with the local concentration gradient on the membrane surface. While these conditions, derived by mass balance in two dimensions, are shown to successfully capture concentration polarization (CP) and the coupling between CP, flow and bulk concentration, the model does not account for unsteady terms and does not include a mechanism to relate CP to fouling.

Recently, Ling and Battiato [46] developed a model that couples the transient Navier–Stokes and the advection–diffusion equations, as well as an adsorption–desorption equation for foulant accumulation. Although they validate it against experimental data and demonstrate that it is able to correctly capture unsteady measurements of permeate flux, its capability of correctly capturing spatial distribution of the foulant in morphologically complex membranes was not evaluated. This is a critical step in assessing the potential of using the model as a virtual laboratory for design and membrane performance optimization purposes. It is worth noticing that Ling and Battiato used an effective-pressure-drop (EPD) model, which couples the flux reduction and fouling accumulation by introducing an additional pressure reducing term. The EPD model varies from the more classical approach of treating the foulant layer as an additional flow resistance, which is often referred to as a resistance-in-series (RIS) model [47]. In this study, we use both approaches and compare them.

Here, all processes are modeled using a 3D fully-coupled system of transient equations: the Navier–Stokes equations for flow, an advection–diffusion equation for the bulk concentration, and an adsorption–desorption equation for fouling. Furthermore, the model allows one to relate concentration polarization, occurring in the bulk solution, with fouling taking place on the membrane modeled as a surface concentration. The flux reduction induced by foulant accumulation is modeled using an adsorption–desorption equation which associates the local bulk concentration, foulant surface concentration, and permeate flux. All equations are implemented through a customized solver SUMS (Stanford University Membrane solver) in the open-source finite-volume framework OpenFoam.

The model is validated by comparing three-dimensional simulations with fouling experiments conducted by Xie and et al. [44], who measured (i) the permeate flux and pressure drop and (ii) the spatial distribution of fouling patterns for different spacer configurations. Such comparisons demonstrate the RIS and EPD models’ capability of capturing both system-scale quantities (i.e., flux and pressure) and local effects (fouling pattern). The spacers studied by Xie et al. do not have conventional geometry; they comprise a set of sinusoidal flow channels that vary in amplitude and frequency. This experimental data set, with its unique design, has a wider variation in geometry (and thus a wider range of flow patterns) and spatial scales than most spacer studies; in addition, the data were readily available to us in raw form, making this a useful data set for testing the effectiveness of the SUMS framework.

The paper is organized with Section 2 introducing the governing equations and simulation scenarios. In Section 3, we present the experimental setup and data post-processing technique to digitize the images of fouling patterns. In Section 4, we compare the simulated permeate flux, pressure drop, and fouling pattern with the corresponding experimental results. We provide concluding remarks in Section 5.

2. Materials and Methods

2.1. Formulation

We are interested in studying fouling accumulation on a flat sheet membrane as a function of time, $\widehat{T}$, and location, $\widehat{\mathbf{X}}=(\widehat{X},\widehat{Y},\widehat{Z})$. The flow field ($\widehat{\mathbf{U}}$) of an incompressible viscous fluid satisfies the Navier–Stokes and continuity equations

$$\begin{array}{cc}\hfill & \frac{\partial \widehat{\mathbf{U}}}{\partial \widehat{T}}+(\widehat{\mathbf{U}}·\widehat{\nabla })\widehat{\mathbf{U}}+\frac{1}{\rho }\widehat{\nabla }\widehat{P}=\widehat{\nabla }·\left(\nu \widehat{\nabla }\widehat{\mathbf{U}}\right)\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill & \widehat{\nabla }·\widehat{\mathbf{U}}=0\hfill \end{array}$$

(1b)

where $\widehat{\mathbf{U}}$ [m/s] is a three-dimensional velocity field $\widehat{\mathbf{U}}=(\widehat{U},\widehat{V},\widehat{W})$, $\widehat{P}$ [kg m⁻¹s⁻²] is the pressure, $\rho$ [kg/m³] is the fluid density, and $\nu$ [m²/s] is the fluid kinematic viscosity. Gravity is neglected in this study. The solute bulk concentration satisfies an advection–diffusion equation

$$\begin{array}{cc}\hfill & \frac{\partial {\widehat{C}}_{\mathrm{b}}}{\partial \widehat{T}}+\widehat{\mathbf{U}}·\widehat{\nabla }{\widehat{C}}_{\mathrm{b}}-D{\widehat{\nabla }}^2{\widehat{C}}_{\mathrm{b}}=0,\hfill \end{array}$$

(2)

where ${\widehat{C}}_{\mathrm{b}}(\widehat{\mathbf{X}},\widehat{T})$ is the solute bulk concentration [mol/m³] in the liquid domain, and D [m²/s] is the molecular diffusion coefficient of the solute in water. A Langmuir adsorption–desorption equation (defined on the membrane surface) is used to model foulant accumulation on the membrane located at $\widehat{Z}=H$ (see Figure 1, i.e., the surface concentration of the foulant ${\widehat{C}}_{\mathrm{s}}$ [mol/m²], defined at $\widehat{Z}=H$, satisfies

$$\begin{array}{c}\hfill \frac{\partial {\widehat{C}}_{\mathrm{s}}}{\partial \widehat{T}}=K_1({\widehat{C}}_{\mathrm{s},\max }-{\widehat{C}}_{\mathrm{s}})\widehat{C}-K_2{\widehat{C}}_{\mathrm{s}},\end{array}$$

(3)

where $K_1$ [$1/(\mathrm{mol}·\mathrm{s})$] is the adsorption coefficient, $K_2$ [$1/\mathrm{s}$] is the desorption coefficient and ${\widehat{C}}_{\mathrm{s},\max }$ is the equilibrium foulant concentration. The adsorption model uses the liquid–domain concentration adjacent to the membrane, ${\widehat{C}}_{\mathrm{b}}$, to determine the driving force for foulant adsorption on the membrane. It is worth emphasizing that the same kinetic equation has been adopted in both organic foulant growth [44,48] and crystal growth [49] modeling, where $K_1$ and $K_2$ can be determined via experiments. Additionally, such a framework allows one to evaluate concentration polarization and foulant accumulation individually. Ion (e.g., Ca${}^{2+}$) transport in solution is modeled by ${\widehat{C}}_{\mathrm{b}}$ and its crystallization (e.g., CaSO${}_4$ or CaCO${}_3$) and accumulation on the membrane is modeled by ${\widehat{C}}_{\mathrm{s}}$.

The previous equations are supported by appropriate boundary conditions at the inlet and outlet for the momentum and mass transport problems. Specifically,

$$\begin{array}{c}\hfill \widehat{U}\left(\widehat{X}=0,\widehat{Y},\widehat{Z}\right)=U_{\mathrm{in}}\mathrm{and}\widehat{C}\left(\widehat{X}=0,\widehat{Y},\widehat{Z}\right)=C_0\end{array}$$

(4)

On the solid walls of the channel, no-slip and no-penetration conditions are employed. On the membrane surface, the velocity components $\widehat{U}$ and $\widehat{V}$ are modeled by the Beavers–Joseph condition [50],

$$\begin{array}{cc}\hfill \widehat{U}=& \frac{\sqrt{{\widehat{K}}_{\mathrm{m}}}}{\widehat{\beta }}\frac{\partial \widehat{U}}{\partial \widehat{Z}},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \widehat{V}=& \frac{\sqrt{{\widehat{K}}_{\mathrm{m}}}}{\widehat{\beta }}\frac{\partial \widehat{V}}{\partial \widehat{Z}},\hfill \end{array}$$

(6)

where $\widehat{\beta }$ is a constant that only depends on the geometry of the membrane porous structure. In addition, the flux balancing boundary condition proposed by Lyster and Cohen [45]

$$\frac{\partial {\widehat{C}}_{\mathrm{b}}}{\partial \widehat{Z}}=-\frac{R_{\mathrm{i}}}{D}\widehat{W}{\widehat{C}}_{\mathrm{b}}$$

(7)

is employed. In (7), $\widehat{W}$ is the permeate water flux, $R_{\mathrm{i}}$ is the intrinsic membrane rejection rate [45], (set to $R_{\mathrm{i}}=100\%$ in this study). The permeate flux across a clean membrane is modeled as:

$$\widehat{W}={\widehat{K}}_m\left(\Delta \widehat{P}-\Delta \widehat{\mathsf{\Pi }}\right),$$

(8)

where ${\widehat{K}}_{\mathrm{m}}$ is the hydraulic membrane water permeability in the absence of fouling (i.e., when ${\widehat{C}}_{\mathrm{s}}=0$), the pressure drop $\Delta \widehat{P}$ is defined as $\Delta \widehat{P}=\widehat{P}-P_{\mathrm{amb}}$, with $\widehat{P}$ the local pressure and $P_{\mathrm{amb}}$ the ambient pressure, here set to zero. $\Delta \widehat{\mathsf{\Pi }}$ is the osmotic pressure difference between the feed and permeate, here we assume concentration at the permeate side is zero, namely:

$$\widehat{W}={\widehat{K}}_{\mathrm{m}}\left(\Delta \widehat{P}-{\widehat{A}}_{\mathrm{o}}{\widehat{C}}_{\mathrm{b}}{|}_{Z=H}\right),$$

(9)

where ${\widehat{A}}_o$ [${\mathrm{m}}^2/(\mathrm{s}·\mathrm{mol})$] is the osmotic coefficient and ${\widehat{C}}_{\mathrm{b}}{|}_{Z=H}$ is the bulk concentration near the membrane surface. When the local concentration increases, the permeate flux decreases due to the osmotic pressure. Additional flux reduction due to fouling can be modeled through (i) a resistance-in-series (RIS) model, and (ii) an effective pressure drop (EPD) model, which are discussed in the following.

2.2. Resistance-in-Series Model

The RIS model treats the foulant layer and the membrane as flow resistors that connect in series such that the fouled membrane permeability is ${\widehat{K}}_{\mathrm{eff}}\left({\widehat{C}}_{\mathrm{s}}\right)$ and is modeled as

$${\widehat{K}}_{\mathrm{eff}}={\displaystyle \frac{1}{R_{\mathrm{m}}+R_{\mathrm{f}}}}.$$

(10)

The former relationship quantifies the combined resistance induced by the membrane and the accumulated foulant. In (10), $R_{\mathrm{m}}$ is the clean membrane resistance,

$$R_{\mathrm{m}}=\frac{1}{{\widehat{K}}_{\mathrm{m}}},$$

(11)

and $R_{\mathrm{f}}$ is the fouled membrane resistance,

$$R_{\mathrm{f}}=\frac{C_{\mathrm{s}}}{{\widehat{K}}_{\mathrm{f}}},$$

(12)

where ${\widehat{K}}_{\mathrm{f}}$ is the fouled membrane permeability and $C_{\mathrm{s}}$ is the normalized surface concentration: $C_{\mathrm{s}}={\widehat{C}}_{\mathrm{s}}/{\widehat{C}}_{\mathrm{s},\max }$. When $C_{\mathrm{s}}=1$, the foulant layer results in the maximum flow resistance. The foulant permeability ${\widehat{K}}_{\mathrm{f}}$ is modeled as a proportion of the clean membrane permeability, i.e.,

$${\widehat{K}}_{\mathrm{f}}=A_{\mathrm{k}}{C}_{\mathrm{s}}{\widehat{K}}_{\mathrm{m}},$$

(13)

where $A_{\mathrm{k}}=(0,1]$ is a dimensionless constant. Combining (9) with (10), while accounting for (11)–(13), the permeate flux across a fouled membrane in the RIS model can be written as

$${\widehat{W}}_{\mathrm{RIS}}=\frac{\Delta \widehat{P}-{\widehat{A}}_{\mathrm{o}}{\widehat{C}}_{\mathrm{b}}}{R_{\mathrm{m}}\left(1+C_{\mathrm{s}}/A_{\mathrm{k}}\right)}=\frac{A_{\mathrm{k}}{\widehat{K}}_{\mathrm{m}}}{C_{\mathrm{s}}+A_{\mathrm{k}}}\left(\Delta \widehat{P}-{\widehat{A}}_{\mathrm{o}}{\widehat{C}}_{\mathrm{b}}\right).$$

(14)

It is worth emphasizing that, when $C_{\mathrm{s}}=0$, then ${\widehat{K}}_{\mathrm{eff}}={\widehat{K}}_{\mathrm{m}}$, then relationship (9) for clean membranes is recovered. However, the model cannot capture local clogging of the membrane (i.e., $W_{\mathrm{RIS}}=0$) when $C_{\mathrm{s}}=1$, since such condition would require $R_{\mathrm{f}}\to \infty$, or $A_{\mathrm{k}}\left(C_{\mathrm{s}}\right)$ such that $A_{\mathrm{k}}(C_{\mathrm{s}}=1)=0$, which contradicts the model formulation where $A_{\mathrm{k}}$ is just a fitting constant different from zero.

2.3. Effective Pressure Drop Model

In the EPD model, Equation (9) is generalized under fouled conditions through a modification of the effective driving pressure drop, $(\Delta \widehat{P}-{\widehat{A}}_{\mathrm{o}}{\widehat{C}}_{\mathrm{b}})$, where a pressure reduction due to local foulant accumulation, ${\widehat{A}}_{\mathrm{p}}{\widehat{C}}_{\mathrm{s}}$ [46], is introduced,

$${\widehat{W}}_{\mathrm{EPD}}={\widehat{K}}_{\mathrm{m}}\left(\Delta \widehat{P}-{\widehat{A}}_{\mathrm{o}}{\widehat{C}}_{\mathrm{b}}-{\widehat{A}}_{\mathrm{p}}{\widehat{C}}_{\mathrm{s}}\right).$$

(15)

In (15), ${\widehat{A}}_{\mathrm{p}}$ is a foulant coefficient. Equation (9) for clean membranes is readily recovered when ${\widehat{C}}_{\mathrm{s}}=0$ and ${\widehat{C}}_{\mathrm{b}}=0$. In addition, (15) is able to capture local blockage (${\widehat{W}}_{\mathrm{EPD}}=0$) when ${\widehat{A}}_{\mathrm{o}}{\widehat{C}}_{\mathrm{b}}+{\widehat{A}}_{\mathrm{p}}=\Delta \widehat{P}$. Additionally, the EPD formulation (15) directly associates the flux reduction with precipitation kinetics. This allows one to achieve the coupling between flow, bulk transport, and foulant deposition exclusively through boundary conditions on the membrane surface, without the need for additional ad hoc parametrization of the fouled membrane resistance. In this study, we will compare these two approaches.

Once the transient, coupled flow and transport problems are solved by using the RIS or the EPD model for fouling, the permeate flow rate $\widehat{Q}$ $[{\mathrm{m}}^3/\mathrm{s}]$ can be calculated as

$$\widehat{Q}=\frac{{\int }_{{\mathsf{\Gamma }}_{\mathrm{n}}}{\widehat{W}}_{i}dA}{{\int }_{{\mathsf{\Gamma }}_{\mathrm{m}}}dA},i=\{\mathrm{RIS},\mathrm{EPD}\}$$

(16)

where $\widehat{W}$ is defined by either (14) or (15), respectively, and ${\mathsf{\Gamma }}_{\mathrm{n}}$ is the non-fouled region of the membrane surface and is defined by using a threshold value of the surface concentration $C_{\mathrm{s}}$, i.e., $\alpha C_{\mathrm{s},\max }$ (with $\alpha =0.7$ in this study), as

$${\mathsf{\Gamma }}_{\mathrm{n}}\in \left\{{\mathsf{\Gamma }}_{\mathrm{m}}|C_{\mathrm{s}}\ge \alpha C_{\mathrm{s},\max },\alpha \in [0,1]\right\}.$$

(17)

The set of Equations (1)–(17) can be cast in dimensionless form. We define the dimensionless quantities

$$\begin{array}{c}\hfill \mathbf{u}=\frac{\widehat{\mathbf{U}}}{U_{\mathrm{in}}},\mathbf{x}=\frac{\widehat{\mathbf{X}}}{B},t=\frac{U_{\mathrm{in}}\widehat{T}}{B},P=\frac{B^2\widehat{P}}{{\nu }^2},h=\frac{H}{B},C_{\mathrm{s}}=\frac{{\widehat{C}}_{\mathrm{s}}}{C_{\mathrm{s},\max }},C_{\mathrm{b}}=\frac{\widehat{C}}{C_0},\end{array}$$

(18)

where $\mathbf{u}=(u,v,w)$ and $\mathbf{x}=(x,y,z)$ are the dimensionless velocity field and coordinate axes, respectively. We also introduce the following dimensionless numbers,

$$\begin{array}{cc}\hfill & \mathit{Re}=\frac{U_{\mathrm{in}}B}{\nu },\mathit{Pe}=\frac{U_{\mathrm{in}}B}{D},{\mathit{Da}}_{\mathrm{I}}=K_1\frac{B{C}_0}{U_{\mathrm{in}}},{\mathit{Da}}_{\mathrm{II}}=K_2\frac{B}{U_{\mathrm{in}}},\hfill \end{array}$$

(19)

where $\mathit{Re}$, $\mathit{Pe}$, ${\mathit{Da}}_i$, $i=\left\{\mathrm{I},\mathrm{II}\right\}$ are the Reynolds, Péclet and Damköhler numbers, respectively. Then, the dimensionless form of Equations (1)–(3) reads as follows:

$$\begin{array}{cc}\hfill & \frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}·\nabla )\mathbf{u}+\nabla P=\frac{1}{\mathit{Re}}{\nabla }^2\mathbf{u},\hfill \end{array}$$

(20a)

$$\begin{array}{cc}\hfill & \nabla ·\mathbf{u}=0,\hfill \end{array}$$

(20b)

for flow, and

$$\begin{array}{cc}\hfill & \mathit{Pe}\left(\frac{\partial C_{\mathrm{b}}}{\partial t}+\mathbf{u}·\nabla C_{\mathrm{b}}\right)-{\nabla }^2{C}_{\mathrm{b}}=0,\hfill \end{array}$$

(21a)

$$\begin{array}{cc}\hfill & \frac{\partial C_{\mathrm{s}}}{\partial t}={\mathit{Da}}_{\mathrm{I}}(1-C_{\mathrm{s}})C-{\mathit{Da}}_{\mathrm{II}}{C}_{\mathrm{s}},\hfill \end{array}$$

(21b)

for transport. On the membrane surface ($z=h$), the dimensionless boundary conditions for flow and transport are slip conditions in the direction parallel to the membrane

$$\begin{array}{cc}\hfill u_{\mathrm{h}}=& \frac{\sqrt{k_{\mathrm{m}}}}{\beta }\frac{\partial u_{\mathrm{h}}}{\partial \mathrm{z}}\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill v_{\mathrm{h}}=& \frac{\sqrt{k_{\mathrm{m}}}}{\beta }\frac{\partial v_{\mathrm{h}}}{\partial \mathrm{z}}\hfill \end{array}$$

(23)

where $\beta$ is the Beavers–Joseph constant and is selected to be 2, and the dimensionless flux balancing condition for mass transport

$$\begin{array}{c}\hfill \frac{\partial C_{\mathrm{b}}}{\partial z}=\mathit{Pe}{w}_{\mathrm{i}}{C}_{\mathrm{b}},i=\left\{\mathrm{RIS},\mathrm{EPD}\right\}\end{array}$$

(24)

where

$$\begin{array}{c}\hfill A_{\mathrm{o}}=\frac{{\widehat{A}}_{\mathrm{o}}{B}^2}{{\nu }^2},w_{\mathrm{h}}=\frac{W_{\mathrm{H}}}{U_{\mathrm{in}}},\Delta P=\frac{\Delta \widehat{P}{B}^2}{{\nu }^2},\end{array}$$

(25)

and

$$\begin{array}{c}\hfill w_{\mathrm{RIS}}=\frac{k_{\mathrm{m}}{A}_{\mathrm{k}}}{\mathit{Re}\left(C_{\mathrm{s}}+A_{\mathrm{k}}\right)}\left(\Delta P-A_{\mathrm{o}}{C}_{\mathrm{b}}\right)\end{array}$$

(26)

for the resistance-in-series model, or

$$\begin{array}{c}\hfill w_{\mathrm{EPD}}=\frac{k_{\mathrm{m}}}{\mathit{Re}}\left(\Delta P-A_{\mathrm{o}}{C}_{\mathrm{b}}-A_{\mathrm{p}}{C}_{\mathrm{s}}\right)\end{array}$$

(27)

for the effective pressure drop model. In (26) and (27), the dimensionless permeability $k_{\mathrm{m}}$ is defined as

$$\begin{array}{c}\hfill k_{\mathrm{m}}=\frac{{\widehat{K}}_{\mathrm{m}}}{B\nu },k_{\mathrm{f}}=\frac{{\widehat{K}}_{\mathrm{f}}}{B\nu }.\end{array}$$

(28)

A complete list of all boundary conditions is provided in

Table 1. Boundary Conditions for the simulation.

Boundary	Flow	Bulk Concentration	Pressure
	$\mathbf{u}$	${\mathbf{C}}_{\mathbf{b}}$	$\mathbf{P}$
Inlet, ${\mathsf{\Gamma }}_{\mathrm{i}}$	$\mathbf{u}=(1,0,0)$	$C_{\mathrm{b}}=1$	$\partial P/\partial \mathbf{n}=0$
Outlet, ${\mathsf{\Gamma }}_o$	$\partial \mathbf{u}/\partial \mathbf{n}=0$	$\partial C_{\mathrm{b}}/\partial \mathbf{n}=0$	$\partial P/\partial \mathbf{n}=0$
Solid Wall, ${\mathsf{\Gamma }}_{\mathrm{w}}$	$\mathbf{u}=0$	$\partial C_{\mathrm{b}}/\partial \mathbf{n}=0$	$\partial P/\partial \mathbf{n}=0$
Membrane, ${\mathsf{\Gamma }}_{\mathrm{m}}$	$\mathbf{u}=(u_{\mathrm{h}},v_{\mathrm{h}},w_{\mathrm{h}})$	$\partial C_{\mathrm{b}}/\partial \mathbf{n}=w_{\mathrm{h}}{C}_{\mathrm{b}}$	${\widehat{P}}_{\mathrm{out}}/\left({\nu }^2/B^2\right)$

. The dimensionless permeate flow rate is

$$q=\frac{{\int }_{{\mathsf{\Gamma }}_{\mathrm{n}}}{w}_{\mathrm{i}}dA}{{\int }_{{\mathsf{\Gamma }}_{\mathrm{m}}}dA},i=\left\{\mathrm{RIS},\mathrm{EPD}\right\}.$$

(29)

The 3D model (20)–(29) is implemented through the customized solver SUMs (Stanford University Membrane solver) in the open source finite volume simulator OpenFOAM, where an implicit time scheme for the transient solver and second order discretization in space are employed. The numerical mesh of the simulation is generated by a built-in OpenFOAM mesh tool, SnappyHexMesh, and the mesh resolution is determined such that the thinnest throat in the channel contains 15 numerical grids.

3. Experimental Data and Image Post-Processing

Experiments were performed with spacers inserted into a flat-sheet crossflow test cell. Each spacer formed ten equivalent flow paths on the membrane, see Figure 2. Each flow path was 6 mm wide and 1.5 mm high. The membrane was on the 6 mm side of the flow path. The straight-line distance between the entrance and exit of each flow path was 130 mm, resulting in an active membrane area of 780 mm² for all configurations.

The experiments involved four sinusoidal spacers with different amplitudes and periods and a straight channel membrane for benchmark, see Figure 3. The data collected involve measurements of steady-state permeate flux and pressure drop [44], as well as spatial distribution of the foulant on the membrane surface after flooding $1L$ concentrated solution. The full description of the setup, data collection procedure, and data type can be found in [24,44]. A list of experimental parameters is provided in

Table 2. Experimental parameters, symbols, values, and corresponding units.

Parameter	Experimental Parameters
	Symbol	Value	$[\mathbf{L},\mathbf{T},\mathbf{M}]$
Width	B	$6\times {10}^{-3}$	m
Height	H	$1.5\times {10}^{-3}$	m
Viscosity	$\nu$	$1\times {10}^{-6}$	${\mathrm{m}}^2/s$
Diffusion Coeff.	D	$2\times {10}^{-11}$	${\mathrm{m}}^2/\mathrm{s}$
Outlet Pressure	$P_{\mathrm{out}}$	4137	Psi
Inlet Velocity	$U_{\mathrm{in}}$	$0.15$	$\mathrm{m}/\mathrm{s}$
Concentration	$C_0$	50	$\mathrm{mmol}/\mathrm{L}$
Permeability	$K_{\mathrm{m}}$	$5\times {10}^{-8}$	$\mathrm{m}/(\mathrm{s}·\mathrm{kPa})$

. The experimental data collected include measured permeate flux, pressure drop, and fouling pattern on the membrane surface.

Images of fouling patterns for different spacer morphologies need to be processed to map color intensity into surface concentration for comparison with numerical simulation. This is achieved in three sequential steps: (i) one flow channel is extracted from the raw image of the membrane, (ii) the image color intensity (in gray scale) is mapped to surface concentration according to

$$\begin{array}{c}\hfill C_{\mathrm{s},\exp }=C_{\mathrm{s},\max }·\frac{I}{I_{\max }},\end{array}$$

(30)

where $C_{\mathrm{s},\exp }$ is the surface concentration from the experiment, $I_{\max }$ is the maximum gray scale intensity and I is the gray scale intensity at a given location; (iii) the experimental fouling patterns for the different spacers morphologies are obtained by thresholding the surface concentration as specified in Equation (17), i.e., surface concentration equal to or higher than the threshold value $\alpha C_{\mathrm{s},\max }$ (with $\alpha =0.7$) is used to represent the experimental fouling pattern. In Figure 4, we show the unprocessed pictures of the fouling patterns in a single channel (top) and the fouling patterns after mapping to concentration fields (bottom) for each spacer morphology. The latter are used for a direct comparison with numerically simulated fouling patterns as discussed in the following section.

4. Results and Discussion

In this section, we present the simulation results from the two fouling models, the RIS and the EPD, defined by Equations (26) and (27), respectively. Both models are used to predict fouling, steady state permeate flux, and pressure drop for all five geometries.

The simulation parameters are set equal to the values reported in the experiments [24], and listed in Table 2. Additionally, studies on membrane adsorption/desorption rates have shown that the ratio between $K_2$ and $K_1$ varies from $0.001$ to 1 [48]. In our study, we set $\theta =K_2/K_1=0.1$. More specifically, since the absolute values of $K_1$ and $K_2$ only affect how fast the foulant reaches equilibrium $C_{\mathrm{s},\max }$, we select $K_1=0.1$ and $K_2=\theta K_1=0.01$. The dimensionless number corresponding to the experimental conditions investigated are reported in

Table 3. The fixed parameters of all the simulations.

Model	Re	${\mathit{Da}}_{\mathbf{I}}$	${\mathit{Da}}_{\mathbf{II}}$	${\mathit{A}}_{\mathbf{o}}$	$\mathit{\beta }$	${\mathit{A}}_{\mathbf{k}}$	${\mathit{A}}_{\mathbf{p}}$
RIS EPD	900	$5\times {10}^7$	$4\times {10}^{-3}$	$0.067$	2	0.07 -	- 3600

, where $\mathit{Re}$ is determined by the experimental parameters, and ${\mathit{Da}}_{\mathrm{I}}$ and ${\mathit{Da}}_{\mathrm{II}}$ are determined by the selection of $K_1$ and $K_2$. Furthermore, $A_{\mathrm{o}}$ is fitted by using the experimental data of R1. We note that, in addition to the parameters listed above, which are shared by both the RIS and EPD models, each model has one undetermined parameter: $A_{\mathrm{k}}$ in the RIS model, and $A_{\mathrm{p}}$ in the EPD model. Such parameters are fitted from experimental flux measurements on the benchmark rectangular geometry, R1, and then kept constant to predict flux, pressure, and fouling pattern for the all other geometries with $A_{\mathrm{k}}=0.067$ and $A_{\mathrm{p}}=3600$, for the RIS and EPD models, respectively. In each simulation, the inlet concentration is set to $C_{\mathrm{b}}=1$ when $t=t_0$, i.e.,

$$C_{\mathrm{b}}{|}_{\mathbf{x}\in {\mathsf{\Gamma }}_{\mathrm{i}}}=\left\{\begin{array}{cc}0\hfill & t<t_0,\hfill \\ 1\hfill & t⩾t_0.\hfill \end{array}\right.$$

(31)

For all simulations $t_0=120\left[\mathrm{s}\right]$ and the total simulated time is 0.5 h.

4.1. Steady-State Flux and Pressure

We compare the steady-state permeate flux measured in the experiments with the simulated flux. The flux is numerically computed from Equation (16) with $\widehat{W}$ defined by (14) or (15) for the RIS and EPD models, respectively. The pressure drop, $\Delta {\widehat{P}}_{\mathrm{L}}$, is calculated as

$$\Delta {\widehat{P}}_{\mathrm{L}}=\frac{{\widehat{P}}_{\mathrm{in}}-{\widehat{P}}_{\mathrm{out}}}{L},$$

(32)

where ${\widehat{P}}_{\mathrm{in}}$ is the average pressure at the inlet, ${\widehat{P}}_{\mathrm{out}}$ is the imposed pressure at the outlet, and L is the total length of the channel. In Figure 5, we plot both the measured and the simulated permeate fluxes from the RIS and EPD models, as well as the pressure drop for all five spacer configurations. The EPD model exhibits better agreement with the experimental flux than the RIS model, which underestimates the experimental flux when the geometry is more torturous (S1–S4). The difference in performance between the two models can be explained as follows. When the geometry is more torturous, the foulant distribution exhibits a more heterogeneous pattern along the membrane (see Figure 4), and it is associated with a less uniform velocity distribution. In the RIS model, the flow resistance is modeled by the combination of the membrane resistance and the foulant resistance, where the membrane permeability is a small value: as a result, the local flux is less sensitive to foulant distribution heterogeneity. On the other hand, in the EPD model, effective pressure loss due to the foulant is calculated by using a linear dependence which allows for accounting for the direct impact of local foulant variations on flux. Both models provided similar longitudinal pressure drop prediction, and the results are in good agreement with the experimental results, except for the S4 geometry. A good match with experiments is expected since the flux in RO systems is small compared to the crossflow velocity and thus flux boundary conditions would not significantly alter the longitudinal flow. Overall, the results suggest that, for the straight spacer (R1), both the RIS and EPD models can match the experimental results, while, for sinusoidal spacers (S1–S4), the EPD model can more accurately predict the measured flux. To estimate the overall accuracy of each model, we define error associated with the prediction of the permeate flux,

$$\begin{array}{c}\hfill Er{r}_{\mathrm{f}}=\frac{\sqrt{\mathsf{\Sigma }{\left[(q_{\mathrm{i},\mathrm{sim}}-q_{\mathrm{i},\exp })/q_{\mathrm{i},\exp }\right]}^2}}{N}i=\{R1,S1–S4\},\end{array}$$

(33)

where $q_{\mathrm{i},\mathrm{sim}}$ and $q_{\mathrm{i},\exp }$ are the numerical and experiment pressure drop, respectively. We also define the error associated with the pressure drop estimation as:

$$\begin{array}{c}\hfill Er{r}_{\mathrm{p}}=\frac{\sqrt{\mathsf{\Sigma }{\left[(\Delta {\widehat{P}}_{\mathrm{i},\mathrm{sim}}-\Delta {\widehat{P}}_{\mathrm{i},\exp })/\Delta {\widehat{P}}_{\mathrm{i},\exp }\right]}^2}}{N}i=\{R1,S1–S4\},\end{array}$$

(34)

where $\Delta {\widehat{P}}_{\mathrm{i},\mathrm{sim}}$ and $\Delta {\widehat{P}}_{\mathrm{i},\exp }$ are the flux results of simulation and experiment, respectively. The error of the RIS is $Er{r}_{\mathrm{f}}=10.78\%$ and $Er{r}_{\mathrm{p}}=11.93\%$, the error of the EPD model is $Er{r}_{\mathrm{f}}=4.51\%$ and $Er{r}_{\mathrm{p}}=11.99\%$, which is consistent with Figure 5.

4.2. Dynamics

Tracking flux decline is an essential component of assessing the filtration process as the decline curve tracks the correlation between foulant accumulation and flux reduction.

In Figure 6, we plot both the average permeate flux normalized by the flux before solute injection begins,

$$\begin{array}{c}\hfill q^{🟉}={\displaystyle \frac{q}{q(t<t_0)}},\end{array}$$

(35)

as well as the average foulant accumulation defined as

$$\begin{array}{c}\hfill 〈C_s〉=\frac{{\int }_{{\mathsf{\Gamma }}_{\mathrm{m}}}{C}_{\mathrm{s}}dA}{{\int }_{{\mathsf{\Gamma }}_{\mathrm{m}}}dA},\end{array}$$

(36)

for the R1 and the S4 geometries and the RIS and EPD models. Both models show transient flux reduction and the flux results are closely coupled with foulant accumulation: as the foulant builds up, permeate flux decreases. Both models predict similar foulant accumulation, although the EPD model shows faster foulant buildup in the initial stage. For the flux reduction curves, the RIS shows little dependence on the two geometries, while the EPD model is able to better capture flux differences between between R1 and S4.

4.3. Fouling Pattern

Once the models have been validated against device-scale measurements, we proceed to test their ability to reproduce the spatial distribution of fouling patterns at steady state. For the RIS and the EPD models, we select the $\alpha ={\alpha }^{\ast }$ individually to plot the foulant distribution, where ${\alpha }^{\ast }$ is determined such that the foulant coverage reaches 50% of the total area of the membrane, i.e.,

$$\begin{array}{c}\hfill \frac{{\int }_{{\mathsf{\Gamma }}_{\mathrm{n}}}dA}{{\int }_{{\mathsf{\Gamma }}_{\mathrm{m}}}dA}=0.5,\mathrm{where}{\mathsf{\Gamma }}_{\mathrm{n}}\in \left\{{\mathsf{\Gamma }}_{\mathrm{m}}|C_{\mathrm{s}}\ge {\alpha }^{\ast }{C}_{\mathrm{s},\max }\right\}\end{array}$$

(37)

In Figure 7, Figure 8 and Figure 9, we compare the experimental and predicted fouling patterns from the two models for all geometries. Overall, the model results show good agreement with data regardless of the flux boundary condition used. Specifically, the models correctly capture a number of features in the experimental fouling patterns: (i) more foulant accumulates near the outlet than at the inlet, (ii) foulant accumulates at the peaks and troughs of the sinusoidal channel, and (iii) for sinusoidal spacers with larger amplitude, the fouling pattern develops an asymmetric shape with not symmetric tails extending upstream.

Overall, the fouling pattern exhibits strong spatial heterogeneity, a result of coupling between adsorption and local flow conditions, which can significantly differ across channel morphologies. A framework coupling between flow, solute transport, and foulant accumulation is robust in modeling heterogeneous spacers and can accurately predict high fouling zones.

5. Conclusions

In this study, we investigate the ability of two different fouling models (RIS and EPD) to correctly capture both system-scale performance quantities, namely permeate flux and pressure drop, as well as fine-scale features, such as high fouling regions. The two models are constructed as boundary conditions on the membrane surface and implemented in the code SUMs within the OpenFOAM framework. Both fouling models have only one fitting parameter, calibrated against the rectangular membrane benchmark geometry. Fit-free predictions are then performed on four membranes with sinusoidal spacers with different amplitudes and frequencies. Model predictions are tested against the experimental data, which included both system scale measurements (flux decline curves and pressure drop) and local measurements (fouling patterns). Both models were overall able to capture both (i) system-scale pressure drop and (ii) spatio-temporal fouling patterns for five different spacer geometries, although the EPD model was more sensitive to the impact of spacer morphologies on flux, and therefore better able to predict both flux decline and steady-state flux for different morphologies. Both the RIS and the EPD models were successful in capturing the spatial distribution of foulant, and its main experimentally observed features. These results suggest that such a framework is able to successfully (i) simulate flow, transport, and fouling process using transient equations; (ii) couple the flow, bulk concentration and surface concentration of the foulant dynamically while elucidating foulant accumulation mechanisms; (iii) associate concentration polarization with fouling by using an adsorption type equation, and (iv) incorporate different flux reduction models such as the RIS model and the EPD model. Future work includes generalization of the code to Membrane Distillation (MD) processes and other filtration techniques as well as to more complex spacer geometries and system-scale (i.e., module scale) domains. Moreover, as a three-dimensional simulator, SUMs are compatible with three-dimensional geometries imported directly from design tools. As a result, the code can be directly used for filtration systems optimization in industrial applications.