On the Formation of Unstirred Layer in Osmotically Driven Flow

Osmotically driven flow across a semi-permeable membrane under a constant static pressure difference is revisited with referring to the previous reports for reverse osmosis. A few mathematical techniques for obtaining the approximate solution, such as that for inverse problems used in the field of heat transfer, are presented with an emphasis on the nonlinear boundary condition and the time-dependent solvent flow-rate. It is concluded that the layer is spontaneously formed by osmosis rapidly in the time scaled by $\sim {\rm O}\bigl(\sqrt{t}\bigr)$, and that the layer thickness grows with no upper limit in an infinite time interval. Based on the obtained solution, we will also discuss the thermodynamical output work in an irreversible process which is extracted from the system as an osmotic engine.


Introduction
Osmosis in membrane transport is a physical mechanism underlying a variety of engineering or biological phenomena [3]. In desalination process [4,5], which is nowadays the essential requisite for supplying water to a rising population globally, fresh water is produced by filtration through a polymeric membrane by a large external hydrostatic pressure applied against "reverse" osmosis. The membrane must suffers not only from hydraulic viscous resistance but also from the osmotic pressure depending on the solute concentration on the membrane. Moreover, beyond an engineering expectation, the highly concentrated solution is accumulated in a thin boundary layer near the membrane at the high pressure side, which induces an unfavourable virtual resistance on the water filtration and eventually gives rise to a membrane fouling, which is of great importance for the recent engineering application.
On the other hand, "forward" osmosis means the spontaneous flow through a membrane mainly driven by concentration difference across the membrane subject to a small external hydrostatic pressure difference, which leads to the concentration relaxation attributed to the physics principle, increase of entropy. The forward osmosis is practically utilised in the wide range of all living creatures, for example, cells and organs, such as epithelia of blood capillary or intestinal membrane in animals, or the translocation of nutrition in phloem network distributed in the whole body of plants [6], where any mechanical pumping system such as a heart in animals is absent. The forward osmosis has recently attracted great attention as a new type of power generation or energy recovery system for future engineering innovation [7,8].
Osmotically driven flow. Here, we focus on osmotically driven flow, an ideal configuration of the forward osmosis so as to measure how quickly the osmosis spontaneously drives the net fluid movement under no external pressure difference across a membrane. Suppose that initially pure solvent is partitioned from a solution of non-electrolyte with bulk concentration, c ∞ , by a semi-permeable membrane, for example, at the centre of a U-tube, which is used often to demonstrate the osmotic driven flow for science education. From the microscopic viewpoint, steric barrier effect of the membrane selectively allows only solvent molecule passage but prevents solute molecules from leaking out through membrane pores (Fig.1).
The solvent spontaneously starts to seepage across the membrane from the pure solvent side (trans-side) into the solution side (cisside) due to the osmotic pressure, which virtually originates to the Brownian motion of molecules from the microscopic viewpoint [9,10,11]. This transmembrane volume flux of solvent from trans-side to cisside is sustained unless the equilibrium is achieved where the static hydraulic pressure difference across the membrane cancels the osmotic pressure. According to classical thermodynamics, the osmotic pressure cis-side trans-side Π of the ideal solution is deduced from the balance of the chemical potentials across the membrane, and is proportional to the product of the temperature and the solute concentration difference between cisand trans-sides, which is known as van't Hoff law. However note that, in reality, osmotic pressure, which is the strength of the macroscopic opposite pressure required to prevent the solvent seepage to the cis-side, is not a stress through a bulk fluid in a mechanical sense. Thermodynamics tells us only the relaxation degree and direction to the equilibrium in terms of free energy, but not the time scale and solvent flow velocity in the primary transient process. The seepage is a non-equilibrium and irreversible process which is still be one of challenging themes in modern physics, and we have to invoke the time-dependent continuum mechanics as well as the classical thermodynamics.
As time elapses without increase of opposing hydraulic pressure, the resultant transmembrane volume flux of solvent advects away solutes from the vicinity of the semi-permeable membrane, so that a certain layer of less solute concentration will be formed locally on the membrane in cis-side. The developing layer spatially localised on the semi-permeable membrane has been named as either concentration boundary layer [12,8] or as unstirred layer [13,14]. The layer could be the origin of the aforementioned fouling in the forward osmosis process as well as in the reverse osmosis. In the forward osmosis, the decrease of solute concentration difference across membrane is inevitably followed by the reduction of the effective osmotic pressure and the transmembrane volume flux, although the development of the layer may be suppressed either by mechanical mixing with crossflow or by supply of solute molecules in the bulk region in cis-side.
Objective. The unstirred layer has been accounted to be one of unpredictable or unobservable factors of engineering problems appeared in desalination process, and thus its thickness is implicitly assumed to be saturated. Actually, the presence of the unstirred layer is apparently plausible in practical cases, partly because the steady equilibrium has been mainly pursued in engineering viewpoint, and partly because the thickness of the unstirred layer is too thin to be observed in experiments and otherwise disturbed by an external flow brought by some bulk crossflow tangential to the membrane. Therefore, the detail discussion on the formation time-scale or the developoing thickness of the layer have hardly appeared in the previous literature with a few exceptions [1,2]; how quickly does it develops ? Does its thickness saturate within a finite timescale ? In the present study, taking into account the unsteadiness of flow and concentration next to the membrane, we will revisit the problem to estimate analytically the formation timescale and thickness of the layer, by solving time-dependent governing equation of layer in the aforementioned ideal configuration.

Formulation
In what follows, the solute (number) concentration and flux and the solvent flow velocity are represented asc,j andũ, respectively, where an tilde symbol indicates a dimensional variable. Hereafter, we shall restrict our attention to the case that these quantities are uniform along the membrane from the macroscopic point of view, so that they depend only on the time,t, and the distance from the membrane surface in cis-side,x. We neglect the effect of gravity for ease of discussion, so that the natural convection due to the inhomogeneity of density distribution does not occur. The macroscopic mass conservation rule is written by ∂c ∂t +∇·j = 0. Following Fick's law, the fluxj is given approximately by the sum of advection and diffusion,j =cũ −D∇c, whereD is the solute diffusivity in solvent. The incompressibility of solvent flow may allow us to presume thatũ is independent ofx. We adopt a couple of boundary conditions,j × n = 0 (the impermeability of the membrane with the normal vector n against solutes atx = 0) andc =c ∞ (the uniform bulk concentration atx → ∞). These conditions mean no solute is supplied x > 0. An initial condition will be required for a fully posed initial-value problem.
Negligible inertia. If the solute concentration is so dilute, van't Hoff law is valid,Π =c 0kBT , whereΠ andc 0 are the osmotic pressure and the concentration in the layer in contact with the membrane in the cis-side,c(t,x = +0), respectively. It is empirically known that the solvent seepage velocity (volumetric flux) is proportional to the sum of the osmotic pressure and the static hydraulic pressure differencep 0 across the membrane, under no electrostatic potential difference across the membrane [15]. This proportionality is known as the classic Starling principle of fluid exchange,ũ·n = σL p (Π−p 0 ), where in the present study the reflection coefficient of the membrane for solutes, σ, is unity and the hydraulic conductance (osmotic permeability, filtration coefficient) of the membrane,L p , is supposed to be a constant independent of concentration or property of impermeable solute. The reference ofp 0 is taken at transside.
In general, the hydraulic pressure differencep 0 may be a timedependent variable, for example, a function of the waterlevel at time t, which is calculated as t 0ũ (τ )dτ in case of U-tube. Hereafter,p 0 is kept to steady, and is zero unless otherwise noted. Starling's relation may be practically satisfied in most of experiments, because the relaxation timeconstant calculated from fluid inertia and hydraulic conductance is, in the most of cases, too small to be observable compared to the experimental timescale. We can estimate the relaxation time-constantt relax =ML p /S from the total fluid mass,M , and the membrane area,S. Although these variables vary, of course, largely depending on the membrane thickness as well as material and size of solvent molecule, the magnitude of hydraulic conductance have been investigated experimentally for a variety of membranes, for example, O(10 −11 )[m/s/Pa] for frog mesentery, O(10 −12 )[m/s/Pa] for plasmalemma of Nitella translucens or Viskingdialysis tubing, O(10 −13 )[m/s/Pa] for toad skin in the literature [12,16]. Even if we roughly overestimate fluid mass and membrane area, we find , the time-constant is of the order of 1[µs] at most. Thus the inertia is negligible in a practical sense, so we hereafter presume that the solvent velocity is determined by the Starling's relation.
Assembling all the aforementioned equations, via nondimensionalising time and length (t t 0 → t,x x 0 → x) and via rescaling of dependent variable (c c ∞ → c), we end up with the following second order nonlinear partial differential equation: with a couple of boundary conditions, where c 0 is defined by c 0 = c(t, x = +0), and p 0 =p 0 /(c ∞kBT ). The nonlinearity in the advection term and the boundary conditions acts an flavour of difficulty in the problem. The units of time and length are taken ast 0 =D/(L pkBTc∞ ) 2 ,x 0 =D/(L pkBTc∞ ), respectively, so that the coefficient of advection and diffusion terms are both unity in the present study. This nondimensionalisation was previously introduced by Nakano [1] or Liu [2], although their focus was mainly in the reverse osmosis. The similar but different nondimensionalisation has been adopted in the literatures [17,14], where nondimensionalisation are determined by some representative length or velocity units given by external conditions in individual cases, such as constant thickness of unstirred layer or hydraulic pressure difference across the membrane.
It should be noted that the present nondimensionalisation is the only possible way in our problem with p 0 = 0, where no external stirring factor is involved, that is, neither time nor length scale are given.

Analysis
pseudo steady solution. The obtained deterministic equation is difficult to be analytically integrated with a uniform initial condition, c(t = 0) = 1. In the literature, the steady equilibrium observable in experiments has been mainly pursued, which would be realised independently of initial conditions uncontrollable in experiments. Following Ref. [13,18], let us consider the case that the thickness of the unstirred layer in the equilibrium, ϑ, saturates to a constant within a finite time interval.
Eliminating the time derivative term, we obtain This condition can be satisfied by the following expression regardless of discontinuity of the first order spatial differential, where Θ is the Heaviside step function. The steady concentration profile is steepest at x = ϑ, which is implausible in reality as previously pointed out by Ref. [14]. Although the magnitude of ϑ is determined by the relation c 0 = e −c 0 ϑ deduced from the boundary condition at x = +0, the arbitrariness of c 0 is still remained, which should be determined by introducing an external mixing effect in the bulk region into account. On the other hand, in the present deterministic system without any external mixing in the bulk region, the magnitude of c 0 should be determined selfconsistently. As pointed out by Liu [2], who have conducted experimental observation of the reverse osmosis, the time-dependency is essential to study the present problem.
Numerical solution by finite-difference method. For the purpose to verify that the steady state is not realised in a strict sense, we examined the numerical time-integration of the original time-dependent equation with the aforementioned initial condition and boundary conditions. Fig.2 shows the time series of nondimensionalised concentration profile, which is numerically integrated under forth-order Runge-Kutta scheme. For the reference, c(t n , x) > 0.999 is satisfied for x > 4.62 at n = 0 and for x > 101 at n = 8, which implies that an artificial boundary at x = 200 given in the simulation affects little to the result for a relatively early stage, n ≤ 8. While the concentration in the layer in contact with the membrane, c 0 (t), is reduced to the half of unity abruptly within one time unit (t/t 0 ≤ 1), the decreasing rate of c 0 in time remarkably slows down over the time unit (tt 0 ≥ 1). In our system where no length scale exists except forx 0 , the steady state does not seemed to be realised. As time elapses further, the front of the unstirred layer proceeds far away from the membrane and the thickness of the layer increases asymptotically with no upper limit, which means the unstirred layer is not a boundary layer but can be thicker than literally would be envisaged [6]. The finiteness of characteristic thickness of the layer for the reverse osmosis was comprehensively discussed in the appendix of Ref. [2]. In other words, the unsteadiness of the system is related to the fact that the thickness of the unstirred layer is intrinsically undetermined. However, of course, a finite length of the chamber or a finite measuring time interval in experimental restriction would make a steady state to be apparently realised.
For comparison, the pseudo steady profile satisfying both conditions, Eq.(3)c 0 = e −c 0 ϑ in case of c 0 = c(t 5 , 0), is additionally indicated by a dashed curve in the figure. The discrepancy between the dashed curve and the solid curve of n = 5 implies that the solute concentration on the membrane can be exaggerated (overestimated) by a simple extrapolation based on the macroscopic concentration measured at a distance far from the membrane under the assumption of steady equilibrium.
Polynomial profile approximation. The boundary condition at x = +0 is not satisfied by the initial condition. The above numerical demonstration suggests that the unstirred layer develops in a relatively short time interval so as to regularise the singularity on the membrane rapidly. In order to refine the above expression of the concentration, suppose that the thickness of the layer, ϑ, and the local concentration on the membrane in cis-side, c 0 , both are dependent on time, t, instead of Eq.(3). Moreover we assume that the profile of concentration, c(t, x), is expressed by the N -th order polynomials of x (N ≥ 1) in the layer in contact with the membrane, 0 < x < ϑ(t) : This expression is a refined expression of Eq.(3). Furthermore, we will assume the profile to be continuous and smooth enough at x = ϑ, if saying in a strict sense, of differentiability class C N −1 . Thus, we obtain r n = (c 0 (t) − 1)ϑ(t) −N for n = N ; otherwise r n = 0 (0 ≤ n ≤ N − 1). Imposing that the governing equation is satisfied at x = 0, we can deduce a first-order ordinary differential equation for c 0 (t), which can be integrated analytically for any N , and thus leads to for p 0 = 0. From the boundary condition at the origin, we can deduce , that is, the thickness of the layer, ϑ(t), will be later determined from c 0 (t). The present problem originates in the quadratic second-order differential equation, Eq.(1), so one would be interested specially in the case of N = 3, under applying the continuity conditions at x = ϑ(t) for the first and second spatial derivatives and c itself for an arbitrary t. The polynomial profile approximation c 0 (t n ) obtained for a finite N provides a qualitatively plausible time series of the concentration at the vicinity of the membrane surface, but quantatively distinct from the numerical result c 0 (t 0 ) = 0.50; for example, c 0 (t 0 ) ≈ 0.44 for N = 3.
Moreover, one would be interested in the fact that the N -th order polynomial form addressed above converges to the simple expression at the limit N → ∞, in spite of ϑ → ∞. Taking into account the original definition of exponential function, lim N →∞ (1 + α/N ) = e α , one will obtain the limit of the profile of the concentration at time t, In the previous studies, the accurate measurement of the solute concentration at the vicinity of the membrane is challenged with the aid of the extrapolation of concentration. As a measure of the thickness of the layer in the literature, ∆x is defined by ∂c ∂x | x=0 · ∆x = 1 − c 0 . In the present limit, ∆x is also the function of time, ∆x = √ 2t + 2t. Fig.3 shows the time series of the nondimensionalised solute concentration c 0 (t) at the origin given by Eq.(5), compared to the numerical solution.
The comparison between the numerical result and Eq.(5) in the figure suggests that a quantitative agreement is limited within the initial stage of the layer formation, t ≪ t 0 . However, we would like to note that Eq.(5) may be a lower bound of the numerical solution at least. Taking into account that the nondimesionalised solvent velocity equals to c 0 − p 0 , the solvent volume seepaging for an infinite time interval may be calculated as Additionally, we note that the curve of c 0 (t n ) obtained from a finite N exists between the numerical simulation and the limit case of N → ∞, for example, c 0 (t 0 ) ≈ 0.44 for N = 3, and 0.46 for N = 2 (not shown in the figure).
Nakano et al presented an analogous approach using a polynomial expression in the context of the reverse osmosis in Ref. [1,2], where the authors approximate the concentration profile as a cubic polynomial with the time-varying thickness of boundary layer. Substituting the polynomial expression Eq.(4) with N = 3 into the integration of Eq. (1) in 0 < x < ϑ(t), they obtained ∂ ∂t so as to deduce another first-order ordinary differential equation for c 0 (t), instead of imposing that the governing equation is satisfied at x = 0. Following their procedure, we obtain the ordinary equation for an arbitrary N , and the implicit solution for p 0 = 0 is For comparison, the above implicit expression of c 0 on t for N = 3 is indicated as the dash-dotted curve in Fig.3. Approach in inverse problem. Finally, we will apply a procedure in the inverse problem to the present problem, Eq.(1) and Eq.(2). With leaving c 0 (t) unsolved provisionally, we transform the frame of reference x → ξ, where x = ξ − ξ 0 (t) and ξ 0 (t) is determined by By the coordinate transformation, Eq.(1) is converted into the standard linear diffusion equation in the ξ−t space (−∞ < ξ < ∞ and 0 < t < ∞), which is solvable analytically for a well-posed initial condition in terms of appropriate Green function. It should be however noted that the initial condition in ξ ≤ 0 is arbitrary for the moment, while c(t, ξ)| t=0 = 1 is given for ξ > 0(see Fig.(4)). Thus, for example, we may assume the initial condition expanded in the −∞ < ξ < ∞ is represented as ∞ n=1 a n ξ n , in terms unknown coefficients, a n , so that lim where (τ, η) = ( √ t, ξ √ 4t ), and I 0 (τ, η) = 1 2 erfc(η) , By the way, taking into account that c 0 (t) is a function of √ t in the approximation discussed before (see Eq.(5)), we may suppose that ξ 0 (t) determined by Eq.(6) is expanded into a series as ξ 0 (t) = which will provide a family of restriction on unknown coefficients of the series, a n (n ≥ 1) and b n (n ≥ 3). The remained boundary condition corresponding to Eq.
(2), which should be satisfied at ξ = ξ 0 (t), is in the (τ, η) space. Eq. (8) and Eq.(9) subject to Eq.(7) constitute a set of algebraic recurrence relations for unknown coefficients of the series, (a n , b n ). The procedure may be classified into the so-called inverse problem, which does not guarantees whether the solution converges properly, in general. The coefficients for the case of p 0 = 0 can be determined as It is numerically confirmed that the obtained form can be convergent only within the relatively small radius of convergence. The convergence may be envisaged from the separation point around (n, c 0 ) = (−2, 0.65), where the upper and lower thick dashed curves obtained by the truncation numbers 10 and 11 detaches each other in Fig.(3). Following a mathematical procedure proposed in Ref. [19], we furthermore perform an analytic continuation by so as to improve the rate of convergence of the obtained alternating series. In fig.3, the intermidiate dashed curve designates the expression obtained from a procedure of inverse problem, which provides an excellent agreement with the numerically integratated solution.

Concluding remarks
In the present study, we revisited a nonlinear partial differential equation describing the osmotically driven flow across a semi-permeable membrane under a constant static pressure difference, which was previously investigated for reverse osmosis [1,2]. Firstly we confirmed the three  Table 1. Coefficient of the series (a n , b n , q n ) solved in case of p 0 = 0. The convergence of a n and b n is not excellent and provides the small radius of convergence. While b n and c n are alternating series, a n is also potentially alternating series because a n is defined in Eq.(7) originating at algebraic series of ξ n for ξ < 0. In the present study, keeping in mind that the phenomenon is essentially time-dependent, we applied a few mathematical procedures to the present problem for the purpose to measure the timescale for concentration boundary layer on the membrane to develop. Specially, the last procedure for the inverse problem is successful to provide an analytical expression of the time series of solute concentration, Eq.(10) with Table 1 for p 0 = 0 case. From the result, we found that the localised concentration on the membrane decaying initially as −t 1 2 and eventually as t − 1 2 (cf. [20,21]). This implies that the unstirred layer thickness ϑ diverges with increase of t with no upper limit and the layer intrinsically never saturates.
Here, we will refer the present system to an osmotic engine that can convert osmotic energy to mechanical work. The solvent volume seepaged into cis-side is calculated as |ξ p 0 (t)| = t 0 (c 0 (τ ) − p 0 )dτ per an area on the membrane. The nondimensionalised work output δW for the initial time interval t is a function of t and p 0 , which equals to the product p 0 and |ξ p 0 (t)| because the present seepages is an isobaric process in thermodynamics. Note that the process under the condition p 0 = 0 corresponds to the irreversible free expansion, and results only in increase of entropy (no work output δW = 0), which may complete in a relatively short time. On the other hand, the flow under the condition p 0 = 0 corresponds to the reversible and quasistatic process, which requires a longer time interval to be completed. From these facts, we may conjecture that the maximum power δW/t is realised at a optimal nondimesionalised pressure p 0 (0 < p 0 < 1). Based on Eq.(6) and c 0 (0) = 1, we obtain for a small t, From the equation, we can conclude that the initially instantaneous output power is the maximum 0.25 at p 0 = 0.5, and that the average output power for a longer t is less than 0.25 because ∂c 0 ∂t | t=0 < 0. The integration of Eq.(10) based on the numerically obtained c n shows that δW/t decreases to 0.144 at t = 1. Moreover, we can predict that the optimal value of p 0 for a longer t is less than 0.5 because ∂ 2 c 0 ∂p 0 ∂t | t=0 < 0.
Finally we note that the initial condition in the present study is not artificial. For instance, the solvent flow can be initially prevented by the static pressure equal to the osmotic pressure, p 0 = c ∞ , then the initial condition in the present study is easily realised by extingushing p 0 . In physiology, it is known that such a sudden change of the solute concentration (or the static pressure) gives rise to a certain dysfunction on the membranes of cells. With the aid of some mechanics, cells or tissue in our body are able to respond such a osmotic shock in the short formation timescale of the unstirred layer neutralizing the influence.