Given its high surface area to volume ratio and desirable mass transfer characteristics, the hollow ﬁber module conﬁguration has been central to the development of RO and UF technologies over the past ﬁve decades. Recent studies have demonstrated the development of a novel class of low-pressure nanoﬁltration (NF)hollowﬁbermembraneswithgreatpromiseforscale-upimplementation. Furtherprogressonlarge-scale deployment, however, has beenrestrainedbythelack ofanaccuratepredictivemodel, toguide module design and operation. Earlier models targeting hollow ﬁber modules are only suitable for RO or UF. In this work, we propose a new modeling approach suitable for NF based on the implementation of mass and momentum balances, coupledwithavalidatedmembranetransportmodelbasedontheextendedNernst-Planckequation to predict module performance at the system-level. Modeling results are validated with respect to synthetic seawater experiments reported in an earlier work. A preliminary module design is proposed, and parametric studies are employed to investigate the eﬀect of varying key system parameters and elucidate the tradeoﬀs available during design. The model has signiﬁcant implications for low-pressure nanoﬁltration, as well as hollow ﬁber NF module design and operation. Keywords: nanoﬁltration; hollow ﬁber module; system-level modeling; seawater pretreatment in desalination; energy consumption Corresponding author: lienhard@mit.edu Preprint submitted to Desalination April 2, 2018O. Labban, T.H. Chong, J.H. Lienhard V, “Design and Modeling of Novel Low-Pressure Nanofiltration Hollow Fiber Modules for Water Softening and Desalination Pretreatment,” Desalination, 439:58-72, 1 August 2018. Nomenclature Roman Symbols A Debye–Hückel constant, m3/2/mol1/2 Ac Flow cross-sectional area, m2 Ak Membrane porosity ci Solute concentration, mol/m3 di Fiber inner diameter, m do Fiber outer diameter, m CPi Concentration polarization factor Di;p Diﬀusion coeﬃcient in the pore, m2/s Di;1 Diﬀusion coeﬃcient in the bulk, m2/s Df Fractal dimension Dh Hydraulic diameter, m Di;module Module inner diameter, m e Speciﬁc energy consumption, kWh/m3 e0 Elementary charge, 1:602181019C fD Darcy friction factor F Farady constant, 96485.3 C/mol I Ionic strength, mol/m3 Jv Permeate ﬂux, m3/m2s k Boltzmann constant, 1:380651023J/K kc;i Solute mass transfer coeﬃcient, m/s Ki;c Convection hindrance factor Ki;d Diﬀusion hindrance factor L Module length, m NA Avogadro’s number, 6:022141023mol-1 Nf Number of ﬁbers in the module _Ni Molar ﬂow rate of species i, mol/s Pw Wetted perimeter, m Q Volumetric ﬂow rate, m3/s rp Eﬀective pore radius, m ri Solute Stokes radius, m R Universal gas constant, 8.31446 J/mol K Rex Experimental solute rejection Ri Modelled solute rejection RR Recovery ratio 2 Re Reynolds number Sc Schmidt number Sh Sherwood number T Temperature, K V Channel bulk velocity, m/s Xd Membrane charge density, mol/m3 zi Ion valency Greek Symbols z Cell thickness, m P Applied pressure diﬀerence across the membrane, Pa  Osmotic pressure diﬀerence across the membrane, Pa WiBorn solvation energy barrier, J x Thickness of membrane active layer, m i Activity coeﬃcient i Ratio of solute Stokes radius to eﬀective pore radius  Solvent viscosity, Pa s Packing density i Ratio of permeate ﬂux to the uncorrected mass transfer coeﬃcient i Steric partitioning coeﬃcient B Born solvation coeﬃcient Electric potential, V f Correction factor for fractal dimension  Solution density, kg/m3 "0 Permittivity of vacuum, 8:854191012F/m "b Relative permittivity/dielectric constant of the bulk "p Relative permittivity/dielectric constant of the pore "r Relative permittivity/dielectric constant  Electric potential gradient at the feed/membrane interface, V/m i Mass transfer coeﬃcient correction factor Subscripts b bulk feed solution D Donnan Dh hydraulic diameter f feed in inlet 3 L hydraulic loss p permeate w condition at the wall Superscripts  Mass transfer correction for the suction eﬀect  optimal feed ﬂow rate for a given pressure 4 1. Introduction Nanoﬁltration (NF), introduced in late 80’s, is a membrane process whose performance falls between ultraﬁltration (UF) and reverse osmosis (RO); and as its name implies, has pore sizes on the order of 1 nm (corresponding to a molecular weight cut-oﬀ, or MWCO, of 300500Da) [1]. With their unique selectivities and high permeabilities relative to RO, NF membranes presented a major milestone in membrane technol- ogy, have attracted considerable attention since their introduction, and have found numerous applications, spanning numerous ﬁelds from water and waste water treatment to biotechnological, pharmaceutical, and food industry applications [2]. A recent review by Mohammad et al. [2] identiﬁes NF as a research domain of surging interest, predicts NF will continue to play a prominent role in membrane technology, and reveals the future prospects and areas of potential growth NF is likely to experience in the long term. Recent studies have demonstrated the development of a novel class of low-pressure NF hollow ﬁber membranes with great promise for scale-up implementation [3, 4, 5]. The design of hollow ﬁber membrane modules based on those membranes is the focus of the current work. Given its high selectivity for monovalent ions, areas in which NF is projected to grow include hard water softening, and more recently, desalination pretreatment. While traditional softening technologies, such as lime softening and pellet softening, have been challenged by their large sludge production and chemical con- sumption, NF has emerged as a viable alternative [6, 7]. Apart from traditional softening, NF desalination pretreatment has lately been gaining attention as evident from a recent review on integrated/hybrid mem- brane processes in desalination and water treatment by Ang et al. [8]. The review cites evidence that NF desalination pretreatment not only improves desalination feeds by reducing scaling, but also allows seawater reverse osmosis (SWRO) to run at lower pressures and potentially achieve higher recoveries [8]. In search of the optimal conditions for NF/SWRO desalination, Park et al. followed with their study on NF pre- treatment and its eﬀect on SWRO recovery under diﬀerent conditions, and concluded that the NF/SWRO conﬁguration featured an improvement in the quality of recovered water, especially at high recoveries [9]. More recent works from Roy et al. [10, 11] explore the eﬀect of temperature on NF pretreatment and scale prevention in thermal desalination. As becomes apparent from Bergman’s cost analysis of softening technologies, however, one primary lim- itation on the economic viability of the NF pretreatment system has been the additional energy penalty incurred by the process [7]. Consequently, the requirement of lower energy consumption has occupied mem- brane researchers and process designers alike in search of better separation eﬃcacy. While most NF mem- branes have been thin ﬁlm composite (TFC) ﬂat sheet membranes, a novel class of low-pressure NF hollow ﬁber membranes developed using layer-by-layer (LbL) polyelectrolyte deposition with chemical cross-linking appears to provide a solution [3, 4, 5]. According to performance tests reported in earlier work [4, 5], the newly developed membranes fared well when compared to state-of-the-art in softening and possessed superior softening performance compared to commercial NF membranes while operating at pressures that did not exceed 5 bar. These results prove the membrane’s potential, raising research interest for potential large-scale 5 implementation. Despite their great potential, further progress on scale-up implementation of these or other NF hollow ﬁber membranes is restrained by the lack of information on large-scale or commercial NF hollow ﬁber modules, which is necessary to evaluate a membrane’s potential for scaling-up. In this work, we develop a mathematical model to our knowledge the ﬁrst to predict the performance of NF hollow ﬁber modules on the system-level, building from experiments run on a bench-scale setup. In an earlier study, we demonstrated the successful implementation of a membrane transport model, introduced by Geraldes and Alves [12] based on the extended Nernst-Planck equation, to the newly developed LbL hollow ﬁber membranes [13]. However, that model focused on coupon-sized systems, and did not take into consideration the streamwise variations that will be inherent in large-scale applications. To address this need, a new model, taking a deeper look into the fundamentals governing ﬂuid ﬂow and mass transfer in hollow ﬁber modules, becomes necessary. Since they were ﬁrst patented in 1966 by DOW Chemical Company [14] followed by DuPont [15], research and development on hollow ﬁber membrane modules has made signiﬁcant progress. These modules became attractive from an application standpoint given their desirable mass transfer characteristics and high surface areatovolumeratio,whichcouldreasonablyreach10timesthatofaspiral-woundmodule[16]. Consequently, hollow ﬁber modules have found applications across a broad spectrum of ﬁelds, ranging from desalination and water treatment to gas separation, hemodialysis, blood fraction, and hydrocarbon separation [17]. The ﬁrst attempts to model these systems were reported in the early 1970’s with the works of Hermans [18] and Gill and Bansal [19] on RO hollow ﬁber modules. These authors developed their models based on ﬁrst principles, namely two dimensional mass and momentum conservation equations. From there, later work focused on incorporating a diﬀusion transport model [20] and accounting for concentration polarization [21]. More recent works have targeted the nonuniformity and randomness that characterize hollow ﬁber modules [22], proposing approaches to modeling them [23, 24], such as the Voronoi tessellation method [25]. While this method presented a signiﬁcant improvement, its underlying assumption of a uniform probability function for all modules makes it diﬃcult to capture the varying degrees of packing irregularity possible under the same packing density [26]. Zhang, more recently, introduced a method based on fractal theory to approach this problem [26]. These hollow ﬁber modeling eﬀorts, nonetheless, still remain unsuitable for NF, whose performance requires radically more challenging approaches. Unlike RO, NF performance is governed by the interplay of diﬀerent separation mechanisms, namely sieving eﬀects, Donnan partitioning, and dielectric exclusion. As noted by Oatley-Radcliﬀe et al. [27], modeling diﬃculty is compounded with the introduction of multi-ionic mixtures, which are of particular interest to softening and desalination applications. Following the approach ﬁrst proposed by Tsuru et al. [28], the majority of NF models in literature are based on the extended Nernst- Planck equation. Having been developed for bench-scale setups, however, these models’ primary limitation has been their inability to predict system-level performance or capture streamwise variations, which can be detrimental to module eﬃcacy. In response, a previous study by Roy et al. [29] targeted the development 6 of comprehensive large-scale models for ﬂat-sheet and spiral-wound NF modules, but did not consider NF hollow ﬁber modules. Furthermore, the feasibility of operating at extremely low applied pressures on the large-scale is not backed by empirical evidence, and still appears questionable. Apart from the NF hollow ﬁber system intro- duced by Pentair for humic substances removal, a commercial scale NF hollow ﬁber module for desalination and softening applications is not currently available, nor is a well-established protocol to design these sys- tems. The problem is further exacerbated by the lack of a practical approach to assess the viability of a speciﬁc hollow ﬁber membrane for such applications, and evaluate the eﬀect diﬀerent module dimensions and operating parameters can have on performance absent the need for expensive experimental procedures. These challenges combined underscore the need for further investigation in this area. In this work, we present the development of a system-level model for NF hollow ﬁber modules, combining state-of-the-art NF modeling, based on the extended Nernst-Planck equation, with fundamental conservation laws that govern ﬂuid ﬂow and mass transfer in hollow ﬁber modules. The newly developed LbL1.5C mem- brane by Liu et al. [4] is adopted in this study, while the model presented is validated against experimental results reported in a previous work [13]. To demonstrate the model’s capabilities, a preliminary module design was proposed, and parametric studies were run to assess the viability of the LbL1.5C membrane for large-scale seawater desalination pretreatment. An approach based on comparing the marginal increase in recovery to the marginal rise in applied pressure is employed to illustrate that raising applied pressure can improve energy eﬃciency in some conditions when operating under extremely low applied pressures, and a preliminary procedure for module sizing, inspired from the insights developed in this work, is proposed. 2. Theoretical Background Generally, hollow ﬁber modules come in two primary conﬁgurations, the shell-side feed and the bore-side feed conﬁgurations [16]. In the shell-side feed systems, similar to those developed by DuPont for RO, the pressurized process stream ﬂows on the shell-side, as the permeate seeps through the membrane and ﬂows on the inside of the ﬁbers. Such a setup requires the incorporation of ﬁbers with a relatively smaller diameter to sustain the large hydrostatic pressures imposed [16]. The bore-side feed conﬁguration, on the other hand, features a process stream ﬂowing inside the ﬁbers with both ends open, while the permeate ﬂows on the shell-side, respectively. Modules of this conﬁguration typically have ﬁbers of larger diameters to minimize pressure drops in the feed channel [16]. This study will focus on modeling the bore-side feed conﬁguration, shown in Fig. 1, which was also the conﬁguration adopted by Liu et al. [3, 4]. According to Baker [16], concentration polarization is better controlled in this conﬁguration as the occurrence of stagnant or dead zones, and feed channeling is avoided. It should be noted, nonetheless, that the small ﬁber diameters encountered in practice could necessitate pretreatment, depending on the feed quality, to avoid ﬁber clogging. Shell-side feed conﬁgurations, on the other hand, are more prone to irreversible fouling as suspended particles ﬁnd themselves trapped in dead 7 zones [16]. Furthermore, shell-side feed conﬁgurations often suﬀer from signiﬁcant parasitic losses in the permeate channel, which could considerably lower the driving force, and hence the permeate produced, across the module [16]. Feed Retentate Permeate Figure 1: Schematic representation of the bore-side feed hollow ﬁber module conﬁguration, featuring an inside-out ﬁltration mode. In this section, the developed model and its underlying assumptions are introduced and justiﬁed. Mass and momentum balances are applied, and a method to calculate the permeate ﬂux as a function of the net driving force is proposed. The result is a system of conservation equations, while being strictly valid under the assumptions invoked, is still under-deﬁned. As a result, we then introduce a membrane transport model to close the system and complete the analysis. Next, metrics to evaluate the module performance at the system-level are discussed. Eventually, a preliminary design for a NF hollow ﬁber module is proposed. 2.1. System-Level Modeling: The ﬁrst step in modeling a NF hollow ﬁber module involves analyzing it at the system-level through the application of mass and momentum balances. For the system at hand, feed water enters the module, where it is guided along the inside of the ﬁbers as Fig. 1 illustrates. The result is an eﬀective splitting of the incoming feed stream among the numerous ﬁbers in the module. Under the assumption that all ﬁbers are identical and independent, it follows that the process stream is divided equally among the ﬁbers. Within each ﬁber, the feed stream is processed in an inside-out ﬁltration mode, such that the permeate is collected on the module’s shell-side, respectively. Once the feed water reaches the end of the module, its concentration would have already increased, leaving the system as retentate. 2.1.1. Conservation Laws: To model the streamwise variations in the system, the ﬁber is divided into cells in the axial direction, as shown in Fig. 2a, and a backward ﬁrst-order ﬁnite-diﬀerence scheme is implemented with a uniform grid. The boundaries of each cell form control volumes, illustrated in Fig. 2b, on which fundamental conservation laws, such as mass conservation, are applied to yield: Qf=Qp=Jv di z (1) 8 where QfandQpare the ﬁnite changes in the feed and permeate volumetric ﬂow rates across a cell, Jv is the permeate ﬂux, diis the ﬁber inner diameter, and zis the step size or cell thickness. 5 Rejection of a Salt and/or Salt Mixture 5.1 Data Requested: The rejection of the solution of modeling interest is required. 5.2 Justiﬁcation: The membrane charge, the only left parameter, is a function of the feed solution chemistry and is generally a result of functional groups’ dissociation and/or ion-adsorption. The DSPM&DE model does not account for these phenomena, respectively. For this reason, the model will need to be calibrated again for diﬀerent ion compositions to ensure accuracy. Althougha single rejection data point is suﬃcient to obtain the charge, providing more points (as a function of ﬂux/∆Pe) is preferable and aids in validation later on. z(2) ∆z(3) 5.3 Sample: For reference, a sample plot is provided in Figure 4. 6 Data Requested for 1D Model Validation To validate the 1D model developed at MIT, experimental results from the actual large scale modulesare requested. In particular, experimental results showing the spatial variation of module perfor-mance (rejection, ﬂux, concentrations) across its length is required. Placing NF modules in stagesmight help make any changes more visible, respectively. 55 Rejection of a Salt and/or Salt Mixture 5.1 Data Requested: The rejection of the solution of modeling interest is required. 5.2 Justiﬁcation: The membrane charge, the only left parameter, is a function of the feed solution chemistry and is generally a result of functional groups’ dissociation and/or ion-adsorption. The DSPM&DE modeldoes not account for these phenomena, respectively. For this reason, the model will need to becalibrated again for diﬀerent ion compositions to ensure accuracy. Althougha single rejection data point is suﬃcient to obtain the charge, providing more points (as a function of ﬂux/∆P e) is preferable and aids in validation later on. z(2) 5.3 Sample: For reference, a sample plot is provided in Figure 4. Figure 4: Rejection of Ternary System 5 FeedRetentate Permeate Feed Permeate(a) (c) (b)Progress Report 4 - Picking up Omar March 25, 2015 1 Questions Posed Last Presentation ci,p (1) 1 Figure 2: System-to-local level modeling: Modeling approach as applied to a single ﬁber to predict its separation performance. (a) The ﬁber is axially divided into cells; (b) conservation laws are applied across these cells; (c) a membrane transport model is coupled to the conservation laws, allowing for the permeate concentrations, ci;pto be determined. Equation 1 has been derived under the assumption of constant density, such that mass conservation becomes equivalent to volume conservation. The validity of this assumption as it pertains to seawater is backed by empirical evidence showing that the density changes by less than 5% for reference salinities varying from 0 g/kg to 50 g/kg [30, 31]. The assumption is even more appropriate for NF membranes, for which the density change is expected to be lower since primarily multivalent ions are retained by the membrane. Similar to Eq. 1, conservation of species is applied to each solute in the system to arrive at: _Ni;f=_Ni;p=Jvci;p di z (2) where _Ni;fand_Ni;pare the ﬁnite changes in the feed and permeate molar ﬂow rates of solute iaross a cell, andci;pis the solute concentration in the permeate. 2.1.2. Permeate Flux and the Driving Force: To successfully apply these equations, however, establishing the link between the permeate ﬂux, Jv, and the net driving force across the membrane becomes necessary. This is accomplished using an empirical expression, introduced earlier by Bowen et al. [32, 33] for NF membranes, that relates the ﬂux to the eﬀective 9 driving force, P, and membrane permeability: Jv=r2 p 8xe P (3) whererpis the membrane eﬀective pore size, is the solution viscosity, xeis the eﬀective membrane thickness deﬁned in line with Geraldes and Alves [12] as the ratio of membrane thickness to porosity x=Ak (Section 2.2), Pis the applied pressure diﬀerence, and is the osmotic pressure diﬀerence across the membrane. The Staverman constant is set to unity following Bowen et al. [33]. A detailed discussion of these parameters and their estimation has been covered in detail in an earlier work [13], as well as previous publications by Bowen and coworkers [32, 33, 34]. One major challenge imposed by Eq. 3 is expressing Pandas functions of position across the module. By accounting for the hydraulic losses, the applied pressure diﬀerence can completely be expressed as: P(z) =Pf(z)Pp(z) (4) Pf(z) =Pf;inPf;L(z) (5) Pp(z) =Pp;inPp;L(z) (6) wherePf;inandPp;inare the inlet pressures, and Pf;L(z) Pp;L(z)are the hydraulic losses on the feed and permeate sides, respectively. Under steady-state operating conditions and system speciﬁcations of practical interest ( di=L 1and L=Vf 1), the ﬂow in the ﬁbers is expected to be laminar, and can be approximated as locally fully- developed [35]. Under this approximation, the pressure drop on the feed side can be expressed according to the Darcy–Weisbach equation [35, 36]: Pf;L(z) =fDz di Vf(z)2 2 (7) fD=64 Re(8) Re =Vfdi (9) withfDbeing the Darcy friction factor, the solution density, Vf(z)the feed channel bulk velocity, and Rethe ﬂow Reynolds number. The pressure drop on the shell-side, on the other hand, is more diﬃcult to account for given the possibility of random ﬁber packing. According to the results reported by Zhang, who based his derivations on fractal theory, the shell-side pressure drop and the Darcy friction factor for modules with non-uniform ﬁber packing is evaluated according to the correlations [26]: Pp;L(z) =fDz Dh Vp(z)2 2 (10) 10 fDReDh= (428:62 3+ 620:49 2304:24 + 107:9) f (11) =Nf di Di;module!2 (12) f= 1:149Df1:0 (13) in whichVp(z)is the permeate ﬂow velocity evaluated as the ratio of the permeate ﬂow rate to the cross- sectional area of the permeate channel, is the packing density, Dfis the fractal dimension, and fis a correction factor that accounts for the randomness in ﬁber distribution at a given packing density. According toZhang,Dfapproaches1forirregularpacking,and2forregularpacking; itrangesfrom 1:61:9forpractical modules, respectively [26]. A value of Df= 1:7has been selected herein. Dhrefers to the hydraulic diameter of the shell-side, and is evaluated in tandem with the deﬁnition used by Lipnizki and Field [37]: Dh=4Ac Pw=D2 i;moduleNfd2 o Di;module +Nfdo(14) withAcbeing the ﬂow cross-sectional area, Pwthe wetted perimeter, Di;modulethe module inner diameter, andNfthe number of ﬁbers in the module. AnotherchallengeimposedbyEq.3isthecouplednatureoftheproblem. Thepermeateﬂux, asexpressed in Eq. 3, is a function of the net driving pressure; and as a result, will vary with the shell-side pressure. On the other hand, the shell-side pressure is also a function of the permeate ﬂux. A greater permeate ﬂux will lead to larger velocities in the permeate channel, resulting in greater pressure drops on the shell-side according to the Darcy–Weisbach equation. By the same token, a greater permeate ﬂux will lead to smaller velocities in the feed channel by conservation of mass, resulting in smaller pressure drops on the feed-side. While it could be argued a priori that the pressure drop on the permeate side will be negligible at low permeateﬂowrates[29], thevalidityofthisassumptionbreaksdownforhollowﬁbersystems, operatingunder considerably larger packing densities or high recovery ratios. In this work, we propose an iterative approach, depicted in Fig. 3, to evaluate the pressure on the shell-side when simulating the module performance. The process starts with an initial guess of the inlet pressure on the shell-side. From there, the model is solved sequentially for the permeate ﬂux and ﬂow rate in each cell, allowing for an estimation of the pressure drop on the shell-side. The estimated pressure drop is then compared to the initial guess, and the process is repeated until convergence. Using this iterative approach, we later show that the shell-side pressure drop is smaller than the feed-side pressure under most operating conditions, and can reasonably be neglected, greatly simplifying the analysis. 11 Pp,in,guess=Pp,inelseexitFindPp,in=fnℓDh,f,ρ,Vp()Pp,in,guessInitializeQpVpFindandJv=fnΔP−ΔΠ()Findε=Pp,in−Pp,in,guess≤εsetFigure 3: Algorithm employed to estimate the shell-side pressure drop. In addition to P, the osmotic pressure diﬀerence, , needs to be evaluated. In a work on modeling RO membranes as mass exchangers, Banchik et al. [38] demonstrated the validity of van ’t Hoﬀ’s equation for salinities in the range of 070g/kg (typical of seawater desalination). The maximum deviation observed from the nonlinear osmotic pressure model was 6.8% [38]. According to this equation, the osmotic pressure may be evaluated as [39]:  =  fp=RTNX i=1(ci;fci;p) (15) where the concentrations ci;fandci;prefer to an ion’s concentration in the feed and the permeate for a given salt mixture. Since the osmotic pressure behaves as a colligative property, the equation should still hold for slight vari- ations from seawater concentrations and compositions expected in nanoﬁltration. NF membranes generally, and the LbL membrane speciﬁcally, are mostly selective towards monovalent ions. Since the concentration of multivalent ions is minor in the solution, the deviation from seawater concentrations caused by their retention should also be minor, justifying the adoption of van ’t Hoﬀ’s equation herein. 2.2. Local-Level Modeling: Membrane-Transport Model While Eq. 15 completes our approach to system-level modeling, it presents yet another challenge by requiring the permeate concentration of each species, ci;p, to be known. Since ci;pis a function of the membrane’s selectivity, a membrane transport model is required. In this section, the choice of the transport model is addressed along with an account of the theory and its underlying assumptions. 2.2.1. Membrane-Transport Model: Model Selection An earlier work by Oatley-Radcliﬀe et al. [27] highlighted state-of-the-art NF modeling as it relates to seawater desalination, and how interest has shifted from models based on irreversible thermodynamics to those based on the extended Nernst-Planck equation. Among those models, Geraldes and Alves [12] introduced Nanoﬁltran , an open-source program for modeling transport of multi-ionic solutions through NF membranesbasedontheextendedNernst-Planckequation. Geraldesetal.[12]validatedtheprogramagainst experimental results previously reported by Bowen and Mohammad [32], and analyzed its stability under a stiﬀ NF system composed of several ions. Several studies in the literature, such as the works of Shah et al. [40] and Deon et al. [41], have successfully implemented a similar model, noting the agreement between experiments and simulations. Furthermore, Roy 12 et al. successfully implemented Nanoﬁltran to predict the performance of NF ﬂat-sheet and spiral-wound conﬁgurations at the system-level [29], and to study the eﬀect of temperature on ion transport [10]. Similarly, Epsztein et al. [42] incorporated Nanoﬁltran in their work investigating the role of ionic charge density on Donnan exclusion. More importantly, however, the program has been demonstrated to successfully predict the performance of the novel LbL NF membranes developed by Liu et al. [5, 13], whose scaling-up is the focus of the current work. Accordingly, Nanoﬁltran is adopted herein as the transport model of choice. In the following section, the theory underlying Nanoﬁltran is brieﬂy presented. For more details, readers are referred to the original work by Geraldes and Alves [12]. 2.2.2. Nanoﬁltran: Theoretical background Generally, the performance of thin ﬁlm composite (TFC) membranes is primarily dictated by the mem- brane active layer, allowing us to neglect the membrane support layer [33]. According to the model, the ﬂux of theithsolute through the membrane active layer is governed by the extended Nernst-Planck equation with hindered transport [43]: Ji=Di;pdci dx+Ki;cciJvziciDi;pF RTd dx(16) Di;p=Ki;dDi;1 (17) In this equation, which accounts for solute transport by diﬀusion, convection, and electromigration, Di;pis the solute diﬀusion coeﬃcient in the pore, Di;1is the solute diﬀusion coeﬃcient in the bulk solution, is the electric potential, ziis the ion valency, Fis the Faraday constant, Ris the universal gas constant, and Tis the temperature. Hindered transport in the pores is accounted for through the incorporation of the hindrance factors for diﬀusion, Ki;d, and convection, Ki;c, expressed as a function of the ratio of solute Stokes radius to pore size, i[44]. A recent work by Oatley-Radcliﬀe et al. [45] reviewed available correlations for the hindrance factors, and experimentally validated them for NF processes. Fori0:95, the result by Dechadilok and Deen can be invoked to determine Ki;d[44]: Ki;d=1:0 + (9=8)iln(i)1:56034i+ 0:5281552 i+ 1:915213 i 2:819034 i+ 0:2707885 i1:101156 i0:4359337 i (1i)2(18) Fori>0:95, we use the result proposed by Mavrovouniotis and Brenner [46]: Ki;d= 0:9841i i5=2 (19) The hindrance factor for convection, Ki;c, is determined in a similar fashion using the result of Ennis et al. [47]: Ki;c=1 + 3:867i1:9072 i0:8343 i 1 + 1:867i0:7412 i(20) 13 Although the membranes modeled in this study are cylindrical, Cartesian coordinates can still be used in Eq. 16 with reasonable accuracy provided that xe=do 1, with xebeing the eﬀective thickness of the membrane active layer and dothe ﬁber outer diameter [48]. Another important constraint to satisfy inside the active layer is charge electroneutrality, taking membrane charge density Xdinto consideration: Xd+NX i=1zici= 0 (21) Apart from modeling solute transport within the membrane active layer, a model is required to de- scribe transport across the concentration polarization layers at the membrane-feed/permeate interfaces. For membranes featuring bulk ﬂuid ﬂow through them, concentration polarization on the permeate side can reasonably be neglected [16]. To model concentration polarization on the feed side, however, the model developed by Geraldes and Afonso has been adopted in this study [49]. This model, which expresses the ﬂux of theithsolute through the concentration polarization layer as the sum of the ﬂuxes due to back diﬀusion, convection, and electromigration as shown in Eq. 22, accounts for concentration polarization through the incorporation of a mass transfer coeﬃcient as follows: Ji=k c;i(ci;wci;b) +Jvci;wzici;wDi;1F RT (22) k c;i=kc;i (23)  = i+ (1 + 0:26 1:4 i)1:7(24) In Eq. 22,ci;wis the solute concentration at the membrane-feed interface just outside the pores, ci;bis the solute bulk concentration, and is the electric potential gradient at the feed-membrane interface. More importantly, k c;irefers to the mass transfer coeﬃcient, estimated based on conventional Sherwood number correlations and corrected for the “suction eﬀect” resulting from the permeation at the membrane interface. This correction is achieved through the incorporation of a ﬂux-dependent correction factor, , deﬁned in Eq. 24 with i=Jv=kc;i[49]. To calculate the mass transfer coeﬃcient, kc;i, a Sherwood number correlation, inferred for the ﬂow inside hollow ﬁber membrane modules and the test conditions of interest, has been adopted [23]: Shi= 1:62Gz1=3 i= 1:62Re1=3Sc1=3 i(di=L)1=3; (Gzi>4) (25) Shi=kc;idi Di;1(26) Sci= Di;1(27) where Gziis the solute Graetz number, Shiis the solute Sherwood number, Reis the ﬂow Reynolds number, Sciis the solute Schmidt number, and Lis the module length, respectively. 14 InadditiontoEq.22, anotherequation, Eq.28, relatingtheﬂuxofthe ithsolutetoitsconcentrationinthe permeate can be obtained from the solute continuity equation under steady state operating conditions. Two more electroneutrality conditions, expressed as Eqs. 29 and 30, should also be satisﬁed on the membrane- feed/permeate interfaces, respectively [12]. Ji=Jvci;p (28) NX i=1zici;w= 0 (29) NX i=1zici;p= 0 (30) Once solute transport in the active layer and the concentration polarization layer has been considered, the ﬁnal step in setting up the membrane transport model is accounting for solute partitioning at the membrane-feed/permeate interfaces at chemical equilibrium. The result is two equations that act as “bound- ary conditions” on the membrane-feed and membrane-permeate interfaces, and are expressed as Eqs. 31 and 32 [12]: i;1ci;1 ( i;wci;w)= iBexp ziF RT D;w (31) i;Nci;N ( i;pci;p)= iBexp ziF RT D;p (32) In Eqs. 31 and 32, Drefers to the Donnan potential forming across the membrane at equilibrium [50, 51], ci;1andci;Nare the concentrations just inside the pores at the feed and permeate interfaces [13], while i refers to the solute activity coeﬃcient, calculated in Nanoﬁltran based on the Davies model, which relates ito the solution ionic strength, I, according to [52, 53]: ln( i) =Az2 ip I 1 +p IbI (33) I=1 2NX i=1z2 ici (34) withbassigned a value of 0:3, andAbeing the temperature-dependent Debye–Hückel constant evaluated according to [53, 54]: A=p2NA ln(10)e2 0 4"0"rkT3=2 (35) whereNAis Avogadro’s number, e0is the elementary charge, "0is the vacuum permittivity, and kis the Boltzmann constant. We note that an assumption concerning the validity of the Debye–Hückel theory inside a nanopore is being made herein. The Debye–Hückel theory assumes a bulk solution with the number of anions and cations being equal, which may be challenged inside a pore of nanometer dimension [33, 55]. Another alternative 15 proposed by Bowen and Welfoot [33] is to assume ideal solution behavior with activity coeﬃcients set to unity. Given the focus of this work on large-scale modeling, however, we choose to retain the assumptions laid previously by Geraldes and Alves [12] for consistency. Furthermore, the term iin Eqs. 31 and 32 refers to the steric partitioning coeﬃcient that accounts for the sieving eﬀects that result from the ﬁnite size of the solute relative to the pore [34]. On the other hand, the term Brefers to the Born solvation coeﬃcient that accounts for the solute partitioning that results from dielectricexclusion, anditsassociatedbarriertoionsolvationarisingfromthediﬀerenceindielectricconstant between the membrane pores and the bulk solution [56, 57]. These two terms are evaluated according to Eqs. 36 and 37 [12]: i= (1i)2(36) B= exp Wi kT (37) Wi=z2 ie2 0 8"0ri1 "p1 "b (38) where Wiis the solvation energy barrier computed based on Born’s model [58], riis the solute’s Stokes radius, and "bis the dielectric constant of the bulk solution. Together, Eqs. 31and32accountforsolutepartitioningoccurringasaresultofsterichindrance, dielectric exclusion, and Donnan partitioning, and completely deﬁne the problem of modeling solute transport. The next steps involve discretizing the membrane’s active layer, applying one extended Nernst-Planck equation for every solute at every discretization point, coupling these equations with the concentration polarization equation and boundary conditions, and ﬁnally solving the resulting system of equations numerically [13, 12]. The outcome of this computation will include the unknown permeate concentrations, ci;p, closing the system of equations as Fig. 2c illustrates. 2.3. System-to-Local Level Modeling Figure 2 summarizes the modeling approach presented in this work, and illustrates how the membrane transport model is coupled to the system-level analysis to provide a comprehensive model that can accurately predict the performance of a NF hollow ﬁber module under diﬀerent operating conditions. The ﬁgure also reveals the necessary implementation of a two-dimensional discretization scheme (parallel + perpendicular to the ﬁber axial direction) for every solute, resulting in a three-dimensional matrix of unknowns, to solve the problem. The solver in this work was implemented using MATLAB. The integration of a membrane transport model over the hollow ﬁber module not only completes the development of a model that can predict streamwise variations in system performance, but also opens the doors for design and optimization studies, and reduces the number of experiments necessary to predict system-level performance. Given the complexity and diﬃculty associated with predicting NF separation, such a model is essential to: (1) evaluate the performance of the hollow ﬁber conﬁguration as applied 16 recently to NF, (2) guide the scale-up implementation of the novel low-pressure NF membranes developed, such as the LbL1.5C membrane, and (3) assess the compatibility of a speciﬁc NF membrane with desalination pretreatment applications, absent the need for expensive large-scale experiments as we show later. 2.4. Assessing Module Performance: Introducing Performance Metrics In this section, we deﬁne the performance metrics used in this study. These metrics allow us to probe and quantify the eﬀects diﬀerent design parameters, and operating conditions could have on the system. Deﬁned in Eq. 39, the membrane rejection ratio, Ri, measures the fraction of a solute that permeates through the membrane. The membrane rejection is a quantitative metric of the membrane’s selectivity, and the quality of permeate produced. The system recovery ratio, RR, as deﬁned in Eq. 40, measures the fraction of the permeate recovered from the process stream. For a given net driving pressure, a higher recovery ratio implies a higher membrane permeability and a greater system yield. The concentration polarization factor, CPi, provides a measure of the signiﬁcance of concentration polarization under given operating conditions, and is deﬁned in Eq. 41. Most importantly, the speciﬁc energy consumption, e, provides a metric to evaluate the energy cost of permeate production in the module, and is deﬁned as in Eq 42. Ri= 1ci;p ci;b(39) RR=Qp;out Qf;in(40) CPi=ci;w ci;b(41) e=Qf;inPf;in Qp;out=Pf;in RR(42) 2.5. Proposed Hollow Fiber Module Conﬁguration: Module Sizing and Design To proceed with system-level evaluation of module performance, a reference large-scale system needs to be adopted. Driven by the lack of commercial NF hollow ﬁber modules for softening and desalination pretreatment applications, we propose in this section a preliminary design inspired by earlier eﬀorts on designing UF hollow ﬁber modules [59]. In spite of expected similarities, NF membranes are expected to depart from UF standards when it comes to the operating conditions, namely the applied pressure and feed ﬂow rate. Since NF membranes have tighter pore sizes relative to UF, the applied pressures are expected to be higher in NF. Driven by diﬀerences in the process streams UF and NF are designed to handle, UF hollow ﬁbers will tend to have higher feed ﬂow rates relative to NF to prevent membrane pore clogging. The speciﬁcations of the proposed design are summarized in Table 1. 17 Table 1: Nominal speciﬁcations of modeled large-scale hollow ﬁber module. Speciﬁcation Nominal Value Membrane Area 50 m2 Module Length ( L) 1.5 m Module Diameter ( Di;module) 165 mm Packing Density ( ) 75% Fiber Inner Diameter ( di) 1 mm Fiber Outer Diameter ( do) 1.38 mm Inlet Flow Rate ( Q) 1000 L/h Similar to commercial UF hollow ﬁber modules, such as those developed by Dow [60] and Koch Membrane Systems [61], a total membrane area of 50 m2with a module length of 1.5 m, and a module diameter of 165 mm were selected. The ﬁber dimensions and operating pressures of the system have been set to match those of the LbL1.5C membranes as reported by Liu et al. [4]. These speciﬁcations in turn set the module’s packing density at 75%. Unless otherwise speciﬁed, the inlet ﬂow rate has been ﬁxed at Q= 1000L/h to minimize the energy consumption of the system while maintaining reasonable rejections as our results will demonstrate later. Given the scope of our work, the model feed was chosen to be synthetic seawater, whose ion composition and properties are provided in Table 2 for reference. Since the module performance is membrane speciﬁc, and the membrane properties are also a function of the feed chemistry, we summarize in Table 3 the rel- evant membrane parameters characterized under these feed conditions. Details concerning the membrane characterization process are provided in the Appendix for completeness. Table 2: Synthetic seawater composition and charged solute properties (pH = 6:5). Na+Ca+2Mg+2SO2 4Cl Ion Concentration (mg/kg) 11122 382 1394 2136 20300 ri(nm) 0.184 0.309 0.347 0.230 0.121 D1(m2/s109) 1.33 0.792 0.706 1.06 2.03 Table 3: LbL1.5C DSPM-DE model parameters [13]. Parameter rp(nm) xe(m)"pXd(mol/m3) 0.5 1.33 41.3 -27 18 3. Results and Discussion In this section, the developed model for predicting the performance of NF hollow ﬁber modules is vali- dated. Parametric studies are run to investigate the eﬀect of diﬀerent parameters on module performance, and oﬀer insights into module design. Speciﬁcally, the signiﬁcance of module pressure drop, including bore- side and shell-side components, is assessed. The eﬀect of feed ﬂow rate on module performance is studied, and an optimal ﬂow rate is observed. The importance of the developed model is underscored by the in- tensity of streamwise variations observed. Concentration polarization is examined under varying operating conditions. More importantly, increasing module length to increase recovery will be shown to be less eﬀective given a decaying driving force, and leading to a tradeoﬀ between CAPEX and OPEX. A comparison of marginal increases in pressure to marginal increases in recovery is employed to show that raising pressure can reduce the speciﬁc energy consumption should the growth of RRoutpace that of P, respectively. This section is concluded with a simple design methodology to guide designers interested in extending the well-established hollow ﬁber conﬁguration to nanoﬁltration. 3.1. Model Validation To assess its accuracy, the proposed hollow ﬁber model was validated with respect to experimental results reported in an earlier work [13] for the LbL1.5C hollow ﬁber membrane developed by Liu et al. [4]. System dimensions and speciﬁcations were scaled down to match those of the experimental setup, and synthetic seawater experiments (for a feed composition matching that of Table 2 and a cross-ﬂow velocity of 0.3 m/s) were adopted as the validation benchmark [13]. The membrane model parameters corresponding to those experiments match those reported in Table 3. Simulation results, which were grid independent beyond 100 cells, are summarized in Table 4. Modeling results were found to agree with experimental observations, and to fall within the experimental error margins (observed deviation was within 10% for multivalent ions) [13]. In addition to experimental error, other sources of computational error responsible for the observed deviation may be attributed to the assumptions made, such as assuming a uniform distribution of ﬁber dimensions, an equal splitting of the feed among the individual ﬁbers, or independence of the ﬁbers from one another. Negative rejection, a phenomenon well reported in literature [43, 62, 63, 64] that refers to an increase in the permeate concentration relative to the feed [13], is also well predicted by the present model. Table 4: Model validation with synthetic seawater (cross-ﬂow velocity was maintained at 0.3 m/s) [13]. Ca+2(%) Mg+2(%) SO2 4(%) RexRiRexRiRexRi 3 bar 71:1 75:790:6 85:485:4 87:4 4 bar 78:9 84:693:6 91:086:8 92:7 19 3.2. Module Pressure Drop Aspartoftheirworkinvestigatingtheeﬀectofshell-sidehydrodynamicsonmoduleperformance, Costello etal.[65]dividedtheshell-sidepressuredropintoalaminarcomponent, whichaccountsforviscousdragalong the ﬁbers, and a turbulent component, which accounts for inertial drag due to stream splitting, channeling, andmixing. Theirresultsrevealatrendfrom“inertial-drag-dominated” to“viscous-drag-dominated” pressure drop with increasing packing density. The authors attribute this shift to the non-uniformity in ﬁber packing, suggesting that uniform packing is more likely at higher packing densities. They concluded that these results fall in line with the view that ﬂow on the shell-side is three-dimensional, featuring bulk ﬂow axially and transverse ﬂow across the ﬁbers due to channeling [65]. Toassesstheimpactthepermeatesidepressuredrophasonmoduleperformanceinthiscase, aparametric study, looking into the module hydraulic losses as a function of the packing density and applied pressure, was run. The results of this study are depicted in Fig. 4, showing the permeate-side and feed-side pressure drops as functions of packing density and applied pressure, respectively. ?0.20.30.40.50.60.70.80.9" P L 00.020.040.060.080.10.120.14 3 bar 4 bar 5 bar Bore-sideShell-side (bar) Figure 4: Bore-side (feed) and shell-side (permeate) hydraulic losses as a function of module packing density and applied pressure. Except for and its corresponding RR, all nominal speciﬁcations were ﬁxed at the values reported in Table 1. As the results show, the hydraulic losses in the feed channel decrease as the packing density increases since the feed ﬂow rate, Q, is now being divided among more ﬁbers, leading to lower velocities in each ﬁber. In contrast, the hydraulic losses on the shell-side increases with increasing packing density as the recovery ratio increases and hydraulic diameter gets smaller. This increasing shell-side pressure drop is similar in behavior to earlier results reported by Costello et al. [65]. Raising pressure has a similar eﬀect by increasing the recovery ratio, which leads to higher hydraulic losses on the shell-side, and lower hydraulic losses on the feed-side. Figure 4 shows that these competing trends reach an optimum at around a packing density of 75%. 20 Comparing the hydraulic losses to the applied pressure applied reveals that the hydraulic losses are always a small fraction in this case. Figure 4 captures the fact that for reasonable packing densities ( 60:75), the permeate-side pressure can be taken to be atmospheric given the low permeate velocities. Subsequently, we adopt this assumption in the model, considerably lowering its computational expense. 3.3. Optimal Flow Rate Solute rejection is the parameter of primary interest from a softening or pretreatment perspective as it quantiﬁes the quality of water recovered from an initial feed stream. At low feed ﬂow rates, the feed concentration rises quickly, resulting in an increase in the osmotic pressure, and a corresponding decrease in permeate ﬂux. Consequently, diﬀusion of solute species from the feed to the permeate becomes more dominant, and the permeate solution increases in concentration, ultimately reducing rejections. To counteract this eﬀect, which is primarily a manifestation of concentration polarization, feed ﬂow rates need to be increased to raise mass transfer coeﬃcients, and alleviate the rise in solute concentration along the module length. As a result, diﬀusive transport of species from feed to permeate is less encouraged, keeping the permeate concentration at a minimum, and ultimately increasing rejections. To quantitatively capture this eﬀect, Fig. 5 plots the observed rejection predicted as a function of the feed ﬂow rate at an applied pressure of 5 bar. As the plot illustrates, the rejection increases with increasing ﬂow rate, until a certain threshold, Q, is reached. Qis deﬁned such that its corresponding rejection ratio is within 97% of the maximum obtainable at a ﬂow rate of 1500 L/h. Increasing the feed ﬂow rate beyond Qdoes not increase rejection signiﬁcantly, but the energy penalty incurred from the increased pumping power continues to rise. To minimize the energy requirements associated with softening or pretreatment, while maintaining high permeate quality, Qrepresents the optimal ﬂow rate at which the system should be run for a given operating pressure. Similar trends were observed for lower operating pressures. Since the LBL1.5C membrane selectively allows monovalent ions to permeate through as reported by Liu et al. [4], their rejection is not sensitive to changes in feed ﬂow rates as Fig. 5 demonstrates. 21 400 600 800 1000 1200 140000.20.40.60.81 QRiCa+2Mg+2 Na+Cl−Q∗ (L/h)SO42−Figure 5: Module average rejection as a function of feed ﬂow rate at an applied pressure of 5 bar. Except for Qand its corresponding RR, all nominal speciﬁcations were ﬁxed at the values reported in Table 1. 3.4. Streamwise Variations Another important aspect of system design is streamwise variations, which are not uncommon given the scale at which these membrane systems are implemented in practice. As the feed stream is processed along the module length, the bulk concentration of retained species in the feed channel rises. The resulting rise in feed osmotic pressure, coupled with the hydraulic losses incurred by the module, lead to decreasing rejections in the ﬂow direction. A well-designed module should in practice be “balanced”, meaning that streamwise variations from inlet to outlet are minimized, and the slopes at which the performance metrics vary per module are not steep. A well-balanced module has a relatively uniform ﬂux, decreases the likelihood of fouling, suﬀerslessfromthermodynamicirreversibilities, andensuresthemembraneareaisutilizedeﬃciently. A parametric study was run by varying the pressure and correspoding recovery ratio, while ﬁxing all other module speciﬁcations per the values reported in Table 1. The results are shown in Figs. 6 and 7. Figure 6 plots the streamwise variations in feed bulk concentration as a function of length across the system. As the ﬁgure illustrates, streamwise variations become more signiﬁcant with increasing applied pressures, or equivalently increasing recovery ratios. 22 00.511.511.21.41.61.8 zcb/cbin(a) (m)Na+Cl−Ca+2Mg+2SO42− 00.511.50.911.21.41.61.8 zcb/cbin(b) (m)Cl−Na+Ca+2Mg+2SO42−Figure 6: Streamwise variations in the feed bulk concentration: (a) applied pressure of 3 bar; (b) applied pressure of 5 bar. Except for Pand its corresponding RR, all nominal speciﬁcations were ﬁxed at the values reported in Table 1. Another point to elucidate is the order by which the concentration of the diﬀerent ion species appears to be rising in the system. The trends observed are primarily a function of the selectivity of the LbL1.5C membrane as we analyzed in a previous study [13]. Since the membrane selectively passes monovalent ions, the concentrations of these ions appear to remain constant throughout the system. According to our previous results, the LbL1.5C membrane, confronted with artiﬁcial seawater, was less selective towards SO2 4given its negative charge under these operating conditions. Given its larger Stokes radius and lower diﬀusivity, Mg+2experiences a higher rejection compared to Ca+2[13]. As expected, the concentrations of SO2 4and Mg+2rise the most according to Fig. 6. As with the feed bulk concentration, ion rejection is also a function of streamwise position. As shown in Fig. 7, the membrane retention of multivalent ions drops in the streamwise direction as a result of increasing ion diﬀusion across the membrane, and decreasing permeate ﬂux. The membrane rejection of monovalent ions stays relatively constant given the membrane’s low rejection towards them. On the other hand, rejection drops for all multivalent ions, decreasing the most for Ca+2since the membrane is more selective towards it relative to the other multivalent ions present in the solution. Similar to Fig. 6, the decrease in multivalent ion rejection becomes more signiﬁcant at higher applied pressures, indicating that a rising bulk concentration, coupled with more prominent concentration polarization, acts to encourage ion transport across the membrane. 23 00.511.500.20.40.60.81 zRi(a) (m)Na+Cl−Ca+2Mg+2SO42− 00.511.500.20.40.60.81 zRi(b) (m)Cl−Na+Ca+2Mg+2SO42−Figure 7: Streamwise variations in module rejection: (a) applied pressure of 3 bar; (b) applied pressure of 5 bar. Except for Pand its corresponding RR, all nominal speciﬁcations were ﬁxed at the values reported in Table 1. Theresultsdiscussedinthissectionunderscoretheimportanceofstreamwisevariationsalongthemodule, andtheirimpactonoverallperformance. Anotherimplicationoftheresultspresentedistheapparenttradeoﬀ a designer would face between module recovery ratio, RR, and rejection, Ri. While an improved system recovery is often desirable, a compromise might be necessary to ensure that an acceptable permeate quality is maintained. In light of this discussion, concentration polarization emerges as an important variable in module design. 3.5. Concentration Polarization Concentration polarization refers to the formation of concentration gradients at the membrane interface. From a design standpoint, concentration polarization is generally undesirable as it reduces permeate ﬂuxes, and aﬀects ion retentions. In this section, we investigate the magnitude of concentration polarization the proposed design is likely to experience during operation. The concentration polarization factor, CPi, will be the primary metric used in this analysis. The results of this analysis are represented in Fig. 8. 24 00.511.50.511.522.533.5 zCPi(a) (m)Na+Cl−Mg+2Ca+2SO42− 00.511.50.511.522.533.5 zCPi(b) (m)Cl−Na+Mg+2SO42− Ca+2Figure 8: Streamwise variations in concentration polarization factor, CPi: (a) applied pressure of 3 bar; (b) applied pressure of 5 bar. Except for Pand its corresponding RR, all nominal speciﬁcations were ﬁxed at the values reported in Table 1. Figure 8 depicts the streamwise variation in the concentration polarization factor across the system. Similar to Fig. 6, the order followed by the diﬀerent ion species in the solution is primary a manifestation of the membrane’s selectivity. The monovalent ions, which are not retained by the membrane, appear to experience very little concentration polarization. In contrast, the multivalent ions SO2 4and Mg+2, which are retained the most by the membrane, appear to experience the highest concentration polarization, followed by Ca+2, as evident from the concentration polarization factor. Another observation to point out is the growing role concentration polarization plays in the system as the applied pressure increases. Interestingly, the concentration polarization factor for all retained species decreases across the length of the module. The decrease in the concentration polarization factor is primarily a result of the driving force, which diminishes across the module length as the osmotic pressure rises and the hydraulic losses accumulate. The result is a diminishing diﬀerence between the bulk and membrane concentration. 3.6. Eﬀect of Module Length For a ﬁxed inlet feed ﬂow rate, Q, two approaches are available to the designer to reach a speciﬁc recovery ratio, increasingtheappliedpressureoravailablemembraneareainthesystem. Inthissection, weinvestigate both options by varying the applied pressure and module length, or equivalently, the number of modules connected in series. The primary limitation on increasing membrane area is the declining driving force, P, along the ﬂow direction, which renders any increase in membrane area less eﬀective at permeate recovery. In this study, the module length and applied pressure were varied for a ﬁxed feed ﬂow rate, and the resulting trends are reported in Fig. 9. 25 0 5 10 1500.20.40.60.81 zRR 3 bar 4 bar 5 bar (m)>50%∆PFigure9: Eﬀectofincreasingmodulelengthandappliedpressureon RR. Exceptforthemodulelength, allnominalspeciﬁcations were ﬁxed at the values reported in Table 1. According to Fig. 9, the recovery ratio increases with increasing applied pressure and membrane area as expected. One important takeaway from this ﬁgure, however, is that the system length and its corresponding recovery ratio are nonlinearly related. Stated diﬀerently, doubling the number of modules present in the system does not double the resulting recovery ratios due to the drop in driving force or ﬂux along the length. These results imply that while it might appear attractive to operate at lower pressures to lower the system’s operational expenditure (OPEX), the additional capital expenditure (CAPEX) and OPEX from maintenance incurred by increasing the number of modules necessary to reach a speciﬁc recovery might be unjustiﬁed. For a ﬁxed recovery ratio of 40%, for instance, our results demonstrate that raising pressure from 3 to 4 bar decreases the length required by more than a factor of two as shown in Fig. 9. Consequently, a designer should take these considerations into account when searching for an optimal operating point for the system. In the next section, we present one approach to arrive at this optimal operating point by minimizing the speciﬁc energy consumed by the system. 3.7. Module Energy Consumption Among all metrics considered in the study, energy consumption is perhaps of primary interest as it provides a handle on the energy penalty and carbon footprint associated with the NF process. In this section, the speciﬁc energy consumption, e, is chosen as the metric to evaluate the energy consumption of the proposed design. A set of parametric studies, whose results are summarized in Figs. 10 and 11, were run to investigate the dependence of the speciﬁc energy consumption on the applied pressure and feed ﬂow rate. 26 400 600 800 1000 1200 14000.20.40.6 Qe 3 bar 4 bar 5 bar (L/h)(kWh/m3) ∆PFigure 10: Speciﬁc energy consumption, e, as a function of inlet feed ﬂow rate and applied pressure. Except for Q, all nominal speciﬁcations were ﬁxed at the values reported in Table 1. Figure 10 depicts the speciﬁc energy consumption of the system as a function of the inlet feed ﬂow rate and applied pressure. The results clearly indicate an increase in the speciﬁc energy consumption with increasing feed ﬂow rate as more pumping work becomes necessary to drive the separation process. More surprisingly, however, these results suggest that the speciﬁc energy consumption decreases with increasing applied pressure for all feed ﬂow rates considered. Another signiﬁcant observation is that the incremental decrease in ediminishes with increasing applied pressures. Based on its deﬁnition, ecould decrease with increasing applied pressures only if the recovery ratio,RR, is growing at a much faster rate relative to the applied pressure. 27 33.544.555.56 00.20.40.60.81 0.250.30.350.40.450.5 33.544.555.56Rie(kWh/m3) ΔP(bar)eMg+2Ca+20!0.2!0.4!0.6!0.8!1! 0.25!0.3!0.35!0.4!0.45!0.5! 3!3.5!4!4.5!5!5.5!6!Ri e (kWh/m3) ΔP (bar) e Mg+2 Ca+2 SO SO42Figure 11: Speciﬁc energy consumption, e, and rejection, Ri, as a function of applied pressure, at the optimal feed ﬂow rate of 1000 L/h. Further analysis at the optimal feed ﬂow rate, Q, reveals that the speciﬁc energy consumption decreases with increasing applied pressure as shown in Fig. 11. In accordance with the results presented in Section 3.3, the rejections of the multivalent ions at Qappear to be constant, nonetheless, indicating that increasing the module recovery ratios by raising the pressure does not impact the permeate quality under these conditions. These results explain the trends observed earlier in Fig. 10, and indicate that the system beneﬁts from raising the pressure since the marginal increase in permeate produced surpasses the marginal increase in applied pressure. Recently, similar trends of decreasing ewith increasing RRin seawater NF have also been observed empirically by Park and coworkers [9]. These results clearly indicate that operating at lower pressures, while it may appear attractive, could be more expensive. The insights presented in this section underscore the importance of running an analysis similar to this work’s to evaluate the feasibility of a given membrane for system-level implementation in a speciﬁc application, or to help the designer locate an optimal point for operating the system. We end this section by proposing a preliminary approach to designing a large-scale NF hollow ﬁber system for a given membrane, based on experiments run on a bench-scale setup. (1) A membrane needs to be characterized according to a speciﬁc membrane transport model. (2) From there, the membrane should be modeled at the system-level in an analysis similar to the one developed in this work to locate an optimal feed ﬂow rate. (3) Once an optimal feed ﬂow rate is established, a set of parametric studies can be executed in search of an optimal operating point, at which the speciﬁc energy consumption is minimized. (4) Eventually, this proposed design can further be reﬁned by running experiments, coupled with sensitivity analysis and optimization studies with the presented model. 28 4. Conclusions In this study, a new system-level model for NF hollow ﬁber membrane modules was developed following an interdisciplinary approach featuring mass and momentum balances, coupled with a transport model based on the extended Nernst-Planck equation. The overarching objective has been to present a model that empowers researchers and designers with: (1) an understanding of the physics that govern mass and momentum transport in the module, as well as the various tradeoﬀs that exist among design parameters and operating conditions; (2) a tool that predicts NF module performance under diﬀerent operating conditions and complex feed chemistries, and enables membrane assessment for a given application, absent the need for costly large-scale experiments; and (3) criteria for module performance assessment along with a preliminary procedure for module sizing. The LbL1.5C membrane, newly developed by Liu et al. [4], was adopted in this study, and the proposed model was validated against experimental results reported previously [13]. A preliminary bore-side feed module design was proposed, and parametric studies were run to study the eﬀect of diﬀerent parameters on moduleperformance. Speciﬁcally, thechangeinhydrauliclossesontheboreandshell-sideofthemodulewith packing density and applied pressure was investigated, as well as the change in module average rejection as a function of feed ﬂow rate. Streamwise variations in feed bulk concentration and rejection were observed, and the occurrence of concentration polarization was assessed. The tradeoﬀs between module rejection, recovery ratio, module ﬂow rate, length, applied pressure, and speciﬁc energy consumption have been studied. The key ﬁndings from this study include: 1. Hydraulic losses increased on the shell-side and decreased on the bore-side, as the recovery ratio was increased, and a crossover point was observed. 2. The magnitude of the hydraulic losses on the shell-side (permeate-side) was shown to be negligible relative to the applied pressure at low to moderate packing densities. 3. Module average rejection increased with increasing inlet feed ﬂow rate as solute transport across the membrane decreased. The decrease in solute transport is a result of the diminished occurence of concentration polarization, and the more spatially-uniform feed concentrations. The trends break down beyond a certain threshold, and an optimal feed ﬂow rate was observed. 4. Streamwise variations in the feed bulk concentrations and rejection ratios were more prominent at higherappliedpressures; andtherelativemagnitudeofthesevariationswasshowntobeamanifestation of the membrane selectivity towards their corresponding species. 5. With a decaying driving force across the module length, the diﬀerence between feed bulk concentra- tion and wall concentration also decays, causing the concentration polarization factor to decrease in magnitude along the ﬂow direction. 6. Employing longer modules to increase module recovery was shown to be less eﬀective at low applied pressures as the driving force decays along the module, making the additional membrane area less eﬀective for permeate recovery. 29 7. Comparing the marginal increase in recovery ratio to the marginal increase in applied pressure shows that raising the applied pressure could lead to lower speciﬁc energy consumption if the recovery ratio grows at a greater rate relative to the applied pressure. Acknowledgements This research is supported by the National Research Foundation (NRF) Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) programme. The Center for Environmental Sensing and Modeling (CENSAM) is an interdisciplinary research group (IRG) of the Singapore MIT Al- liance for Research and Technology (SMART) centre. The authors thank: Yagnaseni Roy for advice on NF modeling; Leonardo D. Banchik for sharing his expertise on operation, design, and modeling of membrane cross-ﬂow systems; and Gregory P. Thiel for oﬀering helpful insights on osmotic pressure modeling. Omar Labban acknowledges the ﬁnancial support received from an MIT Pappalardo fellowship. 30 Appendix A. Characterization of Nanoﬁltration Membranes: The membrane characterization protocol adopted in this study follows that laid out previously in another work [13]. The procedure consists of running four experiments to characterize a membrane’s pore size, eﬀective thickness, charge density, and pore dielectric constant. The membrane pore size and eﬀective thickness are characterized by running pure water permeability (PWP) and uncharged solute experiments, while the charge density and pore dielectric constant are characterized with single salt and synthetic seawater experiments. In this section, a brief overview is provided for completeness. For more details, readers are referred to the original work [13] presented previously. A.1. Characterization of the membrane pore size and eﬀective thickness The ﬁrst step in membrane characterization involves deﬁning a pore size and an eﬀective thickness. The method adopted herein is in tandem to that of Bowen et al. [66], who in an earlier work presented an analytical solution for the uncharged solute rejection as a function of pore size and thickness for a speciﬁed ﬂux,Jv: Ri;real= 1ci;p ci;w= 1Ki;ci 1exp(Pei;w)[1iKi;c](43) Pei;w=Ki;c Ki;dJvxe Di;1(44) where Pei;wis the Péclet number comparing diﬀusive and convective transport in the pores, while the other parameters are deﬁned as in Section 2.2.2. The result accounts for concentration polarization by considering the concentration at the membrane wall ci;w. At very high ﬂuxes, this result approaches a simpler form, termed the limiting rejection [66]: Ri;lim= 1Ki;ci (45) Employing Eqs. 43 - 45 to estimate an eﬀective membrane pore size and thickness requires the mea- surement of the real rejection for uncharged solutes. Accordingly, experiments are set up such that the observed and real rejections are close by limiting concentration polarization. To that end, uncharged solute experiments in this work were operated under low concentrations (200 mg/kg), low applied pressures (1 - 4 bar), and large cross-ﬂow velocities (0.3 m/s). With these experiments, an eﬀective pore size and thickness are then deﬁned for each solute using Eqs. 43 - 44 by regression [13]. Under the assumption that uncharged solute rejection is unaﬀected by Donnan and dielectric eﬀects, Nanoﬁltran was run at exceedingly large ﬂuxes using the parameters ﬁtted to estimate each solute’s limiting rejection(whilesetting Xd= 0and"p="b= 80:4). Theresultswerethenusedtoﬁtforamembraneeﬀective pore size with the help of Eq. 45 for later use in all subsequent analyses. Once the pore size was deﬁned, the membrane eﬀective thickness was then estimated from PWP experiments using the Hagen-Poiseuille equation (given by Eq. 3). 31 A.2. Characterization of the membrane pore dielectric constant and charge density With the eﬀective pore size and thickness deﬁned, the following step in membrane characterization involves deﬁning the membrane’s dielectric constant, "p. While varying pH, single salt experiments (e.g. NaCl) are conducted to locate the membrane’s isoelectric point, where salt rejection is minimized, the membrane is assumed to be neutrally charged and the pore dielectric constant may be ﬁtted by regression. The ﬁnal step in membrane characterization is determining the membrane’s charge density, whose esti- mates are feed-speciﬁc, varying with feed pH and composition. For this reason, experiments are conducted with the feed composition of interest to the work, such as the synthetic seawater provided in Table 2. Similar to the pore dielectric constant, regression is then used in estimating the membrane charge. Once this step is completed, validation against experimental results is necessary to ensure the parameters obtained truly represent the membrane under consideration. A.3. Experimental results employed in membrane characterization In this section of the Appendix, the experimental results obtained and used for the characterization of the LbL1.5C are provided for completeness. These results have been adapted from an earlier work [13], to which the reader is referred for more details. For all experiments, the cross-ﬂow velocity was maintained at 0.3 m/s to mitigate the eﬀect of concentration polarization. For the uncharged solute and single salt experiments, the ﬂux values recorded were close to those reported for the PWP in light of the low concentrations considered. Table 5: LbL1.5C pure water permeability (PWP) measurements. Pressure applied 1 bar 2 bar 3 bar 4 bar Flux (L/m2h) 10 20 30 40 Table 6: LbL1.5C uncharged solute experimental results for a concentration of 200 mg/kg. Pressure Applied 1 bar 2 bar 3 bar 4 bar Glycerol Rejection (%) 11 15 20 29 Glucose Rejection (%) 74 77 83 83 Sucrose Rejection (%) 93 94 93 92 32 Table 7: LbL1.5C single salt experimental rejection (1000 mg/kg NaCl). NaCl rejection (%) Feed pH 1 bar 2 bar 3 bar 4 bar 5 20 31 35 36 6 17 28 31 33 7 14 20 24 25 8 14 17 18 18 9 18 22 32 34 Table 8: LbL1.5C seawater experimental rejection (composition as per Table 2). Synthetic seawater ion rejection (%) Applied Pressure Na+Ca+2Mg+2SO2 4ClFlux (LMH) 2 bar1:454.5 83.4 82.7 11 1.4 3 bar0:671.1 90.6 85.4 13 3.9 4 bar 1.17 78.9 93.6 86.8 15 6.0 33