B Bernales
email: bernales@l3m.univ-mrs.fr

P Haldenwang
email: pierre.haldenwang@univ-amu.fr

P Guichardon
email: pierrette.guichardon@centrale-marseille.fr

N Ibaseta
email: nelson.ibaseta@centrale-marseille.fr

Prandtl model for concentration polarization and osmotic counter-effects in a 2-D membrane channel

Keywords: Numerical modeling, concentration polarization, osmotic

An accurate 2-D numerical model that accounts for concentration polarization and osmotic effects is developed for the cross-flow filtration in a membrane channel. Focused on the coupling between laminar hydrodynamics and mass transfer, the numerical approach solves the solute conservation equation together with the steady Navier-Stokes equations under the Prandtl approximation, which offers a simplified framework to enforce the non-linear coupling between filtration and concentration polarization at the membrane surface. The present approach is first validated thanks to the comparison with classical exact analytical solutions for hydrodynamics and/or mass transfer, as well as with approximated analytical solutions that attempted at coupling the various phenomena. The effects of the main parameters in cross-flow reverse osmosis (RO) or nanofiltration (NF) (feed concentration, axial flow rate, operating pressure and membrane permeability) on streamlines, velocity profile, longitudinal pressure drop, local permeate flux and solute concentration profile are predicted with the present numerical model, and discussed. With the use of data reported from NF and RO experiments, the Prandtl approximation model is shown to accurately correlate both average permeate flux and local solute concentration over a wide range of operating conditions.

pressure, Starling-Darcy law, Prandtl approximation

Introduction

The first paragraph has been cancelled

In cross-flow filtration, species remain in the axial (or feed) flow, while permeate (or purified transverse flow) leaves the duct by leaking transversally through the membrane. The accumulation of rejected species at the membrane inner wall results in a concentration polarization (i.e. a large solute enhancement) which -combined with osmosis-can induce a substantial reduction in permeation. On the one hand, the basic theoretical description of the cross-flow filtration refers to the seminal contribution by [START_REF] Berman | Laminar flow in channels with porous walls[END_REF]) that treats of the channel flow driven by a uniform leakage at the wall. Other contributions on the case of the "pure solvent" flow proposed several analytical solutions that accounts for the pressure dependence on permeate flux (see [START_REF] Regirer | On the approximate theory of the flow of a viscous incompressible liquid in a tube with permeable walls[END_REF], [START_REF] Haldenwang | Laminar flow in a two-dimensional plane channel with local pressure-dependent crossflow[END_REF], [START_REF] Tilton | Incorporating darcy's law for pure solvent flow through porous tubes: Asymptotic solution and numerical simulations[END_REF], [START_REF] Bernales | Laminar flow analysis in a pipe with locally pressure-dependent leakage through the wall[END_REF]). On the other hand, an exact analytical solution for the solute transfer in Berman flow has been derived in [START_REF] Haldenwang | Exact solute polarization profile combined with osmotic effects in berman flow for membrane cross-flow filtration[END_REF] and accounts for concentration polarization and the subsequent hindrance to permeation in the limit of certain RO/NF configurations.

Permeate flux actually results from the coupling of those transport phenomena (hydrodynamics and mass transfer) upstream and through the membrane.

The present contribution can essentially be seen as an improvement in the numerical efficiency for studying the laminar flows involved in filtration. It is however known that the applicability domain of such an approach can be extended by introducing the concept of turbulent viscosity and turbulent diffusivity (as we shall do in the present Subsection 3.3).

This issue has also been recently identified as a challenge for nanofiltration and reverse osmosis processes [ [START_REF] Van Der Bruggen | Drawbacks of applying nanofiltration and how to avoid them: A review[END_REF] and [START_REF] Malaeb | Reverse osmosis technology for water treatment: State of the art review[END_REF]]. The literature is now rich with studies on RO/NF modeling. In the present approach, we hence pay a special attention to develop an efficient predictive model that allows us to accurately describe the full coupling between flow, concentration polarization and hindrance to permeation. As discussed below, efficiency and accuracy come from the fact that the solver is fast and has no limitations in terms of numerical degrees of freedom.

Whereas the conservation laws for diluted solutions of salt and water are well-admitted, solvent and solute transport through a membrane remains a quite controversial issue. A large number of investigators have proposed numerous local models derived from different mechanisms. They are based either on the principles of irreversible thermodynamics models (Kedem-Katchalsky model, Spiegler-Kedem model) or on various homogeneous membrane models (Porous model, Solution-Diffusion model). The oldest of these models are analyzed in the contribution by Soltanieh [START_REF] Soltanieh | Review of reverse osmosis membranes and transport models[END_REF]]. More recent models describing the different local transport phenomena within the membranes with an increasing complexity can be found in [START_REF] Wijmans | The solution-diffusion model: a review[END_REF], [START_REF] Gauwbergen | Modelling revrese osmosis by irreversible thermodynamics[END_REF], [START_REF] Weissbrodt | Separation of aqueous organic multi-compount solutions by reverse osmosis -development of a mass transfer model[END_REF], [START_REF] Moresi | Modelling of ammonium fumarate recovery from model solutions by nanofiltration and reverse osmosis[END_REF], [START_REF] Kahdim | Modeling of reverse osmosis systems[END_REF], [START_REF] Mehdizadeh | Modeling of mass transport of aqueous solutions of multi-solute organics through reverse osmosis membranes in case of solute-membrane affinity part1. model development and simulation[END_REF], [START_REF] Mane | Modeling boron rejection in pilot-and full-scale reverse osmosis desalination processes[END_REF], [START_REF] Malaeb | Reverse osmosis technology for water treatment: State of the art review[END_REF]. Below, we shall nevertheless use a minimal model, since it is essential to investigate the associated limitations.

An evident improvement in modeling consists in coupling a local membrane model with the upstream composition of the solution. Several investigations [START_REF] Alvarez | Permeate flux prediction in apple juice concentration by reverse osmosis[END_REF], [START_REF] Jamal | Mathematical modeling of reverse osmosis systems[END_REF], [START_REF] Prabhavathy | Treatment of fatliquoring effluent from a tannery using membrane separation process: Experimental and modeling[END_REF], [START_REF] Hung | Mass-transfer modeling of reverseosmosis performance on 0.5-2% salty water[END_REF], [START_REF] Choi | Modeling of full-scale reverse osmosis desalination system: Influence of operational parameters[END_REF], [START_REF] Qiu | Concentration polarization model of spiralwound membrane modules with application to batch-mode ro desalination of brackish water[END_REF]] have predicted permeate flux by using the combined solution-diffusion/film model also known as the Kimura-Sourirajan model. The integrated or differential film theory equations (convection due to the pressure difference and back diffusion owing to the concentration gradient) have also been used by Urama [START_REF] Urama | Mechanistic interpretation of solute permeation through a fully aromatic polyamide revrese osmosis membrane[END_REF]] and Ahmad [START_REF] Ahmad | Mathematical modeling of multiple solutes system for reverse osmosi process in pal oil mill effluent (pome) treatment[END_REF]] while they applied the Spiegler-Kedem models in the membrane. It is worth noting that all these studies are based only on mass transfer considerations forgetting the influence of hydrodynamics. Such an approach that involves only averaged data all along the membrane is attractive in view of its simplicity and surely convenient for the description of a complex multistage membrane process in transient condition. However, just using averaged conditions assimilated to feed parameters along the membrane length becomes seemingly unrealistic when an industrial membrane module where a spatial variation of pressure and solute concentration appears in the long cross-flow membrane channel is involved.

We hence turn towards 1-D models that consider the local longitudinal variations along membrane duct. The simplest approach we can mention conceives the feed channel as a black box which can be represented by a series of perfect mixing cells with exchange. No information is needed on hydrodynamics and mechanisms of transport. This is the purpose of [START_REF] Roth | Modelling of stimulus response experiments in the feed channel of spiral-wound reverse osmosis membranes[END_REF] who used the residence time distribution (RTD) method by analyzing in some spiral-wound membranes the response to a stimulus injection of tracer. The main drawback of this very simplified method is to ignore the local solute concentration in boundary layer which has a great influence on permeate flux via the phenomenon of concentration polarization.

In other 1-D studies, the material balance is solved numerically without taking the velocity field into account, the hydrodynamic effects being limited to pressure drop along the membrane, which is generally described by the Hagen-Poiseuille or the Ergun equations [START_REF] Sekino | Precise analytical model of hollow fiber reverse osmosis modules[END_REF], [START_REF] Sekino | Study of an analytical model for hollow fiber reverse osmosis module systems[END_REF], [START_REF] Chatterjee | Modeling of a radial flow hollow fiber module and estimation of model parameters using numerical techniques[END_REF], [START_REF] Senthilmurugan | Modeling of a spiralwound module and estimation of model parameters using numerical techniques[END_REF]]. Material balances have also been considered numerically with the assumption of a negligible pressure drop by [START_REF] Malek | A lumped transport parameter approach in predicting b10 ro permeator performance[END_REF]] who introduced a simple model based on a lumped transport parameter approach combined with the solution-diffusion model. A modified solution-diffusion model [START_REF] Sagne | Modeling permeation of volatile organic molecules through reverse osmosis spiral-wound membranes[END_REF]] was developed to take the sorption pattern into account.

Several approaches investigated the cell model based on single material balances and consisted in dividing the membrane into several parts that were considered as a well-mixed reactors. [START_REF] Voros | Salt ans water permeability in reverse osmosis membranes[END_REF]] neglected pressure drop in comparison with the trans-membrane applied pressure, while others works treated of the influence of pressure drop [START_REF] Costa | Modelling of modules and systems in reverse osmosis. part i. theoretical system design model development[END_REF], [START_REF] Ali | Modeling the transient behavior of an experimental reverse osmosis tubular membrane[END_REF], [START_REF] Fujioka | Modelling the rejection of n-nitrosamines by a spiral-wound reverse osmosis system: Mathematical model development and validation[END_REF]]. MATHCAD8 package was also used to solve the mass transport equations across the membrane and the boundary layer [START_REF] Albastaki | Predicting the performance of ro membranes[END_REF]]. Let us also mention that an analytical model treating the membrane as a heterogeneous system [START_REF] Song | Performance of a long crossflow reverse osmosis membrane channel[END_REF]] was obtained with the assumption of constant driving pressure.

To conclude on the large amount of 1-D models, we stress on the fact that these approaches at least postulate one of the following strong simplifications:

-simplified pressure drop derived from the particular case of flow in an impermeable duct, -permeate flux and concentration polarization are estimated from mass transfer coefficients derived from dimensionless correlations on Sherwood number;

-hydrodynamics is not investigated, only solute mass transfer is considered neglecting the influence of the cross-flow on solution composition.

It is however accepted that permeate flux has a great influence on crossflow, particularly on pressure drop in the situation of highly permeable membrane [START_REF] Mellis | Fluid dynamics in a tubular membrane: Theory and experiment[END_REF], [START_REF] Haldenwang | Pressure runaway in a permeable channel with pressure-dependent leakage[END_REF]]. Concerning the use of mass transfer coefficients, the difficulty lies in the proper choice of the most suitable correlation among the large number of correlations found in literature which are often eclectic. As concentration polarization results from the establishment of a solute boundary layer, it is essential to couple hydrodynamics and mass transport.

We now devote what follows to more accurate models that deal with the 2-D variations of pressure, flow velocities and solute concentration in the feed channel.

We first consider the analytical approaches, which are approximate solutions for the time being, even though some exact solutions have been derived in certain limit cases. The first concern of the approximate analytical approaches considers the various manners of incorporating a satisfactory velocity field into the mass transfer equation. Then, the film theory is often invoked to describe the coupling between permeate flux and concentration polarization. In this way, an analytical expression of the permeate flux can be established. This has firstly been done for the case of a solid suspension [START_REF] Song | Theory of concentration polarization in crossflow filtration[END_REF]]. Song and Elimelech used the excess concentration subtracted by feed concentration to balance the convection and diffusion of solute within the polarization concentration layer. The validation analysis of such an approach can be found in [START_REF] Kim | Modeling concentration polarization in reverse osmosis processes[END_REF]]. Later, this model was extended and applied to salt solutions specific to RO membranes [START_REF] Song | Concentration polarization in cross-flow reverse osmosis[END_REF]]. The specific case of a long-narrow channel in which concentration polarization could develop all along the channel was also investigated [START_REF] Song | Concentration polarization in a narrow reverse osmosis membrane channel[END_REF]]. To take account of the variable cross-flow velocity, a total salt balance model was proposed for different considered kinds of shear flow [START_REF] Song | A total salt balance model for concentration polarization in crossflow reverse osmosis channel with shear flow[END_REF]]. Let us mention that Sundaramoorthy [START_REF] Sundaramoorthy | An analytical model for spiral wound reverse osmosis membrane modules: Part i -model development and parameter estimation[END_REF]] also proposed an analytical model that provides us with explicit expressions for spatial variations of pressure, fluid velocity and solute concentration on the feed channel side of a spiral wound RO/NF membrane module.

If the permeation is supposed to remain more or less uniform all along the channel, the flow field can be described by the Berman exact solution [START_REF] Berman | Laminar flow in channels with porous walls[END_REF]]. This expression is then used for solving the mass transport equation in the diffusion layer [START_REF] Agashichev | Modeling the influence of temperature on gelenhanced concentration polarization in reverse osmosis[END_REF]]. [START_REF] Kim | Permeat flux inflection due to concentration polarization in crossflow membrane filtration: A novel analytic approach[END_REF]] solved the convection-diffusion equation with a Berman flow approximated as a linear shear flow in the close proximity of the membrane surface. In Annex, we briefly recall the exact analytical solution to the mass transfer in a solution carried by a Berman flow, as obtained by [START_REF] Haldenwang | Exact solute polarization profile combined with osmotic effects in berman flow for membrane cross-flow filtration[END_REF]]. Note that the latter solution precisely describes the concentration polarization in the HP-LR limit (HP-LR for High Pressure -Low Recovery).

In addition to these analytical contributions, numerous numerical studies have to be mentioned. Bhattacharyya [START_REF] Bhattacharyya | Prediction of concentration polarization and flux behavior in reverse osmosis by numerical analysis[END_REF]] developed a numerical approach that solves the diffusion-convection equation to compute the concentration profiles throughout a reverse osmosis membrane. Ma and Song [START_REF] Ma | A 2-d streamline upwind petrov/galerkin finite element model for concentration polarization in spiral wound reverse osmosis modules[END_REF]] solved the convection-diffusion equation coupled with Navier-Stokes equations in the feed channel. This work is based on rigorous mass and momentum balances in both the radial and axial membrane dimensions. We now reach an important aspect in the present modeling: it concerns the calculation of the hydrodynamics and mass transfer coupling with high accuracy thanks to computational fluids dynamics (CFD) softwares. Wiley and Fletcher developed a CFD model of concentration polarization and fluid flow in membrane process [see [START_REF] Wiley | Techniques for computational fluid dynamics modelling of flow in membrane channels[END_REF]]. Later, they extended the model to the specific filtration processes [START_REF] Fletcher | A computational fluids dynamics study of buoyancy effects in reverse osmosis[END_REF], [START_REF] Alexiadis | CFD modelling of reverse osmosis membrane flow and validation with experimental results[END_REF]]. A great advantage of the CFD methods lies in their ability to treat of complex system geometries. Certain CFD software are indeed convenient to simulate momentum and mass transport in membrane filtration system with spacer-filled channels [START_REF] Subramani | Pressure, flow and concentration profiles in open and spacer-filled membrane channels[END_REF]] in order to modify the cross-flow. The fully coupled governing equations for fluid dynamics and mass transfer were investigated [START_REF] Lyster | Numerical study of concentration polarization in a rectangular reverse osmosis membrane channel: Permeate flux variation and hydrodynamic end effects[END_REF]], and more recently in [START_REF] Salcedo-Diaz | Visualization and modeling of the polarization layer in crossflow reverse osmosis in a slit-type channel[END_REF]], where the CFD sofware Comsol Multiphysics has been used. A 3-D numerical solution of the coupled fluid dynamics and solute transfer equation [START_REF] Lyster | Coupled 3-d hydrodynamics and mass transfer analysis of mineral scaling-induced flux decline in a laboratory plate-and-frame reverse osmosis membrane module[END_REF]] was obtained with the use of Ansys CFX solver.

It is worthy of note that the numerical approach in cross-flow filtration is faced with an unusual boundary condition that nonlinearly relates permeation flux and concentration at the membrane surface. The general CFD softwares treat this difficulty in an iterative way that is founded on standard methods for solving large non-linear systems; in practice, this limits the number of nodes to get a fast convergence. For instance, in [START_REF] Salcedo-Diaz | Visualization and modeling of the polarization layer in crossflow reverse osmosis in a slit-type channel[END_REF]] the number of nodes is about 3.10 4 . The aim of the present numerical approach is to develop a tailor-made modeling, which does not have such a limitation (our standard runs will be able to involve 2.10 7 nodes and will take one ten of seconds on a standard Intel I5 processor).

Actually, an elegant manner to save computational efforts is to con-sider the laminar Navier-Stokes equations in the limit of the Prandtl approximation. This approach -fully justified in most channel laminar filtration systems-involves a "time marching like" solver along the channel, which makes easier the treatment of the non-linear boundary condition. The Prandtl approximation applies on all conservation laws for a chemical solution flowing in a cross-flow filtration channel. Therefore, the present approach aims at accounting for concentration polarization and osmotic counter-pressure, expecting to predict permeate flux and concentration profile, at any point all along the channel, with a computational cost low enough to develop an exhaustive parametric study. Furthermore, we focus our interest on cross-flow filtration at steady state. The fluid is incompressible and the channel is supposed to be narrow. As a result, this paper takes place within a continuous effort towards predictive filtration models, the efficiency of which increases in terms of precision and in capacities of analysis.

It is organized as follows: In section 2, we present the governing equations and the main hypotheses. We briefly recall why the Prandtl approximation of the conservation laws is justified. The numerical method that solves the Prandtl differential system is described. In section 3, we first compare our numerical results with the Berman exact solution and the associated exact solution for mass transfer [START_REF] Haldenwang | Exact solute polarization profile combined with osmotic effects in berman flow for membrane cross-flow filtration[END_REF]. Second, other approximated models (Elimelech's model and Total Salt balance model) are also compared with the present numerical predictions for the case of a pure solvent and of a single solute solution. Finally, we present the comparison between our model predictions and several experimental data from reverse osmosis experiments of the literature. In Appendix, we briefly recall the spirit of the different analytical results that we need for conducting the validation step [START_REF] Haldenwang | Exact solute polarization profile combined with osmotic effects in berman flow for membrane cross-flow filtration[END_REF], [START_REF] Song | Theory of concentration polarization in crossflow filtration[END_REF], [START_REF] Song | A total salt balance model for concentration polarization in crossflow reverse osmosis channel with shear flow[END_REF]].

Numerical solution of Prandtl system for filtration

Main hypotheses and conservation laws at steady state

The present 2-D numerical development is focused on modelling the concentration polarization and the osmotic (counter-)pressure that occur in the steady cross-flow filtration with one solute diluted solution. The experiments on reverse osmosis show that the rejection of solute is nearly total in most situations. To simplify the approach -and reduce the number of independent parameters-we suppose that the solute rejection is total.

The Prandtl approximation corresponds to a simplification that produces an important gain in numerical efficiency. The price to pay is the following: when we shall decide to cancel the streamwise diffusive terms, we shall renounce to study the details of any recirculating phenomena, as what happens between two spacers.

Two semi-permeable parallel walls compose the 2-D channel as schematically shown in Figure 1. The channel is of length L and spacing 2d. This defines the computational domain as {-d < x < d} × {0 < z < L}. The fluid is supposed to be Newtonian, incompressible and its physical properties, such as dynamical viscosity µ 0 and density ρ 0 , are supposed uniform. Inlet conditions (at z = 0) are the given axial mean velocity W in , the fixed feed concentration C in and the fixed inlet pressure P in . To simplify our development again, the effect of partial fouling or cake formation on permeate flux is supposed to be inexistent or negligible in comparison with the hindrance due to the osmotic (counter-)pressure. In other words, I 0 , the membrane resistance is considered as constant and uniform all along the channel. This allows us to define U in , the permeation velocity in pure solvent cross-flow, as

Solute concentration profile

C in P in W in d z x -U W (z) U W (z) = P W (z)-P osm W (z) I 0
U in = P in I 0 (1)
Accordingly with the Darcy law for porous media and the van't Hoff law for osmotic pressure, and considering a total rejection of the membrane, the Spiegler-Kedem model leads us to the so-called Darcy-Starling law which relates U W (z), the local permeation velocity at the wall, to pressure P W and osmotic pressure P osm W as:

U W (z) ≡ P W (z) -P osm W (z) I 0 = P W (z) -iRT C W (z) I 0 ( 2
)
where i is the number of dissolved ions per molecule of salt, T , the temperature of the solution and R, the ideal gas constant. We consider a twodimensional newtonian flow in cartesian coordinates (x, z) with the velocity vector {U, W }. The classical set of conservation laws reads:

∂U ∂ x + ∂W ∂ z = 0 (3) U ∂U ∂ x + W ∂U ∂ z = - 1 ρ 0 ∂P ∂ x + µ 0 ρ 0 ∂ 2 U ∂ x2 + ∂ 2 U ∂ z2 (4) U ∂W ∂ x + W ∂W ∂ z = - 1 ρ 0 ∂P ∂ z + µ 0 ρ 0 ∂ 2 W ∂ x2 + ∂ 2 W ∂ z2 (5) U ∂ C ∂ x + W ∂ C ∂ z = D 0 ∂ 2 C ∂ x2 + ∂ 2 C ∂ z2 ( 6
)
where C is the field of solute concentration, the molecular diffusivity of which is D 0 . This set of partial differential equations has to be solved with the following boundary conditions: at the porous walls (i.e x = ±d, ∀z)

W = 0 (7) U C -D 0 ∂ C ∂ x = 0 (8) U = P -iRT C I 0 (9)
and with appropriate inlet/outlet conditions at z = 0 and z = L. Note the nonlinearity present in equation ( 8).

Non-dimensioning

At this point, we have several order of magnitude already defined: d, U in , W in , P in , C in determine the order of magnitude for the variations of x, U , W , P and C. To set an order of magnitude for the axial coordinate, we add L de , the so-called dead-end length (or exhaustion length for clear water) [START_REF] Haldenwang | Laminar flow in a two-dimensional plane channel with local pressure-dependent crossflow[END_REF]),

L de = W in d U in = I 0 W in d P in (10)
Let us now define the following dimensionless unknowns and variables, used throughout the article:

u = U U in , w = W W in , p = P P in , x = x d , z = z L de , C = C C in (11)
This process leads us to introduce the six dimensionless numbers that characterize the problem:

R in = ρ 0 U in d µ 0 , λ = L L de , α = µ 0 I 0 W 2 in P 2 in d 1 2 , β = µ 0 I 0 d , P e in = P in d D 0 I 0 and N osm = iRT C in P in (12)
R in is the "pure solvent" transverse Reynolds number, λ is the dimensionless length of the channel reduced by L de , the exhaustion length (or dead-end length). In practice, experimentalists choose λ < 1. α is the square root of the ratio of Hagen-Poiseuille pressure drop throughout the exhaustion length to trans-membrane pressure and β is the dimensionless permeability of the membrane, the typical values of which are much less than unity (for a more detailed discussion on those numbers, refer to [START_REF] Haldenwang | Laminar flow in a two-dimensional plane channel with local pressure-dependent crossflow[END_REF]). P e in is the "pure solvent" transverse Péclet number. N osm , the ratio of the osmotic pressure in absence of polarization to the operating pressure P in , will be denoted the osmotic number. Note that N osm must be less than 1 in pressure-driven filtration.

We are now able to rewrite the whole system, in the domain ]-1, 1[×]0, λ[, as the following non-dimensional form:

0 = ∂u ∂x + ∂w ∂z (13) - ∂p ∂x = β R in u ∂u ∂x + w ∂u ∂z - ∂ 2 u ∂x 2 + β α 2 ∂ 2 u ∂z 2 (14) - ∂p ∂z = α 2 R in u ∂w ∂x + w ∂w ∂z - ∂ 2 w ∂x 2 + β α 2 ∂ 2 w ∂z 2 (15) 0 = P e in u ∂C ∂x + w ∂C ∂z - ∂ 2 C ∂x 2 - β α 2 ∂ 2 C ∂z 2 (16)
with the boundary conditions at the wall (x = ±1):

w(±1, z) = 0, u(±1, z) = p -N osm C W , P e in uC = ∂C ∂x (17)
with the appropriate inlet/outlet conditions at z = 0 and z = λ.

Prandtl system for cross-flow filtration

Note that the ratio β/α 2 reduces to U 2 in /W 2 in , which is by nature of RO/NF cross-flow filtrations a very small quantity (as small as 10 -6 ).

Incompressibility condition (13) confirms that the non-dimensional quantities u and w have the same order of magnitude for their different variations of the RHS of equations ( 14) and ( 15). Therefore, transverse variation of pressure is β/α 2 times as small as its axial variation. It can easily be deduced that pressure is constant within a channel section (i.e p(x, z) ≡ p(z)).

The same rationale applies when comparing transverse diffusion with axial diffusion. It is clear that both magnitudes are in U 2 in /W 2 in ratio. This incites us to cancel the axial diffusion terms in the differential equations.

Consequently, the following set of equations are validated for a large class of filtration processes. This is the so-called Prandtl approximation of the conservation laws. Note that the mathematical nature of these equations has changed, since we switch from a set of elliptic equations to a set of parabolic equations, in which the variable "z" plays the same role as the time in the heat equation.

∂u ∂x + ∂w ∂z = 0 (18) ∂p ∂x = 0 (19) R in w ∂w ∂z - ∂ 2 w ∂x 2 = -R in u ∂w ∂x - 1 α 2 dp dz (20) P e in w ∂C ∂z - ∂ 2 C ∂x 2 = -P e in u ∂C ∂x (21)
In the left-hand-side of equations (20-21), we have gathered the terms of these equations that are mathematically identical to that of heat equation. This Prandtl system for cross-flow filtration calls for the following comments.

Remark 1: If the above inlet data at z = 0 are symmetrical, the laminar solution of the system will develop symmetrically with respect to x = 0. Therefore, assuming symmetrical boundary conditions at x = 0 allows us to gain half of the computational effort, the computational domain reducing to {0 ≤ x ≤ 1} × {0 ≤ z ≤ λ} with the following boundary conditions on the axis (x = 0):

u(0, z) = 0, ∂w ∂x (0, z) = 0, ∂C ∂x (0, z) = 0 (22)
This assumption will hold in all what follows.

Remark 2 : Equation ( 19) expresses the fact that pressure only depends on z, the stream-wise position. Equations ( 20) and ( 21) are clearly of parabolic type, the pressure being -in each section-the Lagrange multiplier that allows us to satisfy the following constraint

d dz 1 0 wdx = -u(1, z) = -[p(z) -N osm C W (z)] ( 23
)
which is obtained owing to transverse integration of equation ( 18). The forthcoming numerical method exploits this mathematical property. Furthermore, the parabolic nature of this system only requires boundary conditions at inlet. Hence, the term "appropriate boundary conditions" previously mentioned after equations ( 17) becomes now defined and simply reduce to entrance conditions at z = 0: only, suitable entrance x-profiles on w and C are required.

Remark 3: Conservation laws are non-linear by nature. Here, their nonlinearity is increased by the boundary condition: P e in uC -∂C/∂x = 0 at x = 1, which non-linearly combines permeation velocity and solute concentration. This coupling is tremendously important since it governs polarization and permeation. Therefore, this nonlinearity needs to be numerically enforced iteratively, as described below.

A new numerical approach

Technically, the solution method for solving the steady conservation laws: equations from eq.( 18) to eq.( 21), coupled with the boundary conditions ( 22) and ( 17) is performed by using finite difference methods (FDM) of order two. The computational domain is discretized using a regular mesh in both transverse and axial direction. A regular mesh is not optimal transversally, but a large discretization is permitted, since the numerical cost is low. As a matter of fact, the conservation laws being conceived in the context of Prandtl approximation, the system is parabolic, the axial coordinate playing the role of the time. Hence, for a given transverse section z, we solve the coupling -which is non-linear-by an iterative process. One iteration is as follows: we first solve the concentration field, then the axial velocity field together with pressure, and finally the transverse velocity field. The solver iterates until a certain convergence criterion is satisfied in what concerns the simultaneous satisfaction of all boundary conditions or couplings, in particular, when the consistency between local permeation velocity U W , local pressure p and local wall concentration C W is achieved. At this convergence stage, we are in condition to proceed with the computation of the various fields in the next transverse section until the whole length of the channel is reached.

As the system of equations [from eq.( 18) to eq.( 21)] is of parabolic type, we now resort to the terminology generally used for the discretization of heat equation with finite differencing. Let us define w

(n) j [resp. C (n) j ]
, the value of unknown w [resp. C] at point {x = j∆x, z = n∆z} for 0 ≤ j ≤ J and 0 ≤ n ≤ N , where ∆x = 1/J and ∆z = λ/N are the mesh sizes in both directions. We suppose that all transverse profiles (velocity field {u

(n) j , w (n) j }, scalar concentration field C (n) j
and local pressure p (n) ) are known up to the section z = n∆z. The purpose of the numerical scheme is to compute the all unknown fields ( {u

(n+1) j , w (n+1) j }, C (n+1) j
and p (n+1) in the next section z = (n + 1)∆z.

In channel section z = (n+1)∆z , the permeation velocity u (n+1) J is a keystone for the numerical method. Since the permeation velocity is involved in the non-linear boundary condition, it assessment will be obtained as the convergence of an iterative process where U (k) , k = 0, 1, 2, • • • , represents a series of local permeation estimates. At convergence of the iterative process, we shall set u (n+1) J = lim U (k) . The final numerical system to be solved incorporates a discretization of the differential operators. More precisely, for the transverse coordinate, the following centered finite difference operators of order two are chosen:

∂w ∂x (n+1) j ≈ w (n+1) j+1 -w (n+1) j-1 2∆x (24) ∂ 2 w ∂x 2 (n+1) j ≈ w (n+1) j+1 -2w (n+1) j + w (n+1) j-1 (∆x) 2 , ( 25
)
As for the axial differential operators, a backward finite difference scheme of order 2 is used:

∂w ∂z (n+1) j ≈ 3w (n+1) j -4w (n) j + w (n-1) j 2∆z ( 26
)
To perform the axial extrapolations, we selected the Adams-Bashforth second order scheme: ŵ(n+1)

j ≈ 2w (n) j -w (n-1) j (27)
The entire procedure starts with the setting of the entrance data. Operating pressure is p (0) = 1. The velocity fields u

(0) 0,•••,j,•••,J , w (0) 0,•••,j,•••,J is either of Berman type, or of Poiseuille type. The entry concentration field C (0) 0,•••,j,•••,J
is supposed transversally uniform and set to C (0) j = 1. Because we use the extrapolation scheme 27 to determine the guessed next fields, for the first (unknown) transverse section (n = 1) we would also need all the fields in n = -1, which can be those of the section (n = 0). Hence, we set {u

(-1) j = u (0) j , w (-1) j = w (0) j and C (-1) j = C (0) j }.
Now we are ready to enter into the main axial loop. The procedure described bellow is the same for all transverse section z = (n + 1)∆z from n = 0 to N -1. For a general node n + 1 the above extrapolation provides us with the initial guesses { ŵ(n+1)

j , û(n+1) j , Ĉ (n+1)
j } for 0 ≤ j ≤ J. At iteration k = 0, we use this guess for the transverse velocity by setting

U k=0 = û(n+1) j
. ¿From the discretized form of equation ( 21), and boundary nodes (symmetry condition at j = 0, and the Robin boundary condition at j = J), we then obtain a tri-diagonal system to inverse for computing the solute concentration profile C (n+1) j . We now proceed with the computation of the axial velocity field w (n+1) . In the same manner as that used for the concentration field, equation ( 20) is developed by its respective discretization. It is worthy noting that the no-slip condition w (n+1) J = 0 (i.e. at node j = J ) reduce the unknown field w (n+1) to J unknowns. However the whole system has always J + 1 unknowns due to additional unknown p (n+1) . Therefore, a complementary constraint is required. The latter unknown depends on boundary condition u = p -N osm C that holds at the membrane (x = 1) and after integration of the incompressibility constraint 18 we obtain another relationship to relate all unknowns (w (n+1) j to p (n+1) ) and complement the tri-diagonal system. From the profile w (n+1) , the incompressibility constraint allows us to compute the transverse velocity profile u (n+1) and the permeate velocity

U (K) At this stage, all fields C (n+1) j , w (n+1) j , u (n+1) j
and p (n+1) are known. The new permeation velocity can be computed thanks to net trans-membrane pressure, as

U k+1 = p (n+1) -N osm C (n+1) J . ( 28
)
We then check the convergence of the procedure by comparing the new filtration with that of the previous step in the loop. More precisely, we

if n < N n = n + 1 Set n = 0
Assign values to u (0) , w (0) , C (0) , p (0) Set values to u (-1) , w (-1) , C (-1) , p (-1) Compute values of û(n+1) , ŵ(n+1) , Ĉ(n+1) , p(n+1)

Set k = 0 Resolve system C (n+1) Resolve system w (n+1) Compute values of u (n+1)
Store values of u (n+1) , w (n+1) , C (n+1) , p (n+1) Set values of u (n-1) , w (n-1) , C (n-1) , p (n-1) Set values of u (n) , w (n)

, C (n) , p (n) End If U (k+1) -U (k) U (k) ≤ Econv k = k + 1 U (k) = U (k+1)
U (k) -U (k+1) U ( k) ≤ E conv (29)
where E conv is a certain criterion of convergence. If convergence is not reached, a new iteration is performed for the same section z = (n + 1)∆z with the new guess of permeation

U k = U k + ω(U k - U k+1
), where ω is a relaxation factor used to speed up the convergence, until an acceptable permeation velocity U K satisfies the convergence criterion. After each iterative loop, all fields

C (n+1) j , w (n+1) j , u (n+1) j
and p (n+1) are stored. If convergence is attained, the computation of the next axial section is prepared by setting all the memories affected to section (n -1) with the corresponding values just obtained for section (n) and fields of section (n) with the values of section (n + 1). Now we are in condition to restart the previous iterative process for the next section z = (n + 1)∆z, following the same calculation sequence until z = λ = N ∆z (i.e. n = N ). All calculation steps are schematically described in figure 2.

A classical study of dependency on mesh size has been performed [START_REF] Bernales | Modélisation de l'hydrodynamique et des transferts dans les procédés de filtration membranaire[END_REF] in both axial and transverse directions; this numerical check confirms the order two expected from the finite difference approximations, that we have used.

Validation and results

Even though we have derived a numerical approach rather light, we need to validate the numerical technique, as well as the spirit of the Prandtl approximation. The best manner is to resort to exact analytical solutions, even when their efficiency with respect to filtration is limited. Once this step is completed, we can tackle the numerical predictions of well-documented experiments.

Validation with exact analytical solutions

First of all, the numerical model is validated by comparison with the Berman flow. This validation implies the simulation of a uniform permeation of pure solvent, namely a pure water flow nearly isobaric. Accordingly with the validity domain of the Berman solution, this is achieved when parameter α becomes as small as possible, since the pressure drop is proportional to α 2 . In Table 1, the numerical data obtained with three studied cases where α rapidly declines are gathered. E p is the maximum of the relative error between the numerical result (average permeate flux) and the solution of reference (here the Berman flow). It is clear that the departure diminishes as the square of α 2 , indicating that the error only comes from the pressure drop, which induces a non-uniform filtration in the numerical simulation.

Test Case α R in number of meshes in the axial direction (N) E p 1 10 -2 0.1 10 4 10 -4 2 10 -3 0.1 10 4 10 -6 3 10 -4 0.1 10 4 10 -8

Table 1: Error with respect to the simple case of uniform permeation (Berman exact solution)

We retrieve the fact that pressure becomes uniform as α 2 diminishes, as well as permeation (according to Darcy's law).

For the sake of illustration, Figure 3 shows the axial and transverse velocity profiles. Both numerical and analytical solutions (Berman flow) are plotted for three axial sections, with no visible departures. We observe that the axial velocity component is close to a parabola (but slightly flatter at the maximum). As expected, the analytical solution well-predicts the velocity (α = 0.01). The reduced pressure p is quietly constant in that case (figure 4d), leading to a constant local permeate flux, which does not depend on the axial position (figure 4b). Hence, the axial flow q decreases linearly as seen in figure 4d. There is consequently a very good agreement with the Berman's theory in this case, as expected.

For the second case, an important pressure drop was tested (α = 0.5). The reduced pressure p decrease all along the channel (figure 5d), leading to a decreasing local permeate flux (figure 5b). Then, the axial flowrate decreases more slowly than expected (figure 5d), and the exhaustion length, z = 1, is reached without completely exhausting the axial flow (figures 5c, 5d). As seen in figure 5d, the reduced flowrate q decreases faster than the reduced pressure p; hence, if the channel would be longer, the axial flow exhaustion (AFE) would occur rather than a cross-flow reversal (CFR). This is consistent with a previous work [START_REF] Haldenwang | Laminar flow in a two-dimensional plane channel with local pressure-dependent crossflow[END_REF] predicting this behaviour for these values of α and R in (R in ≈ 0, α < 1/ √ 3). If the relative pressure drop is further increased (α = 0.75), the reduced pressure p will decrease faster than the flowrate q, as seen in figure 6d. At z = 0.8, the reduced pressure becomes negative (that is, the pressure in the retentate side becomes lower than the pressure in the permeate side), leading to a cross-flow reversal (CFR): water flows from the permeate to the retentate, as expected for R in ≈ 0, α > 1/ √ 3 [START_REF] Haldenwang | Laminar flow in a two-dimensional plane channel with local pressure-dependent crossflow[END_REF]). Hence, the reduced axial flowrate q increases again for z > 0.8. Figure 6c confirms us this behavior when streamlines at z = 0.8 change direction passing through a process of suction to injection.

To investigate the influence of transverse Reynolds number R in , other values of R in (0.2, 0.5) that might correspond to ultrafiltration (UF) or microfiltration (MF) systems have been evaluated; in both cases the same behavior was found for the same values of parameter α.

Finally, let us consider the validation of the solute transfer coupled with the fluid motion. The solute concentration profile will follow the exact analytical solution obtained by [START_REF] Haldenwang | Exact solute polarization profile combined with osmotic effects in berman flow for membrane cross-flow filtration[END_REF], if the following conditions are fulfilled: a) α vanishes (negligible pressure drop to get the Berman flow); b) the exact solute boundary layer and Berman flow are initiated at the duct entrance; c) the recovery is low (i.e. the comparison is carried out for small λ).

In Figure 7, the concentration profiles numerically computed, as well as those corresponding to our analytic studies have been plotted for four cases of increasing operating pressures (i.e. P e in ). As expected, increasing operating pressures (i. e. increasing P e in ) leads to higher permeate fluxes and then to a more important concentration polarization. The discrepancy between the numerical and the approximated analytical solution is hardly discernable in the figure 7.

Comparisons with approximated approaches of the literature

Figures 8 and9 show the solute concentration at the wall and the axial variation of the permeate flux, respectively. These figures conduct the comparisons between the results by [START_REF] Song | Theory of concentration polarization in crossflow filtration[END_REF], TSB-plug [START_REF] Song | Concentration polarization in a narrow reverse osmosis membrane channel[END_REF]) and TSB-shear [START_REF] Song | A total salt balance model for concentration polarization in crossflow reverse osmosis channel with shear flow[END_REF]) models described in Appendix (subsection 5.3), with our numerical model for a RO module of We remark that our numerical model predicts a strong development of concentration polarization layer and a more rapid decay of the permeate flux in the length very close to the entrance. This prediction can be explained because with permeable walls, the transverse advection strongly modifies the layer construction, which is no longer created by diffusion only. Our careful iterative process allowed us to take account of this non-usual balance. Let us compare our model with the TSB models. TSB models uncouple the axial and transverse solute transfer, and its solute concentration profile (50) is then flatter than ours (39). Since the total salt concentration is conserved, the predicted wall solute concentration is weaker than ours, and the local permeate flux is higher. If the membrane is long enough, higher fluxes in the TSB models lead finally to higher solute concentrations, and their predicted local permeate flux decreases then faster than ours for Z > 1 m. Recently, Liu et al. ( 2014) developed a TSB model considering a parabolic velocity profile. Their results show an intermediate behaviour between the TSB-plug and the TSB-shear models with a concentration polarization more intense than predicted by the first, but more attenuated than predicted by Comparison of the numerical and analytical solution shows good agreement for low permeate fluxes and short membranes (figure 10). It should be kept in mind that the analytical model has been obtained for high pressure conditions. When the inlet pressure P in is increased, the value of alpha decreases; however, lambda increases at the same time, from 0.0286 for 1 bar to 0.571 for 20 bar. Hence, the hypothesis of λ ≪ 1 is not valid for high permeate fluxes, which explains the gap between the analytical and the numerical solutions. The analytical solution overestimates the permeate flux because it considers a constant transversal concentration profile (39), whereas the numerical solution considers the more realistic development of the concentration polarization from an homogeneous inlet concentration at z = 0.

For long membranes (0.171 ≤ λ ≤ 3.43, Fig. 11), the numerical solution deviates from the analytical solution (even for low permeate fluxes), since the low recovery hypothesis is not valid. Considering that recovery is low implies an underestimation of the solute concentration, an underestimation of the osmotic pressure and finally an overestimation of the local permeate flux. Furthermore, as TMP increases, the discrepancy between the analytical solution and both numerical models is enhanced. This is of course due to the analytical solution, the hypotheses of which are no longer fulfilled, as pressure increases. More precisely, in Fig. 11, the channel length is large and the hypothesis of low recovery becomes invalid, all the more so that TMP produces a large permeation.

Comparison of our numerical solution with the TSB-plug model shows a quite good agreement for long membranes (figure 11). For short membranes (figure 10), the TSB model predicted higher values than our numerical model; this trend was expected due to a major influence, in our numerical model, of the concentration polarization development at the beginning of the process.

Comparisons with experimental data

At this point, we have demonstrated that the Prandtl approximation offers an accurate alternative with respect to the full conservation laws. We now need to validate the overall approach in terms of basic choice with respect to the mathematical model. This must be established by comparison with experiments. For the sake of conducting this comparison, we choose different experimental contributions that propose careful measurements in both NF and RO measurements. We first compare our averaged permeate fluxes as a function of the operating pressure with the experimental results obtained by [START_REF] Geraldes | The effect on mass transfer of momentum and concentration boundary layers at the entrance region of a slit with a nanofiltration membrane wall[END_REF] in a NF module composed by a CDNF501 commercial membrane with an hydraulic permeability of 1.4 × 10 -11 ms -1 Pa -1 , working at operating pressures that range between 10 and 40 bar and for three feed axial rates. Even though the NF module had a length of 20 cm, the reported data only included average values obtained in the first 6 cm. Table 3 allows the quantitative comparison between the data from the experiment and the predictions of the numerical simulation for three aqueous solutions of N a 2 SO 4 , sucrose and polyethylene glycol 1000 (PEG1000). For characterizing the physical properties of the chemical solution, we recall that the numerical prediction only uses the standard constants taken from the literature.

For both solutions with N a 2 SO 4 and sucrose, the comparison between numerical and experimental results shows an excellent agreement (a discrepancy of the order of 10by contrast this trend is reversed for the solutions of sucrose, the numerical simulation providing us with higher values than the experiments.

As for the PEG1000 solutions, our modeling overpredicts the permeation, in particular when the operating pressure is higher than 20 bars. When the feeding flow is increased (i.e. when the concentration polarization is lesser), Solute the discrepancy diminishes. As a matter of fact, the Schmidt number of PEG 1000 is very large and the polarization at high pressure attains important values, in such a way that the osmotic law increasing linearly with PEG1000 concentration can no longer be appropriate. To conclude this stage of model validation, the numerical predictions can be considered as reliable for both solutions of N a 2 SO 4 and sucrose. For solutes of larger Schmidt numbers (as PEG 1000), the polarization at high trans-membrane pressure leads to very rich concentrations, for which the physical "constants" have to be modified.

P in W in C in Sc J exp J num 10 5 Pa ms -1 kg • m -3 ms -1 ms -1 N a 2 SO
We next turn towards a reverse osmosis experiment conducted with salt. More precisely, we now compare two predicitions (the TSB-shear model and the present numerical model) with the experiments of [START_REF] Zhou | A numerical study on concentration polarization and system performance of spiral wound ro membrane modules[END_REF]. This experiment presents a particular difficulty for our laminar simulation, since the flow is turbulent due to the effects of spacers. Actually, an aqueous solution of NaCl feeds a set of four spiral wound RO modules of a polyamide composite membrane (TFC 2540SW, Koch Membrane Systems). The pilot had a total length L of 4 m with a distance between membranes 2d of 0.6 mm, and a membrane resistance I 0 = 8.41 × 10 10 Pa • s/m. Experiments were carried out with a mean inlet axial velocity W in = 0.075 m/s and for three different feed initial concentrations C in : 500, 1000 and 3000 mg/l.

To overcome the difficulty due to the presence of turbulent mixing, it was proposed in [START_REF] Zhou | A numerical study on concentration polarization and system performance of spiral wound ro membrane modules[END_REF]] to adopt the concept of an effective solute diffusivity, which is suggested at the value of D 0 = 1.81 × 10 -8 m 2 /s, i.e. about ten times the molecular diffusion coefficient of sodium chloride. We hence have used the same turbulent diffusivity in our numerical Prandtl approach. Numerical simulations are depicted by solid lines, while TSB model results and experimental data are plotted as dash-dotted lines and symbols, respectively. For all cases both models agreed very well with experimental fluxes until a pressure of 10 bar but deviated gently as pressure still increases. The maximum difference was about 15% for a feed concentration of 500 mg/l at a pressure of 15 bar. Actually, although the simulations have considered a membrane which was 4 m long, the experiments were carried out in 4 modules 1.016 m long connected in series [START_REF] Zhou | A numerical study on concentration polarization and system performance of spiral wound ro membrane modules[END_REF]. There is no permeate flux along the connection pipes between the modules, thus reducing concentration polarization. This could explain the discrepancy between simulations and experiments, especially for the highest trans-membrane pressures.

We now conduct the last comparison of our numerical predictions with experimental results. Measurements of the concentration polarization profile 2010) and consisted in measuring the concentration polarization layer of a solution of N a 2 SO 4 within a RO module by digital holographic interferometry. The module was composed of a thin film membrane (TFM-50, from Hydro Water S.L.) and sized 10 cm×3 mm×10 mm (L × W × H). In all three cases the operating pressure P in was 7.2 bar and the feed initial concentration C in was 8.5 kg/m 3 , for three different inlet mean velocities W in of 0.2, 0.7 and 1.7 cm/s. Unlike all the previous experimental results which considered two parallel permeable walls, these experiments were carried in a channel with a permeable wall and an impermeable one. Hence, the limit conditions for the axial velocity ( 22) have been modified:

w(0, z) = 0 (30)
The figure 13 shows the concentration polarization profiles obtained by our numerical model and TSB model (solid and dotted lines, respectively) compared with the experimental ones (symbols). Measurements were performed at the distance from the inlet Z = 5cm. A good agreement is found between our numerical predictions and the experimental data. We observe that the modeling slightly overestimates the concentration gradient (and hence the permeate flux). The same trend in discrepancy also appears in [START_REF] Salcedo-Diaz | Visualization and modeling of the polarization layer in crossflow reverse osmosis in a slit-type channel[END_REF]]. Here, we can propose two possible explanations. First, our model considers a total rejection by the membrane, although experiments showed a 97 % rejection for N a 2 SO 4 , only. Second, the experimental channel has a width-to-height ratio of 3.3, and the consideration of a twodimensional geometry might be a slightly inaccurate approximation. Lastly, the TSB model meets difficulties in fitting with the experimental concentration profiles. This can partly be attributed to the fact that the latter predictions have been obtained in the case of a symmetric channel.

To finish this important stage of mathematical model validation with respect to experiments, let us make the following remarks. The Prandtl simulations, in a general manner, have shown an excellent accordance with the analytical and numerical previous results, and a good agreement with the experimental results. This point is valid for both "clean water" filtration and coupling between hydrodynamics and mass transfer. Furthermore, the possible discrepancies we have observed when comparing with experiments are likely to be corrected by using concentration-dependent physical constants. Such an improvement does not affect the present overall numerical approach, and will be the object of future contributions. Another remarkable point is that in all cases the flow always match locally with a pattern that can be identified with the analytical solution (i.e. [START_REF] Berman | Laminar flow in channels with porous walls[END_REF]). In other words, the flow pattern adapts to permeation, while mass transfer rules the whole system (of course, the flow plays an important role in mass transfer...).

Conclusions and final discussion

An efficient two dimensional numerical model has been developed; it solves the solute conservation equation coupled with the steady Navier-Stokes equations under the Prandtl approximation. This approach of crossflow membrane filtration takes the pressure-dependent leakage into account, as well as the effect of concentration polarization on osmotic (counter-)effects in RO/NF filtration. More precisely, the effective trans-membrane pressure depends on osmotic effects due to concentration polarization at the membrane. Under specific operating conditions, this numerical model can predict the local permeate flux and the solute concentration polarization with an excellent accuracy, as shown during the validation steps, where several exact In comparison with other numerical approaches, our model predicts a rapid development of the polarization boundary layer and hence a faster decline in permeation in the vicinity of entrance. This prediction that differs from other numerical models might be explained by the ability of the present model to accurately enforce the non-linear boundary condition at the very entrance.

In the last section, we have compared our numerical predictions with several RO/NF experimental results. A good accordance was generally obtained. This validates the overall modeling approach, from the choice of the model to the numerical solution method. It is furthermore worthy of note that the present numerical tool can be easily adapted to more sophisticate situations in terms of viscous and transport properties. In the same vein, the numerical model can be extended to other geometries, as tubular membranes, or to non-linear laws on osmotic pressure, as realized in [START_REF] Lopes | Prediction of permeate flux and rejection rate in ro and nf membrane processes: Numerical modelling of hydrodynamics and mass transfer coupling[END_REF]]. In the same manner, a model of membrane selectivity with partial solute rejection can be considered as developed in [START_REF] Lopes | Predicting permeate fluxes and rejection rates in reverse osmosis and tightnanofiltration processes[END_REF]]. The present numerical approach has easily been adapted to a situation of fouling, a work which has been left to future publication. Lastly, the filtration of multi-species systems (as in the case of sea water desalination) can easily be envisaged from the present work. For each species, a conservation law can be duplicated from the present solute conservation law. Furthermore, an osmotic number must be introduced for every species. The purpose of this addendum is to present the various main properties of the analytical or semi-analytical results that we have used in the validation step of our numerical approach. We shall particularly stress on the underlying assumptions making feasible the analytical approach.

The Berman flow

The Berman flow is a basis solution for the fluid motion in a leaky channel, when the permeation is uniform and fixed all along the porous walls [START_REF] Berman | Laminar flow in channels with porous walls[END_REF]. Under this hypothesis, Berman derived the flow solution {u(x, z), w(x; z)} to the steady Navier-Stokes equations complemented with the lateral boundary conditions u

(x = ±1, z) = U W /U in = ±1, as u(x, z) = B(x; R in ) ; w(x, z) = (1 -z)B ′ (x; in) (31)
where B(x; in) is the solution of the ordinary differential equation (ODE):

R in (BB ′′ -B ′2 ) -B ′′′ = K(R in ) in 0 < x < 1 (32)
with the boundary conditions:

B(0) = 0, B(1) = 1, B ′ (1) = 0, B ′′ (0) = 0 (33)
and R in is the transverse Reynolds number of permeation. For the standard configurations of filtration, a rich literature has shown that the solution to ODE (32) is unique, stable and attractive in the sense that, if the inlet conditions are different from the Berman x-profiles, Berman solution is recovered shortly after the entrance. Moreover, in the range (i.e. R in ≤ 4), the following Taylor expansion of the Berman function rapidly converges (see for instance [START_REF] Haldenwang | Laminar flow in a two-dimensional plane channel with local pressure-dependent crossflow[END_REF]]) in the standard situations of filtration:

B(x; R in ) = n=∞ n=0 1 n! f n (x)R n in ( 34
)
where the first two coefficients are:

f 0 = 3x 2 1 - x 2 3 ; f 1 = x 280 -x 6 + 3x 2 -2 (35)
The Berman flow for filtration only makes sense when the pressure drop remains weak enough all along the channel, in order to justify that the permeation is uniform. This condition corresponds to the so-called assumption of "High Pressure" (HP).

Exact analytical solution for mass transfer in Berman flow

In a previous study, [START_REF] Haldenwang | Exact solute polarization profile combined with osmotic effects in berman flow for membrane cross-flow filtration[END_REF] obtained an exact analytical solution for the solute conservation law, when the carrying flow is of Berman type. Let us now consider the solute mass conservation law (21), in which we consider the velocity components {u, w} as fixed and given by expression (31) where the transverse Reynolds number is unknown, since the permeation will depend on the hindrance due to concentration polarization. Let us denote this unknown permeation with U 0 as the transverse velocity at the wall, leading us to define Re 0 [resp. P e 0 ] , the (unknown) Reynolds [resp. Péclet] number of permeation, as Re 0 = ρ 0 U 0 d/µ 0 [resp. P e 0 = U 0 d/D 0 ]. An exact solution to this problem has been established in [START_REF] Haldenwang | Exact solute polarization profile combined with osmotic effects in berman flow for membrane cross-flow filtration[END_REF] and reads:

C(x, z) = C(x, z) C in = 1 (1 -z) exp P e 0 x 0 B(x; R 0 )dx (36) with x 0 B(x; R 0 )dx = n=∞ n=0 1 n! F n (x)R n 0 ( 37
)
where the first two F n (x) take the form :

F 0 (x) = 3x 2 4 - x 4 8 ; F 1 (x) = 1 280 - x 8 8 + 3x 4 4 -x 2 (38)
Developing this solution (36) at the lowest order with respect to R 0 is sufficient in most RO/NF situations and leads to the simple expression for the solute concentration field.

C(x, z) = C(x, z) C in = 1 (1 -z) exp P e 0 3x 2 4 - x 4 8 (39)
This relationship gives a mathematical form to the phenomenon of concentration polarization, and allows us to obtain the concentration at the membrane wall as:

C W (z) = CW C in ≈ 1 1 -z exp 5 8 P e 0 (40)
Now, let us consider the osmotic pressure provoked by C W (z). To maintain a uniform leakage, the osmotic pressure must be constant all along the channel. Therefore, we must assume that (1 -z) ≈ 1. In other words, the channel length is limited enough to keep the concentration polarization uniform all along the wall. This condition corresponds to the so-called "Low Recovery" condition (LR). Therefore, under the "HP-LR" assumption, concentration polarization and permeation satisfy

where P e in is the pure solvent transverse Péclet number and P e osm in is a dimensionless form of inlet concentration, defined as:

P e in = P in d I 0 D 0 ; P e osm in = iRT C in d I 0 D 0 (43)
Permeate flux U 0 can be calculated from the last implicit equation ( 42) in P e 0 by using an iterative method. This analytical approach is a useful tool to estimate the permeation velocity U 0 in presence of polarization and osmotic effects. Finally, let us note that, in the limit of the high operating pressures P e in ≫ 1, expression (42) reduces to P e 0 ≈ 1.6ln [P e in /P e osm in ]. The latter expression simply demonstrates the classically observed saturation of the permeation when increasing operating pressure. [START_REF] Song | Theory of concentration polarization in crossflow filtration[END_REF] proposed to consider equation ( 21) under various simplified forms, accordingly with the concerned zone. The first step in the assumption process consists in neglecting the second term in the RHS of equation ( 21) (the mass axial transport) when considering the domain in the very vicinity of the membrane. After integration, this yields the dimensional form

Song-Elimelech model

U 0 ( C -C in ) = -D 0 d C dX (44)
where C in is the feed solute concentration or bulk concentration. This model assume that the concentration polarization layer is smaller than the spacing of the membrane channel and considers at steady state a mass balance relationship in the axial direction through the channel:

∞ 0 γX( C -C in )dX = C in Z 0 U 0 (Z ′ )dZ ′ ( 45
)
where γ is the shear rate, Z the axial distance from the entrance and X the transverse coordinate. In this model the wall concentration could be determine by solving equations ( 44) and ( 45) simultaneously

CW = C in 1 + 1 D 2 γ U 2 0 Z 0 U 0 (Z ′ )dZ ′ (46)
If a constant shear rate γ is assumed (Figure 14) an analytical expression of the permeate velocity can be obtained as follows:

U 0 (Z) = U 0 (0) (1 + AZ/L) 1/3    1 1 + AZ/L + 4 1 2 + 2 1 3 - 1 1 + AZ/L + 4 1 2 -2 1 3    ( 47
)
where U 0 (0) is the permeate velocity at inlet, A is a dimensionless parameter of operation and L is the membrane channel length. Considering that the total solute flux downstream along the channel is constant even though both solute wall concentration and permeate flux changed carried us to the following mass balance relationship:

Cm W m = C in W in (53)
where W in and W m are the axial mean velocity at the entrance and the axial mean velocity at any Z position respectively. Working with an average axial velocity is equivalent to assume that the fluid has a constant axial velocity across the channel height (plug flow profile in Figure 14). ¿From mass conservation principle for the solvent we can establish the following relationship

dW m dZ = - U 0 d ( 54
)
The variation of Cm along the membrane channel can be determine by taking derivate of equation ( 53) with respect to Z and using equation ( 54) we obtained:

d Cm dZ = U 0 W m d Cm ( 55
)
A numerical method can solve equations ( 55),( 54), ( 2) and ( 52). Discretizing the membrane channel length L into n equal segments (∆Z = L/n), the finite difference forms of these equations are:

Cmi+1 = 1 + U 0i ∆Z W mi d Cmi ( 56
)
W mi+1 = W mi - ∆Z d U 0i (57)
U 0i+1 = P -Γ CW i I 0 (58)

CW i+1 = C in + ( Cmi+1 -C in ) f U 0i d D 0 (1 -exp(-U 0i d/D 0 ) (59)
The iterative process begins when the feed solute concentration C in , the axial mean velocity at entrance W in and the driving pressure P are given. The mean solute concentration Cm , the axial mean velocity W m , the permeate velocity U 0 and the solute wall concentration CW , can be then obtained progressively throughout the entire membrane channel using the finite difference equations above.

Total salt balance model with shear flow

In this model a shear flow rate γ for axial velocity was considered, as shown shear flow profile in Figure 14 ¿From Figure 14 we induced that the shear rate γ for an average axial velocity W m was:

γ = 2 d W m (64)
The derivative of γ with respect to Z gave:

dγ dZ = 2 d dW m dZ = - 2 d 2 U 0 (65)
Substituting equations ( 64), (65) into equation ( 63) reached the following relationship:

d Cm dZ = U 0 W m d Cm + U 0 W m d A 1 C in (66)
where A 1 is a new function of the transverse Péclet number:

A 1 = P e 0 (1 -exp(-P e 0 )) 2(1 -exp(-P e 0 ) -exp(-P e 0 )P e 0 )

-1 (67)

In this model, the numerical method solves equations (2), ( 52) and (54) coupled with equation (66) instead of equation ( 55) using the following forward finite difference forms: Let us finally remark that the models described above assume a constant driving pressure P through the membrane channel. The variation in operating pressure is obtained by solving the fluid motion. In the latter case, to analytically study this more complex situation is, to say the least, a hard task. Analytical developments for PDE are never easy, especially when hydrodynamics is strongly coupled with mass transfer equation. This is why a numerical approach seems to be more suitable to analyze this situation.

Cmi+1 = 1 + U 0i ∆Z W mi d Cmi + U 0i ∆Z W mi d A 1 C in (68) W mi+1 = W mi - ∆Z d U 0i (
Prandtl model for concentration polarization and osmotic counter-effects in a 2-D membrane channel

Figure 1 :

 1 Figure 1: Sketch of cross-flow membrane filtration and concentration polarization

Figure 2 :

 2 Figure 2: Calculation sequence of the numerical model

Figure 3 :

 3 Figure 3: Axial and transverse velocity profiles obtained by numerical simulation compared to analytical solutions for a uniform permeation case (α = 0.01, Rin = 0.1, clean water): correspond to symmetric channel (in Fig.3a, note that both curves for z = 1 and vertical axis superimpose).

Figure 4 :

 4 Figure 4: Results obtained by numerical simulation for a pure solvent case with α = 0.01 and Rin = 0.05 R in is reflected in all panels of Figures 4, 5, and 6. Panels a and b illustrate

Figure 5 :

 5 Figure 5: Results obtained by numerical simulation for a pure solvent case with α = 0.5 and Rin = 0.05

Figure 6 :

 6 Figure 6: Results obtained by numerical simulation for a pure solvent case with α = 0.75 and Rin = 0.05

Figure 7 :

 7 Figure 7: Comparisons of concentration profiles obtained by both numerical and by analytical approaches for a uniform permeation case in a symmetric channel for various values of inlet Péclet number (α = 0.001, Nosm = 0.1, z = 0.05).

Figure 8 :

 8 Figure 8: Comparisons of solute wall concentration along the channel between Song-Elimelech, TSB-plug and TSB-shear model with our numerical approach for different operating pressures ranging between 6 and 10 bar

Figure 9 :

 9 Figure 9: Comparisons of permeate flux along the channel between Song-Elimelech, TSBplug and TSB-shear model with our numerical approach for different operating pressures ranging between 6 and 10 bar

Figure 10 :

 10 Figure 10: Channel length L = 1 m

Figure 11 :

 11 Figure 11: Channel length L = 6 m

Figure 12 :

 12 Figure 12: Comparisons of the average permeate flux according to Zhou experiments [Zhou et al. (2006)]and permeate flux obtained numerically for different operating pressures and N aCl initial concentrations.

Figure 13 :

 13 Figure 13: Comparisons of concentration profiles between TSB (dotted lines) and present numerical model (solid lines) with experiments [Fernández-Sempere et al. (2010)] (triangles) for a solution of N a2SO4 operating at Pin = 7.2 bar and a feed initial concentration Cin = 8.5 g/l (a)Win = 0.2 cm/s, (b)Win = 0.7 cm/s and (c)Win = 1.7 cm/s

  ratio between Hagen-Poiseuille pressure drop and transmembrane pressure Γ osmotic factor, Γ ≡ iRT (J) γ shear rate in TBS model (1/s) λ channel length reduced by dead-end length L de , λ ≡ L/L de µ 0 dynamic viscosity (kg/m/s) the membrane wall (moles/m 3 ) CW = C W C in ˜transforms variables or certain unknowns into their dimensional form B(x; R in ) Berman function that gives the dimensionless transverse velocity, i.e. the solution to Eq. (32) combined with boundary conditions (33) C in feed solute concentration (moles/m 3 ) D 0 solute diffusivity through the solvent (m 2 /s) i number of dissociated entities (ionic or neutral) per solute molecule I 0 membrane resistance (kg/m 2 /s) L membrane length (m) L de dead-end length, or location where axial flowrate is exhausted, L de = W in d/U in N osm osmotic number, N osm ≡ P osm in /P in P driving pressure (kg/m/s 2 ) p dimensionless pressure (p = P/P in ) -)pressure due to feed concentration, P osm in ≡ ΓC in P e 0 Actual Péclet number of permeation, P e 0 ≡ U 0 d/D 0 P e in pure solvent transverse Péclet number R perfect gaz constant (J/mol/K) (= W/W in ) W m axial mean velocity at Z (m/s) W in axial mean velocity at entrance (m/s) X dimensional transverse coordinate in TSB model (m) x dimensionless transverse coordinate (x = x/d) Z dimensional axial coordinate in TSB model (m) z dimensionless axial coordinate (z = z/L de ) 5. Appendix: a brief description of the models used for validation

Figure 14 :

 14 Figure 14: Different profiles of a laminar flow in the membrane channel

  Total salt balance model with plug flowThis model takes into account that concentration polarization layer could be developed to the whole height of the membrane channel. Solute concentration at any transverse coordinate can be determine by integrating equation (44) C = C in + ( Cm -C in ) the average solute concentration and P e 0 is the transverse Péclet number:P e 0 = U 0 d D 0(51)We can consequently determine the wall concentration CW by evaluating equation (50) in X = 0: CW = C in + ( Cm -C in ) P e 0 1 -exp(-P e 0 ) (52)

0 γX

 0 therefore the mass balance relationship (53) along the cross-flow direction reads:d CdX = W in C in d (60)Replacing the expression (50) for C in equation (60) above leads to:W in C in d = of equation (61) results in:W in C in d = γd 2 1 P e 0 -exp(-P e 0 ) 1 -exp(-P e 0 ) ( Cm -C in )

  i+1 = C in + ( Cmi+1 -C in ) P e 0i 1 -exp(-P e 0i )(71)

  Table 2 summarizes the dimensionless parameters of three tested cases in "clean water" filtration. The influence of parameters α and Test Case P in (bar) W in (m/s)

				α	R in
	4	100	0.23	0.01 0.05
	5	100	11.2	0.5 0.05
	6	100	16.7	0.75 0.05

Table 2 :

 2 Operating conditions of various test cases in "clean water" filtration

Table 3 :

 3 Average permeate flux of experiments carried out by Geraldes et al. (2002) compared with the same quantity obtained by numerical simulations

  B. Bernales, P. Haldenwang, P. Guichardon, N. Ibaseta M2P2 UMR 7340; Aix Marseille Université, CNRS, Centrale Marseille 38, rue Frédéric Joliot-Curie 13451 Marseille Cedex 20, France Highlights 1. An efficient two dimensional numerical model based on the solute conservation equation coupled with the steady Navier-Stokes equations has been developed. 2. The local permeat flux and the solute concentration polarization are predicted with an excellent accuracy. 3. In comparison with other numerical approaches, our model predicts a rapid development of the polarization boudary layer and a faster decline in permeation. 4. An excellent agreement was obtained between our numerical results and several reverse osmosis experimental results.