A Finite Element Stiffness Method for Modeling Solute Transport in Reverse Osmosis Membranes

Farah Najdawi

Materials Science, Engineering and Commercialization Program, Texas State University, San Marcos, TX, USA.

DOI: 10.4236/eng.2025.1710029                

Abstract

With increasing freshwater demand and rising water scarcity, reverse osmosis (RO) has become a key method for seawater desalination; however, accurately modeling RO remains challenging due to its complex, nonlinear processes. This study develops a new mass transport model by coupling the Advection-Diffusion Model (ADM) with the Solution-Diffusion Model (SDM), and it is solved by using the Finite Element Method (FEM) with the linear stiffness approach. The extended ADM-SDM framework combines both advective and diffusive transport forces and calculates osmotic pressure. The model is validated against COMSOL simulations and commercial RO module data, finding accurate prediction of solute concentration and permeate flux, with errors as low as 0.7% for SW30-2514 and 2.6% for SW30-4021 compared to commercial RO module data. Unlike the traditional SDM, the extended ADM-SDM framework captures advective transport effects, leading to improved permeate flux estimation while maintaining 100% salt rejection. The numerical stability was further validated by mesh sensitivity analysis, which also demonstrated the robustness and effectiveness of the FEM stiffness implementation by showing that outlet flux errors successfully converged to zero across refinements.

Keywords

Advection-Diffusion Model, Nernst-Planck Model, Finite Element Method, Stiffness Matrix Method, Mesh Sensitivity Analysis, Net Driving Force Pressure

1. Introduction

The mass transport through membranes can be described by various mathematical relations. The transport mathematical models are derived from one of the two main group models: irreversible thermodynamic and hydrodynamic model [1]. The irreversible thermodynamic connects the fluxes of solvent in feed water with driving forces by linear relationship such as solution diffusion model. Also, the equation is independent of membrane structure [2]. Over the years, many different salt rejection irreversible thermodynamic and hydrodynamic models have been created to explore the effective rejection of salt ions in RO systems. The hydrodynamic model refers to mass movement carried along by bulk fluid motion [3]. Thermodynamic transport occurs through diffusion, where molecules dissolve into a medium and migrate along concentration gradients [4]. This study utilizes both models.

Irreversible thermodynamics uses linear relationships, such as the solution diffusion model SDM, to relate solvent fluxes in feed water to driving forces. For instance, [5] shows that the Na⁺/K⁺ pump functions irreversibly by employing adenosine triphosphate (ATP) to drive ions against their gradients, underscoring the rationale for neuronal models to emphasize active transport mechanisms over passive diffusion. [6] used dimensionless variables and algebraic equations to create a simplified framework for standardizing flux prediction in membranes. Their work, based on the solution-diffusion model and compared with the solution-friction model, shows both the strengths and limits of assuming constant transport parameters. [7] demonstrated that substrate hydrophilicity affects the structure and efficacy of thin film composite nanofiltration in reverse osmosis membranes. Their results are consistent with the irreversible thermodynamic solution-diffusion model, which posits that transport is governed by sorption and diffusion rather than hydrodynamic pore flow.

It also considers non-linear connections for irreversible thermodynamics models, like the Nernst-Planck Model (NPM), to accurately describe mass transport. The equation is independent of membrane structure [8]. The NPM is derived from the Maxwell Stephan Equation and can be used to observe how transported water is contaminated as it passes through a membrane in reverse osmosis [9]. Also, it accounts for coupled transport phenomena, including the interactions among water, solutes, and the membrane material [10]. This enables precise predictions of water flux and solute rejection.

The osmotic pressure of the treated water becomes the first measurement derived from SDM. In this study, SDM is used to calculate the total net pressure required for system operation. [11] investigated low salt rejection performance in RO membrane systems using module-scale analysis. In addition to determining the system’s specific energy consumption, they examined the impact of feed concentration on brine concentration and permeability flux. [12] compared the pore flow model and SDM to examine transport mechanisms in RO membranes. In order to investigate the existence of a pure water layer at the sodium chloride-cellulose acetate membrane interface and enhance membrane performance and desalination technology, they reviewed the Gibbs Adsorption Isotherm.

[13] validated NPM, by creating an inverse averaging finite element method (IAFEM) to investigate salt permeability in nanofiltration systems under various feed salt concentrations and pressures. In addition to addressing numerical instability and convergence issues related to convection-dominated NPM and large-scale nanopore systems, the study examined pressure drop distributions across membranes, refuting the assumption of no pressure drop in the solution-diffusion model.

In order to investigate the permeate pH across a RO membrane, Zhang, Hamelers, and [14] used the Donnan steric model to model the transport ion mechanism. The results showed that the permeate pH value relied highly on the membrane charge density and less significant effected by salt concentration and feed flux. Also, the permeate pH level might be lower or higher than feed pH solution and this is determined by the polarity of the membrane charged. [15] simulated the impact of initial pH on salt ion transport during the first two hours of saltwater electrolysis using a validated Nernst-Planck ion transport model in COMSOL. In order to overcome a major obstacle and increase the feasibility of green hydrogen production, they discovered that adjusting the initial pH can significantly reduce the transport of chloride ions between the anolyte and catholyte.

While the Simplified SDM assumes constant diffusivity and permeability, which can actually vary, other researchers like [13] [15] and [16] have investigated ion transport using the NPM. However, our study distinguishes FROM existing literature by concentrating on the numerical implementation and validation of mass transport models. It particularly fills in the gaps left by earlier studies, which often lacked comprehensive validation and infrequently used sophisticated numerical techniques like the linear stiffness approach in the FEM. This FEM-based method offers a balance between high accuracy and stability and less computational complexity than traditional Computational Fluid Dynamics (CFD) models.

In order to address the limitations of the basic SDM, this study develops the solution diffusion model and advection diffusion model (SDM-ADM) by adding the advection term from NPM. This expanded model takes into consideration different ion transport forces in addition to osmotic pressure. The FEM-linear stiffness approach is used to solve the SDM-ADM model, and the results are compared to experimental data and CFD-COMSOL simulations.

Following the introduction and literature review, the remainder of the paper is organized as follows. Section 2 presents the methodology, describing the Extended SDM-ADM model, the application of the Finite Element Method (FEM) to solve ADM, and the use of COMSOL Multiphysics. Section 3 reports the results from the FEM stiffness method, including concentration and flux profiles along the membrane length, the effects of varying boundary condition configurations in ADM, and a mesh sensitivity analysis. Section 4 provides the results from COMSOL Multiphysics simulations, while Section 5 examines solvent behavior in the SDM. Section 6 compares the outcomes with commercial RO modules. Finally, Section 7 presents the conclusions of this work. Table 1 summarizes the notation and unit of the model parameters.

Table 1. Notation description.

Symbol	Definition	Unit
$J_i$	Flux of solute i through the membrane	m3/(s·m2)
$D_i$	Diffusion coefficient of solute i	m2/s
$\nabla c_i$	Concentration gradient of solute i	mol/m3·m
$c_i$	Concentration of solute i	mol/m3
$v$	Flow (advective) velocity	m/s
$C\left(x,t\right)$	Solute concentration	mol/m3
$D_{\text{eff}}$	Effective solute diffusivity	m2/s
$D_0$	Bulk diffusivity	m2/s
$\alpha$	Solute-membrane interaction coefficient	-
K	Stiffness matrix in FEM	-
${\varphi }_i$	test function	-
$\left\{F\right\}$	Force vector (boundary conditions and source terms)	-
L	Element length in FEM mesh	m
Aw	Water permeability coefficient	m/(s·Pa)
Q	Solvent (water) flux	m3/day
ΔP	Net driving pressure across membrane	Pa
RR	Salt rejection	%

2. Methodology

2.1. The Extended SDM-ADM Model

The objective of this section is to solve the extended model which is derived from NPM to determine the concentration profile within the RO membrane and other design parameters of the membrane. In fact, NPM describes the ions motion in RO membrane composed of three forces: diffusion, advection and electromigration. To simplify our analysis and also differentiate our contribution from previous research, we focus solely on the diffusion and the advection as the driving forces of salt ions through the RO membrane, disregarding the impact of electromigration force. This is because for typical RO desalination of seawater, the feed solution is electrically neutral [14] In such a condition the contribution of electromigration force is minimal compared to conviction and advection forces. By neglecting the electromigration force, the NPM can be written as:

$J_i=-D_i\nabla c_i-c_{i}v$ (1)

where, $J_i$ is the flux of ion i in m³/s/m², $D_i$ is the diffusion coefficient of ion i in m²/s, $\nabla c_i$ is the concentration gradient of ion i in mole/m³, $c_i$ the concentration of ion i in mole/m³, and $v$ the flow velocity in m/s.

Equation (1) is also referred to as the advanced diffusion model (ADM) and describes the transport of solutes within the membrane. In RO membrane systems, solute transport is governed by both diffusion (random molecular motion) and advection (bulk fluid movement). To solve ADM, we apply FEM using the stiffness approach. We still employ the SDM to capture the osmotic pressure-driven movement of the solvent through the RO membrane. By coupling the ADM and SDM, we are able to investigate the three fundamental forces: advection, osmotic pressure, and diffusion, as they jointly govern solute and solvent transport in RO membranes.

By converting the ADM model to a system of algebraic equations and solving for the unknown parameters, such as $c_i$ and $J_i$ , the FEM provides a reliable numerical approach for obtaining the solution of the new extend model (ADM and SDM) equation. The input parameters for the new extend models are: membrane geometry, feed water, and operating conditions.

The initial step to solve the ADM and SDM by FEM is to define the geometry of the RO membrane, the boundary conditions, and the initial ion concentrations, i.e., feed concentration. The general solute transport equation in one dimension (i.e., x-direction) can be described as follows:

$\frac{\partial C}{\partial t}+v\frac{\partial C}{\partial x}=\frac{\partial }{\partial x }\left(D_{\text{eff}}\frac{\partial C}{\partial x}\right)+R$ (2)

where $C\left(x,t\right)$ is the solute concentration with unit mol/m³, $t$ is the time of solute transport in RO membrane (Second) and $v$ is the advective velocity (m/s), $D_{\text{eff}}$ is the effective solute diffusivity (m²/s), and $R$ is the rejection or adsorption factor.

In Equation (2), the left-hand side represents the rate of change of solute concentration and the advective transport of solute. The right-hand side represents diffusion and any reaction terms. $D_{\text{eff}}$ is the effective solute diffusivity and mathematically can be computed by:

$D_{\text{eff}}=D_0{\text{e}}^{-\alpha C}$ (3)

where $D_0$ is the bulk diffusivity (m²/s), and $\alpha$ is the solute interaction coefficient. Several assumptions are made to solve the Advection-Diffusion Model (ADM), including: 1) steady-state, where $\frac{\partial C}{\partial t}=0$ , and 2) small concentration changes in which ${\text{e}}^{-\alpha C}\approx 1-\alpha C$ , and no adsorption or degradation and mathematically can be described as follows:

$v\frac{\partial C}{\partial x}=\frac{\partial }{\partial x }\left(D_0{\text{e}}^{-\alpha C}\frac{\partial C}{\partial x}\right)$ (4)

By expanding the right-hand side and rearranging Equation (4), it can be written as follows:

$v\frac{\partial C }{\partial x }+\alpha D_0{\text{e}}^{-\alpha C}{\left(\frac{\partial C}{\partial x}\right)}^2=D_0{\text{e}}^{-\alpha C}\frac{{\partial }^{2}C}{\partial x^2}=D_0{\text{e}}^{-\alpha C}\frac{{\partial }^{2}C}{\partial x^2}$ (5)

Since $C$ is small, the term $\alpha C$ is also small, we neglect higher-order terms ${\alpha }^{2}C$ , reducing Equation (5) to:

$v\frac{\partial C}{\partial x }+\alpha D_0{\left(\frac{\partial C}{\partial x}\right)}^2=D_0 \frac{{\partial }^{2}C}{\partial x^2}-\alpha D_{0}C\frac{{\partial }^{2}C}{\partial x^2}$ (6)

Also, since $\alpha C$ and $\alpha D_0$ are small value compared to the leading transport terms, non linear term in Equation (6) such as $\alpha D_0{\left(\partial C/\partial x\right)}^2$ and $\alpha D_{0}C\left(\frac{{\partial }^{2}C}{\partial x^2}\right)$ can be neglected without effecting on the solution accuracy. By linearizing ADM Equation (6) can be simplified to:

$v\frac{\partial C}{\partial x}=D_0\frac{{\text{d}}^{2}C}{\text{d}{x}^2}$ (7)

This linearization is based on the assumption of low solute concentrations and a weak coupling effect, meaning the nonlinear interaction terms are too small to effect on the transport behavior. Under these conditions, Equation (7) suggests a valid approximation. In our analysis, the parameter α < 1, supports the weak coupling assumption. However, if α ≥ 1, the coupling becomes strong, and the nonlinear terms can no longer be neglected.

2.2. Finite Element Method to Solve ADM

In this section, the ADM is first solved (using FEM with the stiffness approach) to obtain the solute concentration profiles and flux distributions along the RO membrane length. The results from this analysis describe how solutes move through the feed channel and are ultimately rejected at the membrane surface. This concentration profile acquired from ADM is used to calculate the osmotic pressure term in SDM.

Finite Element Method (FEM) is the stiffness approach, and the stiffness matrix relates the nodal values of the unknown function to applied forces. The following is the general form:

$K\left\{C\right\}=\left\{F\right\}$ (8)

where K is the stiffness matrix that accounts for diffusion and advection effects, $\left\{C\right\}$ is the vector of unknown nodal values, and $\left\{F\right\}$ is the force vector that incorporates boundary conditions and source terms. For an N-node 1-D finite element mesh, the global stiffness matrix is assembled from local element matrices:

$K=\left[\begin{array}{ccccc}{K}_{11}& K_{12}& 0& \cdots & 0\\ k_{21}& K_{22}& K_{23}& \cdots & 0\\ 0& K_{32}& K_{33}& \cdots & 0\\ ⋮& ⋮& ⋮& \cdots & K_{N-1,N}\\ 0& 0& 0& K_{N,N-1}& K_{NN}\end{array}\right]$ (9)

where matrix K represents the combined effects of diffusion and advection. It has contribution from both the diffusion term $D_0\frac{{\text{d}}^{2}C}{\text{d}{x}^2}$ and the advection term $v\frac{\text{d}C}{\text{d}x}$ .

To solve Equation (7) by FEM using the linear stiffness approach, we solve it over a domain $0<x<L_0$ and apply both the Neumann and Dirichlet boundary conditions, which mathematically can be represented as follows:

$\begin{array}{l}\text{Dirichlet condition at}x=0:C\left(0\right)=C_o\\ \text{Neumann condition captures the zero flux at}x=L:{\frac{\text{d}C}{\text{d}x}|}_{x=L}\approx 0\end{array}$ (10)

To minimize the residual error across the domain, the strong form is converted into a weak form [17]. That is, both sides of Equation (7) is multiplied by a test function ${\varphi }_i$ . We assume two nodes, and integrate over the domain $0<x<L_0$ , the following is obtained:

${\displaystyle {\int }_0^{L}v\frac{\text{d}C}{\text{d}\chi }{\varphi }_i \text{d}x}={\displaystyle {\int }_0^L{D}_0\frac{{\text{d}}^{2}C}{\text{d}{\chi }^2}{\varphi }_i\text{d}x}$ (11)

Solving the diffusion term in Equation (11), by integration by parts, and applying Neumann condition in Equation (10), the weak form can be written as follow:

${\displaystyle {\int }_0^{L}v\frac{\text{d}C}{\text{d}\chi }{\varphi }_i\text{d}x}={\displaystyle {\int }_0^L{D}_0\frac{\text{d}C}{\text{d}x}\frac{\text{d}{\varphi }_i}{\text{d}x}\text{d}x}$ (12)

Also $C\left(x\right)$ is approximated using shape functions as follow:

$C\left(x\right)=C_1{\varphi }_1\left(x\right)+C_2{\varphi }_2\left(x\right)$ (13)

The derivative of Equation (13) is given by,

$\frac{\text{d}C}{\text{d}x}=C_1\frac{\text{d}{\varphi }_1}{\text{d}x}+ C_2\frac{\text{d}{\varphi }_2}{\text{d}x}$ (14)

For a linear shape function:

$\frac{\text{d}C}{\text{d}x}=C_1\left(-\frac{1}{\Delta x}\right)+C_2\frac{1}{\Delta x}$ (15)

where $\Delta x=x_2-x_1$ .

Now substitute $C\left(x\right)$ and $\frac{\text{d}C}{\text{d}x}$ into Equation (12), it can be rewritten as follows:

${\displaystyle {\int }_0^{L}v\left(C_1\frac{\text{d}{\varphi }_1}{\text{d}x}+C_2\frac{\text{d}{\varphi }_2}{dx}\right){\varphi }_i\text{d}x}={\displaystyle {\int }_0^L{D}_0\left(C_1\frac{d{\varphi }_1}{\text{d}x}+C_2\frac{\text{d}{\varphi }_2}{\text{d}x}\right)\frac{\text{d}{\varphi }_{i }}{\text{d}x}\text{d}x}$ (16)

This leads to the final element stiffness matrix to the following:

$K_{ij}={\displaystyle {\int }_{x_i}^{x_j}\left(D_0\frac{\text{d}{\varphi }_i}{\text{d}x}\frac{\text{d}{\varphi }_j}{\text{d}x}+v{\varphi }_i\frac{\text{d}{\varphi }_j}{\text{d}x}\right)\text{d}x}$ (17)

The linear shape functions for a two-node element are [18]:

${\varphi }_1\left(x\right)=\frac{x_2-x}{\Delta x},{\varphi }_2\left(x\right)=\frac{x-x_1}{\Delta x}$ (18)

where $\Delta x=x_2-x_1$ is the element length. By taking the derivative of Equation (18), it can be represented as follows:

${\varphi }_1\left(x\right)=-\frac{1}{\Delta x},{\varphi }_2\left(x\right)=\frac{1}{\Delta x}$ (19)

Substituting the shape function derivative (Equation (19)) into the diffusion term in Equation (16), it will give the following:

${\displaystyle {\int }_{x_i}^{x_j}{D}_0\left(\frac{1}{\Delta x}\frac{1}{\Delta x}\text{d}x\right)}=D_0\frac{1}{\Delta x^2}{\displaystyle {\int }_{x_i}^{x_j}\text{d}x}$ (20)

Since the integral evaluates to $\Delta x$ :

$D_0\frac{1}{\Delta x^2}\times \Delta x=\frac{D_0 }{\Delta x}$ (21)

Thus, for elements stiffness, we have

$K_{11}^{\text{diffusion}}=K_{22}^{\text{diffusion}}=\frac{D_0 }{\Delta x},K_{12}^{\text{diffusion}}=K_{21}^{\text{diffusion}}=-\frac{D_0 }{\Delta x}$

Since the advection terms contain one linear function and one constant derivative, the integral of a linear function over an element gives a factor of $\frac{1}{2}$ , thus the advection form in Equation (17) can be evaluated as follows:

${\displaystyle {\int }_{x_1}^{x_2}v{\varphi }_1\frac{\text{d}{\varphi }_1}{\text{d}x}\text{d}x}\approx v\times \frac{1}{2}$ (22)

This result in:

$K_{11}^{\text{advection}}=\frac{v}{2}$ ,

$K_{22}^{\text{advection}}=-\frac{v}{2}$ ,

$K_{12}^{\text{advection}}=-\frac{v}{2}$ ,

$K_{21}^{\text{advection}}=\frac{v}{2}$

Adding the diffusion and advection contributions, the stiffness of two nodes can be represented as follows:

$K_e=\left[\begin{array}{cc}\frac{D_0}{\Delta x}+\frac{v}{2}& -\frac{D_0}{\Delta x}-\frac{v}{2}\\ -\frac{D_0}{\Delta x}+\frac{v}{2}& \frac{D_0}{\Delta x}-\frac{v}{2}\end{array}\right]$ (23)

The discretization is applied along the membrane’s length, which is 0.356 meters. The domain is divided into N = 50, 100, 500, and 1000 mesh elements, respectively. As shown in Figure 1, the membrane is treated as a sheet and meshed into uniform triangular elements. Since transport occurs within the membrane, it is modeled as a 2-D field divided into equal triangular grids. FEM is then used to solve for the permeate concentration based on Equation (7). This is followed by calculating the solvent flux using SDM, and finally, the transmembrane pressure is determined from the SDM’s solvent equation.

Figure 1. Membrane flat sheet divided up to uniform triangle mechs.

The stiffness matrix represents the stiffness of an individual finite element in the mesh [19]. In our analysis it is a correlation between the forces imposed on the solute transport and displacement at the nodes in the systems. To initialize the stiffness matrix in MATLAB, we create a global matrix of size (N + 1) × (N + 1). In FEM, Equation (7) is approximated by a system of algebraic equation using polynomial approximation over a number of elements.

To solve those algebraic equations, we need to assemble the stiffness matrix in Equation (8). The process iterates through each element in the mesh, starting from $i=1$ to the total number of the elements, as seen in Figure 2. The element length, $\text{d}x$ , is calculated as the distance between nodes $i$ and $i+1$ , and it defines the spatial step size used in the FEM. A diffusivity value is then selected from the predefined matrix $D$ . Next, the local stiffness matrix ($K_e$ ) is assembled using Equation (23) and added to the corresponding section of the global stiffness matrix $K$ . After that, the local mass matrix is assembled, which accounts for the time accumulation effects in each element. This local mass matrix is then added to the global mass matrix M [20].

The loop index $i$ is incremented, and the process repeats for the next element. Once all elements have been processed, the global matrices K and M are fully assembled and ready to be used in solving the FEM system of equations. We loop over the nodes, compute the integrals over the nodes and then solve for the linear system to calculate for the overall permeate solute flux, as illustrated in Figure 2. In this analysis, we find the values of solute concentrations at the nodes of the finite element mesh and then interpolate the concentration profile across the domain of the ADM, as shown in Figure 3. Finally, To verify the accuracy and reliability of the solution, the ADM and SDM data are compared with the experimental findings.

Figure 2. Flowchart for solving the stiffness matrices in ADM (N = 100).

Figure 3. Flowchart for solving the ADM using FEM with the stiffness approach (N = 100).

By performing the analysis mentioned previously, ADM captures the solute-side transport dynamics, while SDM captures the solvent-side flow driven by osmotic gradients.

2.3. COMSOL Multiphysics for Solving ADM

In Section 2.1, we have explained how ADM can be solved using FEM through the stiffness matrix approach. In this section, we take a different approach by using Computational Fluid Dynamics (CFD) in COMSOL Multiphysics to analyze and obtain the concentration profile.

To begin with, the model geometry is divided into three separate rectangular sections, each representing a specific region: the feed, the active membrane (transport zone), and the permeate side. These sections are modeled as solid rectangular sheets. The dimensions for each region are illustrated in Figure 4 and the values are given in Table 2.

Figure 4. Membrane geometry in COMSOL multiphysics.

Table 2. Geometric dimensions and offsets of RO membrane regions.

Region	Width (m)	Height (m)	Offset (m)
Feed	0.0508	0.254	0
Membrane (active transport area)	0.0254	0.254	0.0508
Permeate	0.0508	0.254	0.0762

The material properties are then defined using seawater for the feed region, which has a given dynamic density and viscosity. Then, under Material Transport, two physics modules are added:

1) Diluted species transport, accounting for convective transport mechanisms.

2) Laminar Flow, considering steady-state fluid behavior.

In the diluted species transport physics, the velocity field is specified to incorporate convective transport of solutes across the RO membrane. The diffusion coefficient for the membrane region is set to 1 × 10⁻⁶ m²/s, and for the feed region, it is larger, at 1 × 10⁻⁵ m²/s.

The initial concentration at the feed inlet boundary ($x=0$ ) is set to match the conditions previously used in the FEM model. At the membrane outlet boundary ($x=L$ ) the solute flux is assumed to be zero, as described by Equation (10). A partition condition of $K_c=1$ is applied, indicating that the solute concentration remains continuous across the membrane interface, with no additional partitioning effects. In the laminar flow physics, an inlet velocity of 0.01 m/s is set. Also, the outlet flow condition incorporating osmotic pressure effects.

A finer mesh is created to match the resolution used in the earlier MATLAB simulations, as shown in Figure 5. The concentration profile is then extracted along the active boundary to evaluate the interface behavior [21].

Figure 5. Finite element mesh of the RO membrane geometry in COMSOL.

3. Results from the FEM Stiffness Method

3.1. Concentration and Flux Profiles along the Membrane Length

Figure 6 displays the concentration and flux behavior along the membrane length during solute transport at n = 100 based on the ADM in Equation (7). This analysis is important because separating between diffusive and advective fluxes along $x$ which reveals the dominant transport mechanism that drives solute and solvent inside the RO membrane.

In Figure 6, the blue dashed line represents the solute concentration $C\left(x\right)$ , which starts high near the feed side (at $x=0$ ) and rapidly declines along the membrane, indicating strong solute rejection and minimal penetration through the membrane. The diffusive flux (brown dashed line) and advective flux (yellow dashed line) reaches the highest near the inlet and diminish quickly, showing that most of the transport activity occurs within a short region near the membrane surface. Along the membrane, the green dashed line, which represents the total flux or the sum of diffusion and advection, likewise drops off abruptly. The steep concentration gradient close to the inlet and the high advective velocity, which cause the solute to move quickly at the entrance before rapidly declining downstream, are the main causes of this decline.

Figure 6. Concentration and flux profiles along the membrane length using the FEM stiffness approach (Feed Salinity = 32 k ppm).

Overall, the results show that solute movement is highly localized near the membrane’s inlet, and the fluxes approach zero beyond a certain distance, indicating a dominant boundary layer effect in the transport process.

The ADM result in Figure 6 shows a very sharp concentration drop near the inlet along the membrane length, followed by near-zero concentration for most of the domain. Both advection and diffusion forces contribute to this sharp drop; the inlet flow quickly pushes solute downstream, resulting in the concentration decreasing more abruptly.

Overall, the findings indicate a dominant boundary layer effect in the transport process, with solute movement being highly localized close to the membrane’s inlet and fluxes approaching zero beyond a certain distance.

The ADM result in Figure 6 shows a very sharp concentration drop near the inlet along the membrane length, followed by near-zero concentration for most of the domain. This steep decline is influenced by both advection and diffusion forces, where the inlet flow rapidly pushes solute downstream, causing the concentration to decrease much more abruptly.

3.2. Effect of Changing Boundary Condition Configurations in ADM

Since boundary conditions can highly influence the behavior and the accuracy of numerical solutions, this section investigates the concentration and solute flux along the membrane length under varying inlet and outlet boundary conditions. Different scenarios of experiments are performed to help us better understand how concentration gradients and transport dynamics respond to changes in physical constraints.

Figure 7 shows the concentration profiles $C\left(x\right)$ and solute fluxes $q\left(x\right)$ across the membrane length under different Dirichlet and Neumann conditions as described in Equation (10). By exploring multiple outlet flux settings, it guarantees that the model will always be adaptable and representative of real-world membrane operation circumstances. This figure highlights the boundary-layer location and mass-transfer direction that determine the outlet concentration.

In all four subplots, the blue solid lines represent the solute concentration, while the orange dashed lines show the corresponding flux. Despite the changes in the outlet flux condition $q_L$ ($J_{\text{solute}}$ ), the unit of $q_L$ is (mol/m³/m), the concentration profile consistently and quickly drops at the membrane inlet ($x=0$ ) and flattens out. However, the flux behavior is highly sensitive to the boundary condition. For example, $q_L=0$ , the flux approaches zero gradually; when $q_L$ is positive or negative, the flux either maintains a steady value or shifts in response to the imposed boundary influence.

A negative outlet flux means the solutes are being pushed into the membrane from the right side instead of exiting, which creates a reversed concentration gradient. This backward flow behavior is reflected in the shifted flux profile, even though the concentration still drops from the inlet. This creates an unusual situation where diffusion and possibly advection are influenced in the opposite direction.

A positive $q_L$ , $C\left(x\right)$ drops sharply near the inlet. This is indicating a thin depletion (boundary) layer where most transport occurs, while $q\left(x\right)$ stays nearly constant along $x$ .

Figure 7. Effect of changing boundary conditions on concentration and flux profiles along the membrane length.

3.3. Mesh Sensitivity Analysis

Figure 8 illustrates how the concentration profile evolves along the membrane length when different mesh densities (50, 100, 500, and 1000 elements) are used in the finite element model. Figure 8 and Figure 9 validate grid-independence to guarantee that outlet concentration and flux are accurate and not artifacts of an inadequately refined mesh, as the abrupt boundary layer requires appropriate mesh resolution.

Effective solute rejection is demonstrated in every case by the concentration $C\left(x\right)$ rapidly approaching zero after dropping sharply at the membrane inlet. As the mesh becomes finer, the transition in the concentration profile appears smoother and more precisely defined, particularly near the inlet where the gradient is steepest. This illustrates the importance of mesh refinement in enhancing numerical stability and capturing abrupt concentration changes, particularly in regions with rapid variation.

Figure 9 shows that in every case, the solute flux is highest at the beginning of the membrane and quickly drops and approaches zero. This sharp decline highlights a strong concentration gradient at the inlet, with much less transport happening farther along the membrane. As the number of elements increases from 50 to 1000, the plots become more detailed and consistent, especially around the inlet where the flux changes rapidly. This suggests that using a finer mesh provides a more accurate representation of the diffusion process [22]. In this case, a mesh size of 500 to 1000 elements is sufficient to capture the true behavior.

Figure 8. Effect of mesh refinement on solute concentration profile along the membrane length.

Figure 9. Effect of mesh refinement on flux profile along the membrane length.

In a discretized domain, the concentration gradient $\frac{\text{d}c}{\text{d}x}$ is approximated by finite differences between adjacent nodes. When the number of mesh 50 elements, the spacing between nodes is relatively wide. This makes it challenging to depict abrupt concentration changes accurately, particularly near the membrane inlet where gradients are steep. As a result, the computed flux at the inlet tends to be lower than it should be because the steep drop in concentration is numerically smoothed out. In contrast, using a finer mesh with 500 or 1000 elements provides much better resolution. Sharper gradients can be more accurately resolved by the model because nodes are spaced closer together. This improves the accuracy of the flux calculation near the inlet and results in a higher estimate of solute transport at the feed side, better capturing the physics of diffusion in that region. However, employing a finer mesh also increases the computational time required to generate results [23].

By changing the mesh size, Figure 10 shows how good our numerical approach is in predicting concentration and flux at the membrane outlet. The Concentration Error Convergence plot shows the error in our calculated concentration. The resulting error is extremely small, around 10⁻¹². As the grid size (h) increases, the error decreases slightly. This indicates our numerical method is quite stable and produces accurate results, even when fewer computational points are used. The Outlet Flux Error plot shows the error in the calculated outlet solute flux. Here, the error can be effectively pushed to zero for all tested mesh sizes. That means the model perfectly matches our boundary conditions at the outlet. Overall, these plots confirm that our numerical method is reliable and robust.

Figure 10. Mesh convergence of the solute concentration and outlet flux.

4. Results from COMSOL Multiphysics Simulation

Figure 11 shows a 3-D concentration distribution across RO membrane module at 32 k ppm feed salinity. It shows the solute concentration distribution within a cylindrical geometry, likely simulating RO membrane module. The simulation uses mol/m³ as the unit for concentration. Areas with high solute buildup are visually indicated by the red zone in the center. (above 500 mol/m³) while the blue regions near the outer boundaries represent zones of low concentration, approaching that of pure solvent. This gradient suggests that solutes accumulate near the feed side and are progressively rejected or diluted when moving outward, which is a typical behavior in RO membranes.

Figure 11. Concentration distribution across RO membrane at 32 k ppm feed salinity.

The model is sliced open to reveal the internal concentration structure across the membrane’s length, providing a clear view of how solute transport behaves inside the domain. This helps evaluate membrane performance and supports the validation of numerical models such as FEM. The sharp gradient seen inside the module is further supported by the color bar on the right, which quantifies the concentration field and displays a range of 0 to 550 mol/m³.

Figure 12 presents the concentration distribution within a rectangular section of a membrane system, likely from a simulation of diffusion process. The color map displays solute concentration in mol/m³, with red areas indicating high concentration (above 500 mol/m³) and blue areas indicating low concentration (approaching zero). The active membrane zone where mass transfer is most intense is indicated by the sharp gradient between red and blue, particularly in the narrow vertical band. Streamlines and arrows overlaid on the plot illustrate the direction and pattern of solute flux.

Figure 12. Streamline flow distribution across RO membrane in COMSOL.

On the left side, the vertical arrows show a fairly uniform upward movement of concentrated solute, indicating flow parallel to the membrane surface probably representing feed flow. As the solute approaches the membrane surface (the red to blue transition), the direction shifts, and some of the solute begins to diffuse or be drawn perpendicularly across the membrane. The streamlines curve and disperse in the blue area on the right, illustrating how the concentration decreases on the permeate side. The contours in this region represent equal concentration levels and help illustrates the diffusion path of the solute. Overall, Figure 12 makes it abundantly evident how the solute is moved from the high concentration feed region to the low concentration permeate side across the membrane, with the narrow transition zone exhibiting the strongest mass transfer.

The results of the COMSOL simulation and the FEM-stiffness method agree well, especially when it comes to the evolution of concentration and flux across the membrane. In the COMSOL analysis, the steep concentration transition at the membrane interface is clearly visible, with high solute levels gradually dropping from the feed side (red region) to the permeate side (blue side). This matches the pattern observed in the FEM output, especially when using finer mesh resolutions. Also, the streamline flows in the COMSOL model, which illustrates the direction and strength of the flux. This shows the same trend as the FEM results. That is, there exists a strong flux at the inlet, and it weakens along the membrane length. This consistency between the two approaches reinforces the accuracy of the FEM-stiffness method in modeling solute transport in RO systems.

Table 3 compares COMSOL Multiphysics and a stiffness-based FEM for three randomly selected boundary-layer thicknesses. For each case, the outlet concentration was calculated. The methods show excellent agreement, with differences of only ~0.4% - 0.7%

Table 3. Comparison between COMSOL multiphysics and the FEM-stiffness approach.

Boundary layer thickness (m)	COMSOL multiphysics (ppm (×10^3))	FEM-stiffness approach (ppm (×10^3))	Error (%)
0.0245	31.800	31.960	0.503%
0.0310	30.610	30.815	0.670%
0.0445	30.410	30.536	0.414%

5. Solvent Analysis of SDM

In this section, we calculate the design parameters of the RO membrane under varying feed pressures and water permeability values. Water permeability, which indicates the rate of water flux through the membrane into the permeate, typically ranges from 2.78 × 10⁻¹² to 5.56 × 10⁻¹² m/(s·Pa) (or 1 - 2 LMH/bar) as noted by [24]. Table 3 presents the net driving pressure, solvent flow rate, and salt rejection at different salinity levels and membrane permeability values. Table 4 also explains the influence of feed salinity and water permeability on reverse osmosis (RO) membrane performance in terms of net driving pressure, solvent (water) flow rate, and salt rejection efficiency. Four salinity levels from 32,000 to 35,000 ppm were evaluated across five different water permeability values ranging from 2.78 × 10⁻¹² to 5.56 × 10⁻¹² m/(s·Pa).

As feedwater salinity increases, the net driving pressure decreases slightly from 2.79 MPa at 32,000 ppm to 2.53 MPa at 35,000 ppm. This trend reflects the increasing osmotic pressure that offsets the applied hydraulic pressure. As a result, even when permeability stays constant, water flux (Q) through the membrane likewise decreases with increasing salinity. For example, at a permeability of 2.78 × 10⁻¹² m/(s∙Pa), the flow rate decreases from 0.405 m³/day at 32,000 ppm to 0.368 m³/day at 35,000 ppm. Furthermore, at each salinity level, increasing membrane water permeability results in a proportional increase in flow rate. This confirms that membrane permeability is a key design parameter in achieving higher water productivity. For instance, at 33,000 ppm, increasing permeability from 2.78 × 10⁻¹² to 5.56 × 10⁻¹² m/(s∙Pa) increases the flow rate from 0.392 to 0.785 m³/day, effectively doubling the water output.

Table 4. Effect of salinity and water permeability on RO membrane performance.

Salinity (ppm)	Net driving pressure ΔPnet × 10−6 (Pa)	Water permeability Aw × 10−12 (m/(s·Pa))	Flow rate (Q) (m3/day)	Salt rejection (%)
32,000	2.79	2.78	0.405	100
		3.47	0.506	100
		4.17	0.607	100
		4.86	0.708	100
		5.56	0.809	100
33,000	2.71	2.78	0.392	100
		3.47	0.49	100
		4.17	0.589	100
		4.86	0.687	100
		5.56	0.785	100
34,000	2.62	2.78	0.38	100
		3.47	0.475	100
		4.17	0.57	100
		4.86	0.665	100
		5.5	0.76	100
35,000	2.53	2.78	0.368	100
		3.47	0.46	100
		4.17	0.552	100
		4.86	0.644	100
		5.56	0.735	100

In all scenarios, the salt rejection remains constant at 100%, indicating that changes in water permeability and salinity within this tested range do not compromise membrane selectivity. While higher permeability improves water flux, elevated salinity imposes limits due to osmotic pressure buildup, illustrated in Figure 13 and Figure 14.

Figure 13. Effect of water permeability and salinity on permeate flow rate.

Figure 14. Relationship between net driving pressure, salinity, and permeate flow rate.

Figure 15 explains the effect of net driving pressure and membrane permeability on permeate flow rate. Also, it Validates the linear dependency of flow rate on both membrane permeability and net driving pressure. Obviously selecting membranes with higher values of $A_w$ and operating at greater net pressures can significantly improve water production in RO systems. However, such improvements must be balanced against energy costs and system limits. Also, the x-intercepts in left plot (~0 flow rate) are around 0 Pa, consistent with the physical meaning that there is no flow without pressure difference. The right plot 3D view reinforces the same trend seen in the 2D plot but adds depth by showing how flow rate depends on the combined effect of pressure and permeability. The surface smoothly rises, confirming that both increasing permeability and pressure increase the water production.

Figure 15. Effect of net driving pressure and membrane permeability on permeate flow rate.

6. Comparing with Commercial RO Modules

The SW30-2514 and SW30-4021 are the two commercial RM modules manufactured by DuPont (previously called Dow FilmTec). Both are conventional spiral-wound RO membranes often found in small- to medium-scale seawater desalination facilities. The SW30-4021 is larger than SW30-2514 and has area of around 7.9 m², while the SW30-2514 has a membrane area of roughly 2.6 m². Both modules run best at the feed pressure of 5500 kPa. Table 5 and Table 6 present a comparison between the outcomes of the extended ADM-SDM model, solved using the finite element stiffness method, and the manufacturer-reported data at a feed pressure of 5500 kPa.

We validate the extended model (ADM-SDM) by comparing its results with experimental data from a commercial seawater reverse osmosis (SWRO) module. Permeate flux and salt rejection are the two main performance metrics used to evaluate the accuracy of membrane models. The results demonstrate strong agreement between the simulated and manufacturer’s reported data, confirming the accuracy of the FEM model. By including advection effects, the extended model improves flow rate predictions over the traditional SDM and lowers calculation errors by more accurately representing the physical transport processes inside the membrane.

Table 5. The comparison between ADM-SDM and SW30-2514.

Design parameter	SW30-2514	ADM-SDM	Error (%)
Permeate flux (m3/day)	0.60	0.607	0.7
Salt rejection (%)	99.40	100	0.60

Table 6. The comparison between SDM and SW30-4021.

Design parameter	SW30-4021	ADM-SDM	Error (%)
Permeate flux (m3/day)	7.4	7.204	2.6
Salt rejection (%)	99.40	100	0.60

Due to idealized assumptions in diffusivity and membrane homogeneity, the model’s constant salt rejection prediction of 100% may somewhat exceed real-world performance. This validation confirms that by using FEM, the expanded ADM-SDM framework offers a strong tool for simulating solute and solvent transport in RO membranes. Its precision in anticipating commercial module performance supports its application in membrane selection, and performance analysis.

To validate our model against experimental results, we compared our salt rejection results with findings from two independent studies. For instance, [25] reported a salt rejection of 97.5% at a feed pressure of 1100 kPa and a feed concentration of 2 kg/m³. Our model predicted 96.8% under the same conditions, which is only 0.72% off. This indicates strong agreement with the report.

Also, [26] reported a salt rejection of approximately 99% at 40 bar (4000 kPa) and 30,000 ppm. Our model predicted 98.78% under the same conditions, which is only 0.22% off. This indicates strong agreement with the report.

Like other models, the ADM-SDM framework has limitations in fully capturing transport mechanisms that occur within the thin boundary layers near the membrane surface because the flow behavior becomes complex and less predictable. In addition, mathematical models like the ADM-SDM framework does not fully reflect events such as membrane fouling, in which solute particles build up and clog membrane pores over time, hence gradually reducing performance.

7. Conclusions

This study proposes an extended mass transport model for reverse osmosis (RO) membranes by coupling the ADM with SDM and solving it using FEM with a linear stiffness approach. The results demonstrate good agreement between ADM-SDM via COMSOL and experimental data from commercial RO modules, validating the model’s accuracy. The FEM stiffness method showed that solute concentration and flux are strongly localized near the membrane inlet. Both advective and diffusive fluxes peak at the feeding side and decrease rapidly downstream, confirming the dominant role of boundary layer effects in solute rejection.

Varying boundary conditions highly affected solute flux behavior while concentration profiles remained steep at the inlet and flattened along the membrane. Reversed gradients were produced by negative outlet flux conditions, demonstrating how sensitive flux predictions are to boundary conditions. Changing the mesh size improved the accuracy of concentration and flux profiles near the inlet, where gradients are steepest. A mesh resolution of 500 - 1000 elements was sufficient to capture sharp variations, while error convergence plots confirmed that the numerical method is both stable and reliable.

The extended Advection-Diffusion solution (ADM-SDM) model showed close agreement with manufacturer-reported data for commercial modules (SW30-2514 and SW30-4021). The model’s predictions were accurate because the errors in permeate flux forecasts stayed below 3% and the estimates for salt rejection were slightly higher than the reported values.

The salt rejection is calculated using following formula:

$\text{Salt Rejection}\left(\%\right)=\left(1-\frac{C_p}{C_f}\right)\times 100\%$

where $C_p$ and $C_f$ are the solute concentrations in the permeate and feed, respectively.

Because $C_p\ll C_f$ the calculated salt rejection was 99.99%, which was rounded to 100%. While the extended ADM-SDM framework provides reliable predictions of solute transport and permeate flux, certain model assumptions contribute to differences between the estimated salt rejection and the manufacturer specifications of commercial RO membranes. These assumptions include the linearization based on low solute concentrations and weak coupling effects, as well as the exclusion of electromigration.

Overall, including the advection effects improved the flux predictions compared to traditional SDM approaches. It is precisely calculating salt rejection and permeate flux. The model offers a scalable and computationally effective tool for designing, optimizing, and assessing the performance of reverse osmosis (RO) systems in practical desalination applications. In addition to accurately forecasting performance metrics, the extended ADM-SDM framework offers recommendations for creating anti-fouling and energy-saving techniques for RO desalination systems.

The accumulation of solute at the membrane surface relative to the bulk solution is known as concentration polarization [27]. There are significant design and operational implications to our discovery that flux is highly localized at the membrane inlet. High concentrations in this area raise the local osmotic pressure, which decreases flux efficiency and the effective driving force for water transport. By determining the concentration behavior at the inlet, as we performed in our analysis, strategies like staged feed distribution or inlet pressure adjustment can be used to more evenly redistribute flux along the membrane. As a result, fouling is minimized, which lowers water flow, boosts pressure drop, shortens membrane life, and raises operating expenses [28]. Although the ADM-SDM model is quite effective and accurate, it is still unable to account for membrane fouling and electromigration force. Therefore, the integration of fouling behavior and transient transport processes should be the main focus of future research. Predictive simulations for RO systems would become much more realistic as a result.

Acknowledgements

The leading author would like to thanks Dr. Tongdan Jin, Professor of Ingram School of Engineering at Texas State University, for reviewing and revising the manuscript as well as providing valuable comments during the writing and submission of this manuscript.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Ang, W.L. and Mohammad, A.W. (2015) Mathematical Modeling of Membrane Operations for Water Treatment. In: Basile, A., Cassano, A. and Rastogi, N.K., Eds., Advances in Membrane Technologies for Water Treatment, Elsevier, 379-407.[CrossRef]
[2]	Giorno, L., Strathmann, H. and Drioli, E. (2015) Mathematical Description of Mass Transport in Membranes. In: Drioli, E. and Giorno, L., Eds., Encyclopedia of Membranes, Springer, 1-4.[CrossRef]
[3]	Liu, L. (2018) Application of a Hydrodynamic and Water Quality Model for Inland Surface Water Systems. In: Malcangio, D., Ed., Applications in Water Systems Management and Modeling, InTech.[CrossRef]
[4]	Ault, J.T. and Shin, S. (2025) Physicochemical Hydrodynamics of Particle Diffusiophoresis Driven by Chemical Gradients. Annual Review of Fluid Mechanics, 57, 227-255.[CrossRef]
[5]	Ayrapetyan, S. (2020) The Quantum-Mechanical Sensitive Na/k Pump Is a Key Mechanism for the Metabolic Control of Neuronal Membrane Function. Open Journal of Biophysics, 10, 59-83.[CrossRef]
[6]	Tiraferri, A., Malaguti, M., Mohamed, M., Giagnorio, M. and Aschmoneit, F.J. (2023) Standardizing Practices and Flux Predictions in Membrane Science via Simplified Equations and Membrane Characterization. npj Clean Water, 6, Article No. 58.[CrossRef]
[7]	Dai, R., Yang, Z., Qiu, Z., Long, L., Tang, C.Y. and Wang, Z. (2022) Distinct Impact of Substrate Hydrophilicity on Performance and Structure of TFC NF and RO Polyamide Membranes. Journal of Membrane Science, 662, Article ID: 120966.[CrossRef]
[8]	Giorno, L. and Drioli, E. (2000) Biocatalytic Membrane Reactors: Applications and Perspectives. Trends in Biotechnology, 18, 339-349.[CrossRef] [PubMed]
[9]	Marx, D., Joss, L., Hefti, M. and Mazzotti, M. (2016) Temperature Swing Adsorption for Postcombustion CO2 Capture: Single-and Multicolumn Experiments and Simulations. Industrial & Engineering Chemistry Research, 55, 1401-1412.[CrossRef]
[10]	Tsirlin, A.M. (2023) Thermodynamic Analysis of Irreversible Processes of the Reverse Osmosis and Diffusion Machines. Journal of Engineering Physics and Thermophysics, 96, 555-567.[CrossRef]
[11]	Du, Y., Wang, Z., Cooper, N.J., Gilron, J. and Elimelech, M. (2022) Module-scale Analysis of Low-Salt-Rejection Reverse Osmosis: Design Guidelines and System Performance. Water Research, 209, Article ID: 117936.[CrossRef] [PubMed]
[12]	Matsuura, T., Guiver, M.D., Lau, W.J. and Ismail, A.F. (2025) Transport in Reverse Osmosis Membranes: Observations and Comments on the Pore Flow Model versus the Solution-Diffusion Model. Frontiers in Membrane Science and Technology, 4, Article 1513591.[CrossRef]
[13]	Wang, L., Cao, T., Dykstra, J.E., Porada, S., Biesheuvel, P.M. and Elimelech, M. (2021) Salt and Water Transport in Reverse Osmosis Membranes: Beyond the Solution-Diffusion Model. Environmental Science & Technology, 55, 16665-16675.[CrossRef] [PubMed]
[14]	Biesheuvel, P.M., Porada, S., Elimelech, M. and Dykstra, J.E. (2022) Tutorial Review of Reverse Osmosis and Electrodialysis. Journal of Membrane Science, 647, Article ID: 120221.[CrossRef]
[15]	Taylor, R.F., Martinez-Jimenez, F. and Logan, B.E. (2025) Tuning Initial pH to Decrease Salt Ion Transport in Saltwater Electrolysis. Electrochemistry Communications, 171, Article ID: 107858.[CrossRef]
[16]	Zhang, L., Hamelers, H.V.M. and Biesheuvel, P.M. (2020) Modeling Permeate Ph in RO Membranes by the Extended Donnan Steric Partitioning Pore Model. Journal of Membrane Science, 613, Article ID: 118511.[CrossRef]
[17]	Siva Prasad Reddy, B.V. and Rajesh Babu, K. (2016) Application of Variational Methods and Galerkin Method in Solving Engineering Problems Represented by Ordinary Differential Equations. Proceedings of IRF International Conference, Goa, 28th February 2016, 36-37. https://www.digitalxplore.org/up_proc/pdf/206-145672797836-41.pdf
[18]	Nguyen, V.P., Rabczuk, T., Bordas, S. and Duflot, M. (2008) Meshless Methods: A Review and Computer Implementation Aspects. Mathematics and Computers in Simulation, 79, 763-813.[CrossRef]
[19]	Hansen, M.S. (2013) Boundary Conditions for 3D Fluid-Structure Interaction Simulations of Compliant Vessels (By Norwegian University of Science and Technology). https://leifh.folk.ntnu.no/alumni/2013_Merethe_Sirevaag_Hansen.pdf
[20]	Dougaj4 (2017) Stepping through 3DFrame. Newton Excel Bach, not (Just) an Excel Blog. https://newtonexcelbach.com/2017/09/16/stepping-through-3dframe/
[21]	Moradi, N., Liu, B., Iversen, M., Kuypers, M.M., Ploug, H. and Khalili, A. (2018) A New Mathematical Model to Explore Microbial Processes and Their Constraints in Phytoplankton Colonies and Sinking Marine Aggregates. Science Advances, 4, eaat1991.[CrossRef] [PubMed]
[22]	Lee, S. (2025) Numerical Solution of PDEs for Chemical Engineers. Number Analytics. https://www.numberanalytics.com/blog/numerical-solution-of-pdes-for-chemical-engineers
[23]	Kolar-Požun, A. AND Kosec, G. (2025) A Meshless Method for Computational Electromagnetics with Improved Dispersion Properties. arXiv: 2508.18205. https://arxiv.org/html/2508.18205v1
[24]	Okamoto, Y. and Lienhard, J.H. (2019) How RO Membrane Permeability and Other Performance Factors Affect Process Cost and Energy Use: A Review. Desalination, 470, Article ID: 114064.[CrossRef]
[25]	Hadadian, Z., Zahmatkesh, S., Ansari, M., Haghighi, A. and Moghimipour, E. (2021) Mathematical and Experimental Modeling of Reverse Osmosis (RO) Process. Korean Journal of Chemical Engineering, 38, 366-379.[CrossRef]
[26]	El-Brol, A., Ahmed, E., EL-Samadony, Y., Hammad, F. and Ismaeil, K. (2023) Experimental and Theoretical Study on Reverse Osmosis Based Water Desalination. Journal of Engineering Research, 7, Article 26.[CrossRef]
[27]	Jung, O., Saravia, F., Wagner, M., Heißler, S. and Horn, H. (2019) Quantifying Concentration Polarization—Raman Microspectroscopy for In-Situ Measurement in a Flat Sheet Cross-Flow Nanofiltration Membrane Unit. Scientific Reports, 9, Article No. 15885.[CrossRef] [PubMed]
[28]	Tayeh, Y.A., Alazaiza, M.Y.D., Alzghoul, T.M. and Bashir, M.J. (2025) A Comprehensive Review of RO Membrane Fouling: Mechanisms, Categories, Cleaning Methods and Pretreatment Technologies. Journal of Hazardous Materials Advances, 18, Article ID: 100684.[CrossRef]