Runge–Kutta Numerical Method Followed by Richardson’s Extrapolation for Efficient Ion Rejection Reassessment of a Novel Defect-Free Synthesized Nanofiltration Membrane

A defect-free, loose, and strong layer consisting of zirconium (Zr) nanoparticles (NPs) has been successfully established on a polyacrylonitrile (PAN) ultrafiltration substrate by an in-situ formation process. The resulting organic–inorganic nanofiltration (NF) membrane, NF-PANZr, has been accurately characterized not only with regard to its properties but also its structure by the atomic force microscopy, field emission scanning electron microscopy, and energy dispersive spectroscopy. A sophisticated computing model consisting of the Runge–Kutta method followed by Richardson extrapolation was applied in this investigation to solve the extended Nernst–Planck equations, which govern the solute particles’ transport across the active layer of NF-PANZr. A smart, adaptive step-size routine is chosen for this simple and robust method, also known as RK4 (fourth-order Runge–Kutta). The NF-PANZr membrane was less performant toward monovalent ions, and its rejection rate for multivalent ions reached 99.3%. The water flux of the NF-PANZr membrane was as high as 58 L · m−2 · h−1. Richardson’s extrapolation was then used to get a better approximation of Cl− and Mg2+ rejection, the relative errors were, respectively, 0.09% and 0.01% for Cl− and Mg2+. While waiting for the rise and expansion of machine learning in the prediction of rejection performance, we strongly recommend the development of better NF models and further validation of existing ones.


Introduction
With an ever-growing population and an increase in their standard of living and needs, as well as the expansion of industrial and agricultural activities, there is still an increasing demand for good quality water around the world. Moreover, worldwide, water scarcity is recognized as not only a present but also a future threat to human activities. To meet this increase in demand, with no loss being tolerable, a water treatment plant, in all its aspects and all its forms, needs to be optimized and sophisticated. One of the major technological challenges nowadays is the development of sustainable processes for water desalination, water reuse, wastewater treatment, and recovery of valuable chemicals from water [1][2][3]. Among membrane technologies, the pressure-driven process was reported to be effective for the separation of multivalent ions such as Mg 2+ , Ca 2+ , SO 2− 4 , etc. [3,4]. A nanofiltration (NF) membrane is a pressure-driven process exhibiting properties that lie betwixt ultrafiltration (UF) and reverse osmosis (RO) [5]. NF membranes have proven their effectiveness in several sectors such as water reclamation [6] and dye separation [7], heavy metals [8,9], pesticides [10], viruses and bacteria [11,12], beverages [13], natural organic matter [14], food [15], dairy processing [16], hardness removing [17], taste and odors [18], and even water softening [19]. Since there are several applications for NF membranes, there is, therefore, a compelling need to understand its separation behavior and how to improve the solute particles' transport mechanisms through its active selective layer, especially in water desalination. The NF process is reported to be extremely complex and dependent on both interfacial and micro-hydrodynamic events that occur at its surface, even within its nanopores. Donnan, steric, transport, and dielectric effects, combined together, determine the NF membrane's removal performance [20].
At the start, the dielectric exclusion phenomena were not understood very well with two competing hypotheses that tried obscurely to explain the exact nature of interaction that takes place. The first hypothesis was called the "image forces" phenomenon [21], while the second was called the "solvation energy barrier" mechanism [22]. These two exclusion phenomena have already been investigated in detail [23]. The selectivity of solute particles and the permeate flux are the two main factors that determine the membrane performance. Recent investigations demonstrated that the membrane must be highly chemical resistant, thermally stable, and loose to exhibit good separation results. NF membranes can be either symmetric or asymmetric in structure, homogenous or heterogeneous, and neutral or positively or negatively charged. The NF membrane separation process is not only faster but also more efficient and cost-effective than conventional separation techniques. The separation process using NF membranes has the following advantages: (i) energy consumption is low, (ii) separation can be performed under mild conditions, (iii) no additives are required, (iv) separation occurs continuously, (v) it is possible to combine NF process with other separation processes, (vi) the separation process can be up-scaled easily, and (vii) NF membrane properties can be adjusted.
In the present work, an organic-inorganic NF membrane has been synthesized to conciliate the tremendous advantages of both organic and inorganic membranes. A polyacrylonitrile (PAN) ultrafiltration (UF) membrane has been used as a substrate, and zirconia (ZrO 2 ) nanoparticles (NPs) have been chosen for deposition under the in situ formation process. Dopamine hydrochloric/sodium bicarbonate buffers are co-deposited on the PAN platform to induce the further growth of dioxide zirconium (ZrO 2 ) NPs on the membrane surface. The resultant organic-inorganic NF membrane, NF-PANZr, has been depicted as to its structure by various tools including atomic force microscopy (AFM), field emission scanning electron microscopy (FESEM), and energy dispersive spectrometry (EDS).
NF models based on the extended Nernst-Planck equations that govern the transport of solute particles through the NF membrane active layer have been imagined so far. A good model should include more of the complex phenomena that govern separation mechanisms to improve not only the physical robustness but also the relevance of the description of the process. Runge-Kutta method, a reasonably simple and robust method [24], has been used for ion rejection reassessment followed by Richardson's extrapolation to get a better approximation of solute particle rejection.

Mathematical Modeling
In numerical analysis, the Runge-Kutta methods include both implicit and explicit iterative methods. Runge-Kutta methods are based on the reputed routine called the Euler method, used for approximate solutions of ordinary differential equations. Also known as RK4 (fourth-order Runge-Kutta), this method is reasonably simple and robust and is highly recommended for all differential equations, provided that a smart, adaptive step-size routine is chosen [24].

(i) Boundary conditions:
It is assumed that Equations (5) and (9) can be solved over the following conditions: The solution understudy is ideal.
(iii) Each solute particle is subjected to an extended Nernst-Planck equation and could therefore be transportable. (iv) The effective charge density of NF-PANZr membranes does not change from one point to another on its surface. (v) The layer thickness of nanoparticles is assumed to be negligible toward the platform thickness. (vi) The NF membrane consists of an identical bundle of straight cylindrical pores, with each pore displaying a uniform depth and radius r p ∆x. (vii) The electric potentials inside the membrane and the Na 2 SO 4 , MgSO 4 , NaCl, CaCl 2 , and MgCl 2 solutions are all defined in terms of averaged quantities. (viii) The Donnan equilibrium is applied at both the interface of feed solution-membrane and the interface of membrane-permeate solution.

Focus on Model Equations
Runge-Kutta (RK) techniques were first introduced at the beginning of the 19th century by C. Runge and M. W. Kutta. Then, shortly after, this method took a major role in the study of iterative methods based on explicit, partial, and implicit methods that were applied to solve the ordinary differential equations (ODEs) using a temporal discretization.
Through his famous book, C. Runge had the ingenious idea to develop and make more accurate Euler's approximation method by proposing a scheme that could offer very high precision. Runge's work was first developed by Heun around 1900 and then by Kutta in 1901. The base of this scheme is expressed as follows: , w i are constants, and k i are given by This is the consistency relation of the Runge-Kutta (RK) method, ∑ s j=0 a ij = c i and ∑ s j=0 b j = 1.
The main equation to solve, which governs the solute particles' transport across the active NF membrane active layer, is known as the extended Nernst-Planck equation and is given as The transport of ions through the membrane can be achieved by applying the defined boundary conditions. It is easier to assess solute particle rejection by writing the Nernst-Planck equation as concentration and potential gradients. For the concentration gradient determination, the relation between the ion flux and its concentration is depicted as where C i,p is the ion [i] concentration in the permeate mol·m −3 . By substituting Equation (2) into Equation (1) and rewriting, the concentration gradient is given as Several conditions were involved in the potential gradient obtention. The electroneutrality within the pore, the permeate solution, and feed solution are governed respectively by Equations (6)-(8): where X d is the effective charge density of the membrane (in mol·m −3 ). The potential gradient is obtained taking into account the above-defined conditions in Equations (6)-(8) for the concentration gradient depicted in Equation (5). The potential gradient is then obtained by The equilibrium suggested by Donnan is naturally ensured by applying Equation (10) to the two important interfaces, NF-PANZr/permeate solution and feed solution/NF-PANZr. This equilibrium is ensured by Equation (10) below: z i , the valence of ion [i]; γ 0 i , bulk activity coefficient of ion [i], dimensionless; ∆Ψ D , Donnan potential variation (V); γ i , activity coefficient of ion [i] within pore; T, temperature; Φ i , steric partition coefficient. Furthermore, the steric partition coefficient of ion i is obtained by Equation (11): λ i , the ratio of ionic solute radius r i to NF-PANZr pore radius, r p . Considering an ideal condition, the steric partition has been removed from the Donnan equation. Assuming also that the solution is very dilute, the activity coefficient of ion [i] within the pore, to be taken into account inside the membrane by the effective charge density of the membrane, would be equal to 1. Equation (11) then becomes c i , ionic concentration within pore, mol·m −3 ; ∆Ψ D , Donnan potential variation (V); R, universal gas constant; F, Faraday constant; C i , ionic concentration in the solution mol·m −3 and z i , valence of ion i; The ion-i rejection (R) is given as follows: where C i, f , the concentration of ion [i] in the feed solution mol·m −3 , C i,p , ion-i concentration in the permeate mol·m −3 . Moreover, the ionic pore diffusion coefficient D i,p and the ionic hindrance factor for convection K i,c in the extended Nernst-Planck equation could be obtained respectively by Equations (14) and (15). The ionic pore diffusion coefficient, D i,p , is then obtained as below: K i,d , ionic hindrance factor for diffusion, dimensionless; D i,p , the pore diffusion coefficient of ion [i], in m 2 ·s −1 ; D i,∞ , ionic bulk diffusion coefficient (m 2 ·s −1 ). The hindrance factor for convection K i,c is given, taking into account the ionic velocity in the nanofiltration membrane (NFM) pores [22]: K i,c , ionic hindrance factor for convection; Φ i , ionic steric partition coefficient (Equation (11)); and G, hydrodynamic drag coefficient. The ionic hindrance factor for diffusion K i,d [25] is defined as follows: λ i , the ratio of the Stokes radius of element [i] to the pore radius of the membrane; G and K represent the hydrodynamic drag coefficients. For any solute particle [i], the ratio of ionic or uncharged solute to the pore radius, which is dimensionless can be obtained by Equation (19) below.
where r i is the Stokes radius and r p is the effective pore radius (membrane) of ion i. Equations (20) and (21) were finally obtained through a substitution of Equations (11) and (15)-(18) into Equations (5) and (9).
These last two equations define the variation in concentration inside the active layer of NF-PANZr as a function of the effective pore-radius, r p , and Stokes radius, r i , of the ion (i).

Description of the Computation Procedure
The internal concentration, c i,1 , of ion [i] is relative to the feed-solution concentration, C i, f , at the interface feed solution/membrane while the internal concentration, c i,N , of ion (i) is relative to the permeate concentration, C i,p , at the interface membrane/permeate solution. Runge-Kutta method was used to integrate both Equations (5) and (9) through the membrane active layer thickness Based on Equation (12), the knowledge of the value of C i, f , makes possible the integration of both Equations (5) and (9) after the determination of the initial concentration inside the NF-PANZr membrane c i,1 . (ii) Based on the Runge-Kutta numerical method, k 1 , k 2 , k 3 and k 4 and then c i,1 , c i,2 , c i,3 , c i,4 , . . . , c i,N could be well estimated (Equations (24)-(28); Equations (5) and (9)). (iii) Since the c i,N value is obtained, the permeate concentration, C i,p , was then computed. (iv) Lastly, Equation (13) was used to evaluate the ion [i] rejection.
The initial value for permeate solution, C i,p , is assumed equal to the feed-solution concentration, C i, f , which is the same as assuming zero rejection. The parameters such as hindered diffusivity, D i,p ; solute particle diffusivities; Stokes radii; the Donnan potential, ∆Ψ D ; partial molar volumes; and the hindrance factor for convection inside the membrane, K i,c can be found in the literature [26]. The membrane thickness, ∆x; molecular weight cutoff (MWCO); and the membrane pore size, r p ; were available from the nanofiltration membrane synthesized in the present study. Figure 1 shows the flowchart describing the Runge-Kutta modeling of NF membrane transport equation.

Ion Transport across NF-PANZr Membrane
Equation (5) according to the Runge-Kutta numerical method was given as follows: The step size is equal to the membrane thickness over the number (n) of steps, n = 100.
The active layer thickness (∆x = x 2 − x 1 ) of the novel organic-inorganic NF-PANZr is depicted in Figure 2 while Table 1 showed the typical model parameters required, h is the step-size, n is equal to 100. The Runge-Kutta fourth-order formula, sometimes known as RK4, is given as  The solute particle [i] concentration inside the membrane active layer varies from c i,1 at the feed solution/membrane interface side to c i,100 at the membrane interface/permeate solution side. Thereafter, c i,100 was used to calculate the permeate concentration, C i,p .
The program, in this case, was kept running until the deviation between the initial and final permeate concentration would be inferior to 10 −6 , we can finally evaluate this difference, for verification, by using the last formula:

Novel Organic-Inorganic Nanofiltration Membrane NF-PANZr Preparation
We assume the different steps below for the synthesis of the NF-PANZr membrane. The main steps for NF-PANZr preparation are well depicted in Figure 3. Step-1: Hydrolysis of polyacrylonitrile (PAN) membrane UF membranes of polyacrylonitrile were hydrolyzed in a sodium hydroxide,NaOH, solution (2.0 mol·L −1 ) for 2 h at 50 • C. The resultant membrane is H-PAN.
Step-2: Co-deposition of dopamine hydrochloride (DA) and sodium bicarbonate buffer (Buffer) A solution (S 1 ) was obtained from dopamine hydrochloride (DA) dissolved in a solution (pH = 8.0, 50 mmol·L −1 of sodium bicarbonate buffer (Buffer) for deposition. S 1 solution preparation is a result of in situ formation technology of DA-Buffer deposition on the PAN membrane surface in order to obtain a thin-film composite (TFC) layer acting like bio-glue. H-PAN membranes were pre-wetted using ethanol solution for 30 min before immersion into the solution (S 1 ) and stirred at 25 • C for 1 h 30 min. The resultant membranes (BG-PAN) were rinsed with pure water and dried in an oven set at 25 ± • C overnight.
Step-3: Deposition of zirconium (Zr) nanoparticles In this last step of novel organic-inorganic thin-film composite nanofiltration synthesis, another fresh solution, S 2 , is prepared by dissolving in hydrochloric solution (50 mmol·L −1 ), zirconium sulfate tetrahydrate with a concentration of 10 mmol·L −1 . Thereafter, the DA-Buffer coated H-PAN membrane pieces were immersed in the solution (S 2 ) at room temperature for 15 h. Following the deprotonation of the pyrocatechol groups ( C 6 H 4 (OH) 2 → C 6 H 4 (O−) 2 ) on the "bio-glue" layer, a covalent bond is established with Zr nanoparticles (as depicted in Figure 3) to generate a rigid layer of a complex compound of Zr. Finally, the resultant organic-inorganic NF-PANZr membranes were rinsed and dried in an ambient environment for later use and to be the subject of characterizations (properties, structures) and verification of performance.  The pore size distribution, effective mean pore size, and molecular weight cutoff (MWCO) of PAN, H-PAN, BG-PAN, and NF-PANZr membranes were measured by the solute particle rejection experiments using neutral rejection probes [27]. The feeds were glucose, sucrose, EG, α-CD, DEG, and PEG solutions of 150 ppm. Concentrations of organic solute particles were measured by a conductivity meter (Metrohm AG) and inductively coupled plasma optical emission spectroscopy (ICP-OES, Optima 7300 DV, PerkinElmer). The MWCO was calculated when the rejection reached 90%. The probability density function curve of the membrane was determined using the following Equation (31):

NF-PANZr Membrane Properties' Characterization
d p is the pore diameter. The mean pore size, µ p , is the pore diameter at which rejection, R, is 50% and the geometric standard deviation, σ p = d p (R = 84.13%)/d p (R = 50%).

NF-PANZr Membrane Structure Characterization
Field emission scanning electron microscopy (FESEM) has been used for investigation on the surface morphology of the novel thin-film composite NF membrane, NF-PANZr created in this work.
Atomic force microscopy (AFM, Multi-Mode VECCO, Denton, TX, USA) was used to observe the morphology and roughness of H-PAN, BG-PAN, and NF-PANZr membranes, and the results are shown in Figure 4.
Energy dispersive spectrometer (EDS) was used together with FESEM in order to get more insights on the NF-PANZr surface elements' arrangement, charge, and particularly about the phase state of the zirconium film. Zirconium (Zr), carbon (C), nitrogen (N), and oxygen (O) were among the different elements observed on the NF-PANZr surface.

Filtration Performance of Organic-Inorganic NF-PANZr Membrane
The performance of the novel organic-inorganic NF-PANZr has been investigated by a flat membrane module ensuring transverse flow at the laboratory scale under a pressure of 0.6 MPa, the temperature was set at 30 • C. Various classical salts Na 2 SO 4 , NaCl, MgSO 4 , CaCl 2 , and MgCl 2 solutions were prepared at a concentration of 1000 mg·L −1 , using DI water, and used as feed solutions at a constant cross-flow rate of about 30 L·h −1 and the effective surface area of NF-PANZr samples was 29.22 cm 2 . The permeate flux (J v , L·m −2 ·h −1 ) and chemical rejection (R, %) were obtained using Equations (32) and (33): A is the NP-PANZr effective surface, Q is the permeate-solution volume, and t is the time of permeation across the NP-PANZr layer.
C f and C p are, respectively, the solute particle concentration in the feed side and permeate side, C p and C f were determined by Metrohm AG, which is a conductivity meter, and another instrument, ICP-OES-Optima 7300 DV, PerkinElmer. All results presented were repeated at least three times, and the average values have been recorded and plotted.

Long-Term Stability of NF-PANZr Membrane
The NF-PANZr membrane was tested continuously for a whole week, and water flux and rejection readings were taken every 12 h and written down.

Richardson Extrapolation
Richardson's extrapolation was used to get a better approximation of Cl − and Mg 2+ rejection. The Richardson extrapolation of order n is given by the formula:

Statistical Error Analysis
Error analysis was performed through the agreement between experimental and model data. In the n data points for each solute, the least-squares fitting objective function RMSE is defined as follows:

NF-PANZr Properties' Characterization
The evaluation of the wettability of PAN, H-PAN, BG-PAN, and NF-PANZr membranes, performed with a Drop-Meter A200 contact angle system purchased from MAIST Vision Inspection & Measurement Co., Ltd., China, is shown in Figure 5a. The PAN membrane exhibited the best hydrophilicity since it possesses the most porous structure. The water contact angles of H-PAN and BG-PANZr were quite similar mainly from 60 s. The WCA has increased with the deposition of Zr nanoparticles (NPs). Although the novel organic-inorganic NF-PANZr membrane is the least hydrophilic of the four membranes tested in this section, it exhibits excellent hydrophilicity since the water contact angle is less than 72 • . Therefore, the water molecules could quickly spread over the NF-PANZr surface before seeping in. The pseudo-measurement of the Donnan potential was used to determine the membrane surface zeta potential (ζ) as shown in Figure 5b [28][29][30]. The state of PAN, H-PAN, BG-PAN, and NF-PANZr membranes' surface charge was determined for pH ranged between 3 and 10. BG-PAN membrane (bio-glue coated the PAN substrate) exhibited a negatively charged surface in contact with a basic solution (pH > 6), due to the deprotonation of the pyrocatechol groups ( C 6 H 4 (OH) 2 → C 6 H 4 (O−) 2 ) on its top layer surface. The NF-PANZr membrane surface is positively charged unlike the BG-PAN membrane surface as a result of the deposition of Zr NPs that in solution generate positive ions Zr 4+ . The NF-PANZr membrane is positively charged in the pH range used to carry out the experiments, namely 6 ± 0.5. Furthermore, the more or less attractive or repulsive effect of the NF-PANZr membrane on ion rejection was highlighted in the analysis of the results produced by the Runge-Kutta model. The charge of the membrane plays a very important role in the transport of ions through the membrane active layer, and several recent studies confirm this propensity of the membrane to attract and allow itself to be more easily crossed by ions of a charge contrary to its charge [8,[31][32][33][34].
The four membranes showed a negative charge above pH 9. The negative ions could then be easily and efficiently rejected for pH > 9 according to the Donnan effect. If the NF-PANZr membrane was used to remove ions at high pH (pH > 9), its behavior depicted above toward negatively charged ions will change. Instead of an attraction, we will expect a repulsion of these ions. Thus, for two ions, taken under the same stoichiometric conditions and of the same valence, the organic-inorganic NF-PANZr will exhibit better rejection toward the negatively charged ion. However, several parameters govern and condition the transport of ions.

NF-PANZr Structures' Characterization
The H-PAN, BG-PAN, and NF-PANZr membranes' surface morphology and roughness, which are valuable features in NF are depicted in Figure 4. The H-PAN membrane surface is porous ( Figure  4a) and accordingly, its roughness is the highest as shown in the AFM image (Figure 4d). The "bioglue" (dopamine hydrochloric-sodium bicarbonate buffer) coated substrate generated a smoother surface with almost invisible pores as depicted in Figure 4b,e. The hydrolyzed polyacrylonitrile membrane (H-PAN), once modified by deposition of zirconium nanoparticles displays the most smooth surface without visible pores (Figure 4c), while the membrane roughness is very small according to the AFM image shown in Figure 4e.
The membrane roughness decreased, to reach Ra = 5.17 nm, for the novel thin composite NF-PANZr membrane. The DA-Buffer coated the polyacrylonitrile substrate that acts like "bio-glue" has played an important role in zirconium NPs' deposition, having made it possible to obtain a more uniform and thinner membrane surface. Lv et al. have concluded from a recent study performed on a TFC NF membrane that a dense and smooth selective layer is beneficial for excellent rejection performance [35]. The NF-PANZr membrane may, therefore, provide an excellent solute particle removal and high permeation performance. This assertion will be verified later in this investigation.
To hold more insights into the arrangement of the Zr nanoparticle (NP) layer and its potential function on the membrane platform, the energy dispersive spectrometer (EDS) was associated with the field emission scanning electron microscopy (FESEM). Figure 6 shows the energy spectrum irradiation of elements on the NF-PANZr membrane surface. Four different elements, including nitrogen (N), carbon (C), oxygen (O), and zirconium (Zr) were displayed by the energy dispersive spectrometer image results, due to their unique X-ray signals. Zr, N, O, and C elements were distributed uniformly on the surface of the novel synthesized NF-PANZr membrane as shown in Figure 6b-e, respectively. It is therefore possible now to modify and characterize the materials at the atomic scale, providing unparalleled insight into the behavior of nanomaterials and particles, since each atomic position can be distinctly distinguished by its specific chemical signal (Figure 6a). The polyacrylonitrile used for the preparation of the NP-PANZr membrane is an organic membrane; carbon is, therefore, the most dominant element of the new membrane and consequently, a very strong signal is observed at the element C position. The sequence of the presence of the elements is C > O > Zr > N.   Table 2 as a function of copolymer content. The PAN substrate exhibits the largest pore diameter of about 20 nm. After the in situ formation process, the pore diameter of the resultant membrane (NF-PANZr) decreases dramatically to 0.4 nm. The hydrolysis step influenced the polymer very slightly since its pore diameter was around 18 nm. The observation is the same for the step of covering the membrane with DA/Buffer. On the other hand, the Zr nanoparticles significantly modified the membrane in terms of its structure. Figures 5 and 7 are a good match. The nanoparticles have contributed to significantly reduce the pores of the PAN membranes, which served as a platform (Figure 7). These nanoparticles are fine elements that have helped to close the "voids" of the PAN membrane, which, being an ultrafiltration membrane, is more porous. This resulted in a smoother structure with invisible pores (Figure 4).

Experimental Salt Rejection
The water flux and salt rejection of the NF-PANZr membrane are shown in Figure 8 for various salt solutions including NaCl, MgSO 4 , Na 2 SO 4 , CaCl 2 , and MgCl 2 at different applied transmembrane pressures, ∆Pe = 0.2, 0.3, 0.4, 0.5, 0.6, and 0.7 MPa. The water flux only increases with the pressure until it stabilizes around 58 L·m −2 ·h −1 under a pressure of 0.5 MPa. The thin-film composite NF-PANZr exhibited an excellent rejection performance for MgCl 2 , CaCl 2 , and MgSO 4 salts. These salts' rejection reached, respectively, 89.2%, 99.3%, and 95.3%. Na 2 SO 4 and NaCl salts were not efficiently removed by the novel synthesized membrane NF-PANZr. The salts' rejection sequence can be summarized as CaCl 2 > MgSO 4 > MgCl 2 > NaCl > Na 2 SO 4 . As a result, NF-PANZr has shown its ability to effectively reject multivalent ions such as Mg 2+ , Ca 2+ , and SO 2− 4 from water and its incapability to remove monovalent ions such as Cl − . This point of view is widely shared by a large number of publications recently performed [36][37][38]. Table 3 below provides an overview of the data used in the present model.   In both figures, the concentration of the divalent cation Mg 2+ in the permeate decreases dramatically for the highest water flux ( J v > 40 L·m −2 ·h −1 ), especially in the typical cases of number of steps = 200. Globally, the divalent cation Mg 2+ ions in the permeate decrease with increasing permeate flux J v (L·m −2 ·h −1 ), and this occurs independently of the number of incrementation steps. On the contrary, the rejection of Mg 2+ ions increased as the water flux (J v ) through the NF membrane increased. NF-PANZr rejection of Mg 2+ in the typical case of n = 200 was higher than that of n = 100. A smaller and more refined increment step, therefore, promotes the rejection efficiency. In the hypothetic case of J v = 58 L·m −2 ·h −1 (salt MgSO 4 ) the theoretical rejection of the Mg 2+ ion is, respectively, 95.1% and 91.8% for n = 200 and n = 100. Since the experimental rejection of Mg 2+ is 95.3%, the relative error made in the typical case where a high number of increments has been chosen is very closed to 0.02% (Table 4).

Runge-Kutta Model Reevaluation of Cl − and Mg 2+ Rejection
Moreover, Figures S1 and S2 depict the SO 2− 4 ion concentration inside the membrane active layer versus the step size for different water flux, J v .
NF-PANZr rejection of Mg 2+ was higher than the rejection of SO 2− 4 . Such rejection behavior is related to the membrane charge, which is of a positive charge at pH 6 ( Figure 5b). Figure 3 shows the deposition of the ion Zr 4+ on the membrane surface.   (6) and (8), the NF-PANZr charge effect appeared in the module in the electrical potential gradient and resulted in a disparate influence on SO 2− 4 and Mg 2+ removal. Repulsion between the NF-PANZr membrane charge and the Mg 2+ ions occurred while the attraction between the NF-PANZr membrane charge and the SO 2− 4 ions took place; in other words, SO 2− 4 ions would pass more freely across the membrane active layer, while the Mg 2+ ions would be pushed back. For all plots, it was noticed that the concentration of SO 2− 4 and Mg 2+ ions inside the membrane active layer decreased as the ions moved via the membrane active layer from the feed side to the permeate side.
Another parameter that affected SO 2− 4 and Mg 2+ ions' rejection is the membrane active layer pore size, r p , which appears in the hindrance factor for diffusion, K i,d , and the hindrance factor for convection, K i,c , as clarified in Equations (9), (14), and (15). In the typical case of a water flux, J v , of 58 L·m −2 ·h −1 , the concentrations of Mg 2+ ions across the NF-PANZr membrane can be determined for number steps (100 and 200) by the respective equations: Globally, the Cl − ions' rejection increased as the water flux, J v , via the NF-PANZr membrane increased. Since the thickness (x 2 − x 1 ) of the membrane has been divided by 100 (Figure 10b) and then by 200, the xstep has been set at 11.8 nm in the first case and 5.9 nm in the second case. NF-PANZr rejection of Cl − was higher for a small step size (Figure 10a) for all water fluxes, J v . In the hypothetic case of J v = 58 L·m −2 ·h −1 , the theoretical rejection of the Cl − ion is, respectively, 32.2% and 36.4% for n = 200 and n = 100. Since the experimental rejection is 31.9%, the relative error made in the typical case where a high number of increments has been chosen is small. To better analyze the behavior of the NF-PANZr membrane concerning the monovalent ions Cl − and Na + taken under the same stoichiometric conditions, Figures S3 and S4 show the plots of ion Na + model rejection for n = 200 and n = 100 respectively. NF-PANZr rejection of Na + was higher than the rejection of Cl − . Such rejection behavior is related to the membrane charge, which is of a positive charge at pH 6 ( Figure 5b). Figure 3 shows the deposition of the ion Zr 4+ on the membrane surface.
From Equations (6) and (8), the NF-PANZr charge effect appeared in the module in the electrical potential gradient and resulted in a disparate influence on Cl − and Na + removal. Repulsion between the NF-PANZr membrane charge and the Na + ions occurred, while the attraction between the NF-PANZr membrane charge and the Cl − ions took place; in other words, Cl − ions would pass more freely across the membrane active layer, while the Na + ions would be pushed back. For all plots shown in Figure 10a,b and Figures S3 and S4, it was noticed that the concentration of Cl − and Na + ions inside the membrane active layer decreased as the ions moved via the membrane active layer from the feed side to the permeate side. The Table 5 depicts the properties such as Stokes radii, diffusivities, and partial molar volumes of solute particles used in this investigation. Another parameter that affected Cl − and Na + ions' rejection is the membrane active layer pore size, r p , which appears in the hindrance factor for diffusion, K i,d , and the hindrance factor for convection, K i,c , as clarified in Equations (9), (14) and (15).
In the typical case of a water flux, J v , of 58 L·m −2 ·h −1 , the concentrations of Cl − ions across the NF-PANZr membrane can be determined for number steps = 100 and 200 by the respective equations: c i Cl − = − 1.0632 × x step + 5.2632; n = 100; c i Cl − = − 1.0706 × x step + 5.2632; n = 200.

NF-PANZr Long-Term Stability
The long-term stability of the novel organic-inorganic NF-PANZr membrane was investigated under a continuous filtration test for one week, and the results are shown in Figure 11. Of the various salts studied in this work, the long-term stability test was performed with MgCl 2 . The conditions under which the experiment was carried out were as follows: pH = 6.0, concentration [MgCl] 1000 mg/L, temperature = 30 • C, transmembrane pressure ∆Pe = 0.6 MPa, cross − flow rate = 30 L·h −1 . Globally, the thin-film composite NF membrane NF-PANZr has shown excellent long-term stability for permeate flux and salt removal. The water flux through the NF-PANZr membrane was almost constant over time for a slight decrease of less than 1.1%. The water flux was about 58 L·m −2 ·h −1 during the long-term operation and did not change remarkably till the end of the test. During the 168 h test, the novel organic-inorganic NF membrane exhibited an excellent rejection performance of 88.9% toward MgCl 2 salt. The rejection, like water flux, hardly varied over time, which is the advantage of this newly synthesized membrane. This one-week continuous test made it possible to appreciate the long-term effectiveness of the NF-PANZr membrane synthesized in this work. The good durability of the NF-PANZr membrane is closely related to the interfacial compatibility between the ZrO 2 nanoparticles (NPs) as a selective layer, the support surface polyacrylonitrile (PAN) membrane through the robust and multiple binding forces between the DA-Buffer ("bio-glue") coating and the H-PAN membrane.

Richardson Extrapolation and Statistical Error Analysis
The results predicted by the Runge-Kutta method are very far from the experimental results if the number of incrementations chosen is n = 100. However, these results are very good in the case where the number of incrementations is large, n = 200, and consequently produces a very small step size, Q app (h/2). Using Richardson's extrapolation, the theoretical results are found to be very refined and practically the same as the experimental ones.

Comparison of the Model Implemented in This Study with Other Previous Membrane Models
Most of the modeling-based research has been carried out on commercial membranes. The only reported work that combines a synthesized membrane and modeling, to our knowledge, is that of Farci et al. [39]. This famous work reported a Donnan steric-pore model with dielectric exclusion (DSPM-DE) model for both water flux and salt rejection prediction of microporous organosilica and mesoporous γ-alumina membranes, and scanning electron microscopy (SEM) was used to check the structure of the modified membrane. Table 6 provides an overview of some work on modeling. This work reports a novel NF membrane obtained from the deposition of Zr nanoparticles; the NF-PANZr membrane is fully characterized from a structural point of view (FESEM, EDS, AFM images), and its properties have been amply elucidated (water contact angle and zeta potential). The Runge-Kutta method supplemented by Richardson's extrapolation has been shown to be very effective in predicting ion rejection. This fact is understandable when we compare the numerical method of Euler (RE = 0.9%) to that of Runge-Kutta (the numerical method of Euler is the numerical method of first-order Runge-Kutta). The fourth order Runge-Kutta has been refined by the Richardson extrapolation. The nature and specificity of the project on which we are working require it.

Conclusions
The NF-PANZr membrane, like the NF membranes, has demonstrated excellent rejection performance against multivalent ions, and the ion rejection prediction model based on the Runge-Kutta method has been shown to be effective. In addition, the long-term stability carried out on NF-PANZr reassured that the membrane can be used for a long time before being replaced, especially if the membrane cleaning work consisting of its bombardment with drafts is periodically undertaken. The models are set up and calibrated for predefined input parameters, and they are not subject to change as chemical agents are changed on an experimental scale. We strongly recommend, therefore, that the scientific world study the distinctive characteristics of the next-generation NF membranes by incorporating the use of simple but effective models for more reliable results.
Further validation of existing models and the development of better NF models is a requirement for better NF membrane characterization, perfect chemical rejection reassessment, and even more insights on additional experimental available data. Beyond physical models, the rise of machine learning in the prediction of the rejection performance of NF membranes is also a horizon to explore.