# Critical Review of Fluid Flow Physics at Micro- to Nano‐scale Porous Media Applications in the Energy Sector

Authors: Harpreet Singh; Rho Shin Myong
Journal: Advances in Materials Science and Engineering
 (2018)

---

## Abstract

While there is a consensus in the literature that embracing nanodevices and nanomaterials helps in improving the efficiency and performance, the reason for the better performance is mostly subscribed to the nanosized material/structure of the system without sufficiently acknowledging the role of fluid flow mechanisms in these systems. This is evident from the literature review of fluid flow modeling in various energy-related applications, which reveals that the fundamental understanding of fluid transport at micro- and nanoscale is not adequately adapted in models. Incomplete or insufficient physics for the fluid flow can lead to untapped potential of these applications that can be used to increase their performance. This paper reviews the current state of research for the physics of gas and liquid flow at micro- and nanoscale and identified critical gaps to improve fluid flow modeling in four different applications related to the energy sector. The review for gas flow focuses on fundamentals of gas flow at rarefied conditions, the velocity slip, and temperature jump conditions. The review for liquid flow provides fundamental flow regimes of liquid flow, and liquid slip models as a function of key modeling parameters. The four porous media applications from energy sector considered in this review are (i) electrokinetic energy conversion devices, (ii) membrane-based water desalination through reverse osmosis, (iii) shale reservoirs, and (iv) hydrogen storage, respectively. Review of fluid flow modeling literature from these applications reveals that further improvements can be made by (i) modeling slip length as a function of key parameters, (ii) coupling the dependency of wettability and slip, (iii) using a reservoir-on-chip approach that can enable capturing the subcontinuum effects contributing to fluid flow in shale reservoirs, and (iv) including Knudsen diffusion and slip in the governing equations of hydrogen gas storage.


---

## Body

## 1. Introduction

World’s energy demand is expected to grow approximately by 25–37% between 2014 and 2040 per some estimates [1, 2]. In order to secure the future energy supply while also reducing the global carbon footprints, new avenues are being explored to produce energy that is efficient and clean. Some major technologies that are either already contributing to this objective or can potentially contribute in the future include (i) electrokinetic energy conversion devices that power wide variety of applications, (ii) water desalination using synthetically created nanosized membranes, (iii) extraction of hydrocarbons from ultra-tight shale formations, and (iv) storage of hydrogen. Optimum utilization of these technologies requires understanding of fluid transport at micro- to nanosized pores and developing theoretical models that capture the physics of fluid flow at these scales to predict their performance. Understanding the material of the medium together with the physics of fluid flow is important because they are often coupled together in addition to the importance of their role independently; for instance, properties of a rock matrix (material) are important in shale to locate regions with hydrocarbons as well as in designing a hydraulic fracturing job that is strongly affected by the rock’s material properties. As another example, properties of the material are used to control the performance of electrokinetic devices. While it is evident from the literature that more emphasis is given to studying the material of these systems, similar emphasis must be given to study the physics of fluid flow. Incomplete or insufficient physics for the fluid flow can lead to untapped potential of these applications that can be used to increase their performance.This paper reviews the current state of research for the physics of gas and liquid flow at micro- and nanoscale and identified critical gaps to improve fluid flow modeling in four different applications related to the energy sector. The review for gas flow focuses on fundamentals of gas flow at rarefied conditions, the velocity slip, and temperature jump conditions. The review for liquid flow provides fundamental flow regimes of liquid flow, and liquid slip models as a function of key modeling parameters. The four applications from energy sector considered in this review are (i) electrokinetic energy conversion devices, (ii) membrane-based water desalination through reverse osmosis, (iii) shale reservoirs, and (iv) hydrogen storage, respectively. The structure of this paper is divided into three sections: gas flow, liquid flow, and relevant applications in the energy sector, respectively.


## 2. Gas Flow

The flow of gas in micro/nanoscale medium is usually a function of the ratio of the molecular mean free path of a gas molecule to characteristic length of the medium; this function is referred to as theKnudsen number. In such mediums, collisions between a molecule and wall dominate over intermolecule collisions, which causes each molecule to act independently and control the gas properties [3], a condition often referred as rarefied gas flow. As a result, gas slips over the wall surface in case of flow, whereas in case of heat transfer this leads to jump in temperature. Knudsen’s pioneering experimental and theoretical work [4, 5] on rarefied gas flows showed that the compressible form of the Navier–Stokes–Fourier (NSF) model [6–8] is inadequate in describing rarefied gas flows subject to the no-slip boundary condition. One of the challenges in modeling these systems is that the assumption of continuum scale does not apply because the characteristic length scale of the medium approaches the mean free path (MFP) of the molecules.


### 2.1. Fundamentals

#### 2.1.1. Rarefaction and Slip Effect

Rarefaction effect in microsystems is a function of mean free path (MFP) of the gas. Usually the degree of rarefaction effect is characterized by the Knudsen number, which is defined as the ratio of the molecular mean free path (MFP) of a gas molecule to characteristic length of the medium:(1)Kn=λs,λ=MFP=kBT2πσ2p,where s is characteristic length of the medium, σ is diameter of a molecule, kB is the Boltzmann constant, T is absolute temperature, p is pressure, and λ is MFP.TheKn provides a direct indication whether the continuum approach for a flow of gas in a medium is applicable or not. Beyond the continuum approach, the flow is no longer near equilibrium; the assumptions of the no-slip boundary condition, thermodynamic equilibrium, and linear relation between the stress and the rate of strain fail. As the Kn increases, the rarefaction effect becomes more prominent which results in the following:(i)
Flow slippage at the wall(ii)
Flow is no long in thermodynamic equilibrium(iii)
Transfer of momentum from the bulk fluid region to the walls(iv)
Increase in mass flow rate(v)
The appropriate flow and heat transfer models depend on the range ofKn


#### 2.1.2. Compressibility Effect

A change in pressure causes the gas properties to change in order to satisfy mass balance. For example, if the density of the gas decreases due to decrease in pressure, then it causes the gas to accelerate. Colin [9] recommended that compressibility effect is important for gas microflows when Mach number (Ma) is greater than 0.2, while Li and Guo [10] demonstrated experimentally that for an average Ma<0.3 the effect of compressibility is negligible. The compressibility and rarefaction phenomena in gas microflows are generally coupled together, and they can have an opposite effect on each other in the slip flow and early transition regimes based on the magnitude of Reynolds number (Re):(2)Kn=cπγ2MaRe,where γ is the ratio of specific heats of a gas. In the above expression, c=1 as given by Colin [9], whereas c=λ/s in Zhang et al. [11].Tang et al. [12] demonstrated experimentally that the compressibility effect causes the gas flow characteristics to deviate from conventional laminar flow at high Re regime (about Re=2000). In general, it has been demonstrated [13, 14] that the compressibility effect results in negative curvature of the pressure distribution, while the rarefaction effect weakens the negative curvature.


#### 2.1.3. Flow Regimes

On the basis of the Knudsen number, different flow regimes with their physics and applicable models can be defined as shown in Table1. In the case of rarefied gas flow, continuum models are applicable for Kn≤0.01, and for Kn≥10 free-molecular models should be used. However, neither continuum models nor free-molecular models can be used for Kn lying in the intermediate range between the continuum and free-molecular flow regimes.Table 1
Different flow regimes as a function of Knudsen number.
Kn
Flow regime
Physics
Applicable models
Lower bound
Upper bound
0
10–2
Continuum (no slip flow)
Thermodynamic equilibrium and no-slip at the boundary
(1) Navier–Stokes equations, (2) Euler equations, (3) LBM
10–2
10–1
Slip flow
Nonequilibrium effects dominate near the walls. Assumptions of no-slip boundary condition, thermodynamic equilibrium, and linear stress-strain relationship fail
(1) Navier–Stokes equations with slip velocity and temperature jump boundary conditions, (2) Boltzmann gas-kinetic equation, (3) DSMC, (4) Burnett equations, (5) LBM, (6) gas-kinetic scheme, (7) method of moments
10–1
10
Transition flow
Rarefaction effects dominate and slip models become more complex. Assumptions of no-slip boundary condition, thermodynamic equilibrium, and linear stress-strain relationship fail
(1) Navier–Stokes equations with slip velocity and temperature jump boundary conditions, (2) Boltzmann gas-kinetic equation, (3) DSMC, (4) Burnett equations, (5) LBM, (6) gas-kinetic scheme, (7) method of moments
10
∞
Free-molecule flow
Collisions between gas molecules and boundary surface become dominant compared to intermolecular collisions. Assumptions of no-slip boundary condition, thermodynamic equilibrium, and linear stress-strain relationship fail
(1) Boltzmann gas-kinetic equation, (2) DSMC method
Flow regimes as a function of medium and flow parameters.Using the general definition of the Knudsen number that is defined as the ratio of the MFP of a gas molecule to characteristic length of the medium, and the mathematical formulation of MFP, we can writeKn as a function of the characteristic length of the flowing medium (s) and the average flowing pressure in that medium (p) as follows:(3)Kn=kBTs2πσ2p=fs,p.TheKn is estimated from the above equation for various characteristic lengths and average flowing pressures, and the corresponding flow regimes using the values of Kn are plotted as a function of characteristic lengths and flowing pressures as shown in Figure 1.Figure 1
(a) Variation ofKn values as a function of channel height and pressure, and (b) graphical representation of corresponding gas flow regimes. Here, Cont., Slip, Trans., and Mol. represent continuum flow regime, slip flow regime, transition flow regime, and molecular flow regime, respectively.
(a)
(b)Another measure ofKn was proposed by Tsien [15] that was based on thickness (δ) of the boundary layer on the flowing medium as Kn=λ/δ, where δ∝s/Re. Here Re is Reynold’s number for the characteristic length s and is defined as Re=ρus/μ, where ρ,u,andμ are fluid’s density, velocity, and viscosity, respectively. Using these relationships, Kn for the boundary layer thickness (δ) can be defined in terms of the freestream Mach number (M) and Reynold’s number as follows:(4)Knδ=λδ=MReδ∼MRe.TheKn is estimated from the above equation for various M and Re, and the corresponding flow regimes using the values of Kn are plotted as a function of M and Re as shown in Figure 2.Figure 2
(a) Variation ofMa/Re values as a function of Re and Ma, and (b) graphical representation of corresponding gas flow regimes. Here, Cont., Slip, Trans., and Mol. represent continuum flow regime, slip flow regime, transition flow regime, and molecular flow regime, respectively.
(a)
(b)It must be noted that Tsien’s parameter (Knδ) reduces to the ordinary Kn defined for the characteristic length in the free-molecular flow regime where Re≪1, such that δ≈s that results in Knδ=Kn.For flows with large amounts of nonequilibrium, modifying boundary conditions of continuum equations to account for slip will not eliminate all sources of errors. For example, the pressure tensor in more rarefied flows turns anisotropic, whereas the derivation of continuum equations assumes an isotropic pressure field [16, 17]. In those cases, it is recommended to use higher-order equations such as Burnett [18] equations and super-Burnett equations [19]. However, in some specific cases N-S equations with slight correction can be used to match data for a high Kn number regime (0.04–8.3) as shown by some authors [20, 21].


#### 2.1.4. Accommodation Coefficients

The accommodation coefficients describe the macroscopic estimate of the fraction of momentum and energy transfer between the impinging/incident gas molecules and the boundary wall surface; therefore, their values vary between unity (complete accommodation of gas molecules with no reflection) and zero (complete reflection of gas molecules with no accommodation on the solid surface). These two accommodation coefficients for momentum and energy transfer [22] are key parameters in the estimation of slip velocity and temperature jump, respectively. A smaller value of accommodation coefficient would lead to higher velocity slip and temperature jump, because in that case the gas molecules are less affected by the solid surface.(1) Tangential Momentum Accommodation Coefficient. This coefficient defines the fraction of momentum transfer between the incident gas molecules and the boundary wall surface:(5)σu=τi−τrτi=1−τrτi,where σu is tangential momentum accommodation coefficient (TMAC), τi is the tangential momentum of the incident gas molecules, and τr is the tangential momentum of the reflected gas molecules. It can be seen that σu=1 if the reflected component of the tangential momentum of the gas molecules is zero (or fully diffusive), and σu=0 if the tangential momentum of the reflected molecules is equal to the tangential momentum of the incident molecules (or fully reflective/specular).Agrawal and Prabhu [23] presented the values of TMAC as a function of gas rarefaction magnitude (Kn) by using the data reported by various researchers in the literature for both monoatomic and diatomic gases. By fitting a curve to this data, they obtained and proposed the following expression for the TMAC:(6)σu=1−log1+Kn0.7.(2) Thermal Accommodation Coefficient. This coefficient defines the fraction of kinetic energy transfer between the incident gas molecules and the boundary wall surface. For a molecular beam type measurement, thermal accommodation coefficient is defined as follows:(7)σT=Ei−ErEi−Es,where σT is thermal accommodation coefficient (TAC), Ei and Er are the kinetic energies (for the normal component of velocity) of incoming and reflected gas molecules, respectively, and Es is the kinetic energy that would be produced by a diffusive reflection at the temperature of the solid surface Ts. Another definition of TAC that is directly related to measurement of temperatures is defined as follows:(8)σT=Ti−TrTi−Ts,where Ti, Tr, and Ts are the temperatures of the incident gas molecules, reflected gas molecules, and solid surface, respectively.Song and Yovanovich [24] presented the values of TAC as a function of gas molecular weight (Mg) and Ts by using the data reported by various researchers in the literature for both monoatomic, diatomic, and polyatomic gases. By fitting a curve to this data, they obtained and proposed the following expression for the TAC:(9)σT=expC0Ts−T0T0Mg∗C1+Mg∗+1−expC0Ts−T0T02.4μ1+μ2,where C0=−0.57 (dimensionless), C1=6.8g/mole, Ms is molecular weight of solid, μ=Mg/Ms, T0=273K, and Mg∗ is defined as follows:(10)Mg∗=Mg,for  monoatomic  gas,1.4Mg,for  dia/polyatomic  gases.


#### 2.1.5. Knudsen Layer

As the flow deviates from the continuum regime, a sublayer with thickness of few mean free paths (MFPs) starts to develop from the wall of the medium such that it starts to become dominant between the bulk of the fluid and the wall of the medium (Figure3). The molecules within the KL collide more with the boundary wall than with other molecules [25], and this results in fluid (i) having nonzero velocity at the boundary surface and (ii) exhibiting non-Newtonian properties, thus introducing discontinuity in fluid properties between the bulk fluid and the KL.Figure 3
Schematic of the Knudsen layer near a boundary wall. The property in KL is discontinuous across which certain slip and jump conditions expressing conservation of momentum, mass, and energy are applied. The KL region is in nonequilibrium, whereas the bulk fluid flow region is in equilibrium.The flow in theKnudsen layer (KL) cannot be analyzed with conventional N-S equations, which does not account for certain slip and jump conditions expressing conservation of momentum, mass, and energy. The flow within the KL is typically regarded as being governed by the linearized Boltzmann gas-kinetic equation (LBE) [26–29] or by direct simulation Monte Carlo (DSMC) method [30–32]. Velocity slip results in increasing the mass flow rate by approximately 70%, and about one-third of this increase in flow rate is contributed by the non-Newtonian structure of the KL [33]. Therefore, accurate characterization of the KL is important in modeling of gas flow through microsystems.(1) Description of Flow in Knudsen Layer. Accurate simulation of rarefied flow effects requires accurate description and modeling of the KL. For this, it is important to know the performance of parameters at the nonequilibrium state within the KL. The flow within the KL can be described by three different approaches: (i) wall-function method, (ii) higher-order continuum method, and (iii) power-law method.The wall-function [33–35] is independent of the accommodation coefficient and is valid for planar surfaces up to Kn of 0.1. This function predicts the velocity gradient at the wall to be 1.7.Higher-order continuum models [33, 36] derived from the linear Boltzmann equation (LBE) [18, 19, 33] and truncated at disparate orders have been used to describe the flow in KL. One of the motivations for the development of these higher-order hydrodynamic models to simulate rarefied flows was to reduce the prohibitive computational demand of the DSMC method [37, 38]. However, there are questions regarding validity of these models [39] and in accurately capturing the KL [33, 40]. Additionally, these models cannot treat the boundary conditions at the wall [41, 42].Power-law model proposed by Lilley and Sader [40, 43] can be used to accurately describe the flow within the KL, and they suggest that this method describes a general physical phenomenon. The power-law model was derived by using the solution of LBE and DSMC, and it depicts the physics of velocity gradient not captured by the wall-function method [44] and the higher-order continuum method [33].(2) Thickness of Knudsen Layer. Lockerby et al. [33] presented the thicknesses of KL calculated from eight different models that used three different solution approaches of kinetic theory, higher-order continuum models, and molecular dynamics (MD) simulation. The thickness of KL obtained from these eight models varied between 0.9 λ and 4.9 λ. Gusarov and Smurov [45] presented the following expression to estimate an approximate value of KL thickness lc:(11)lc=kBTπd2p,where kB is the Boltzmann constant, T is the absolute temperature, and p is the gas pressure.


### 2.2. Slip and Jump

Knudsen’s pioneering experimental and theoretical work [4, 5] on rarefied gas flows showed that the compressible form of the Navier–Stokes–Fourier model [6–8] is inadequate in describing rarefied fluid flows subject to the no-slip boundary condition. A common way to model the flow of rarefied gas is to modify the no-slip and no-jump boundary conditions in continuum models by revising the basic equations of momentum, mass, and heat transfer. There have been a number of approaches to account for slip boundary conditions through the Maxwell slip (1878), second-order slip [46], and Langmuir slip [47–49]. Maxwell boundary slip condition is the first-order approximation from kinetic theory of gases that has been widely used along with its derived forms in the past 135 years. One challenge in most gas flows at micro- and nanometer scale is the wide range of Knudsen numbers and the lack of a general slip model that can be applicable from the slip regime to the transition regime. Under this scenario, a different approach using molecular models that simulate the gas flow at the molecular level can be adopted. For example, one of most commonly used molecular simulation models is the Boltzmann gas-kinetic model [26–28], as well as many approximate methods based on particle and moments [50] like direct simulation Monte Carlo (DSMC) method and Burnett equations set [18, 26], respectively.Several studies on theory [16, 46, 51–53], Zhang et al. [11], and experimental data [54–57] of gaseous slip flow have cumulatively added to the knowledge base of gaseous slip. Below, we briefly review the fundamentals of velocity slip and temperature jump in terms of their governing equations and provide references for further discussions on these topics.


#### 2.2.1. Velocity Slip

Although the concept of slip was first proposed by Navier in his model in which the magnitude of the fluid velocity is proportional to the magnitude of the flow shear rate, Maxwell [58] was the first person to quantify the slip length of a gas flowing over a solid boundary surface. According to Maxwell [58], for a dilute, ideal monoatomic gas, the first-order slip condition can be described as follows:(12)uslip=ut−uw=2−σuσuλ∂ut∂zw+34μρT∂T∂xw,where σu is TMAC, uslip is the gas slip velocity at the wall, ut is the gas velocity in the tangential direction, λ is the mean free path of the gas molecules, uw is the velocity of the solid wall, ρ and T are the density and temperature of the gas, respectively, and subscript w represents wall surface. The directional axes are represented as x,y,andz with respect to the direction of the flow. The second term on the right hand side represents the effect of thermal creep due to an axial temperature gradient, and it provides slip to the velocity at the wall in the presence of varying temperature.In case of a medium with curvature, the gas slip velocity at the wall can be expressed by adding an additional derivative of the gas velocity inz-direction as follows [59, 60]:(13)uslip=ut−uw=2−σuσuλ∂ut∂z+∂uz∂xw+34μρT∂T∂xw.Colin [61] suggested that for significantly rough surfaces, expression for slip velocity on curved surfaces is equally applicable. Other authors who studied the effect of roughness on rarified gas flow through microchannels [62, 63] made three important conclusions regarding the impact of surface roughness on the flow behavior of gas: (i) the flow behavior in the transition regime is more sensitive to the roughness height than that in the slip flow regime, (ii) influence of roughness height on the flow behavior is more significant than that of fractal dimension, and (iii) presence of surface roughness reduces the gas slip at the boundary wall irrespective of the flow regime. Pelević and van der Meer [64] studied the effect of surface roughness on heat transfer in microflows and reported that the roughness has a marginal effect in increasing heat transfer in microflows. For further readings on various approaches to model gaseous slip flow, readers can refer to the reviews by Cao et al. [16] and Zhang et al. [11].(1) Measurement of Gas Slip Velocity. Despite the presence of different correction coefficients [9] to account for slip flow at the boundary, there is no consensus on the most appropriate equation that can account for gas slip [65]. Additionally, accurate values of TMAC are generally not well known in the literature as it is a function of type of gas, material, and degree of rarefaction. These limitations are usually accounted by experimental data on slip length of rarefied gas flows that is defined by Maxwell’s equation in its simplest form for an isothermal surface as follows:(14)uslip=2−σuσuλ∂ut∂zw=b∂ut∂zw,where b=2−σu/σuλ, which is called as the slip length of rarefied gas flow. Maali et al. [65] reviewed the experimental techniques used to measure the slip length and TMAC for gas flow on isothermal surface. These techniques are Millikan experiment, rotating body techniques, mass flow rate measurements, and atomic force microscopy method, respectively. The gas slip is measured indirectly from these experimental techniques by using slip length, which is deduced from the measurement of the velocity of the drop in the Millikan method, the deflection of the cylinder in rotating body technique, the mass flow rate in rarefied gas flow through microchannels, and the force in atomic force microscopy method, respectively. For detailed discussion on these experimental techniques and their use in measuring gas slip, readers can refer to the review by Maali et al. [65].


#### 2.2.2. Temperature Jump

Similar to the velocity slip expressions presented above, temperature jump boundary conditions can be written as follows:(15)Tjump=Tg−Tw=2−σTσT2γγ+1λPr∂T∂zw,where σT is thermal accommodation coefficient, Tjump is the temperature jump at the wall, Tg is the gas temperature, Tw is the temperature of the solid wall, γ is specific heat ratio, and Pr is the Prandtl number. For further discussions on temperature jump, readers can refer to the reviews by Sharipov [56] and Colin [61].


#### 2.2.3. Second-Order Velocity Slip and Temperature Jump

The linear velocity slip equation (12) and temperature jump equation (15) were obtained from the original (nonlinear) models in which the velocity slip and temperature jump are assumed directly proportional to the degree of nonequilibrium near the wall surface [60]:(16)uslip=ut−uw=2−σuσuλΠμw+34γ−1γ/PrQxpw,Tjump=Tg−Tw=2−σTσT2γγ+1λPrQzkw.Only after the degree of nonequilibrium is taken as linear with the first-order accuracy (that is, shear stressΠ∼∂ut/∂z, tangential heat flux Qx∼∂T/∂x, and normal heat flux Qz∼∂T/∂z), the models are reduced to the well-known first-order models (Equations (12) and (15)). Therefore, if the linear constitutive relations are replaced by higher-order relations like the second-order constitutive relations, the resulting models become a second-order velocity slip and temperature jump. Such models were recently developed by Myong [66] and were applied to the theoretical analysis of the gaseous Knudsen layer in Couette flow. It was shown that the velocity gradient singularity in the Knudsen layer can be explained within the continuum framework, when the nonlinearity of the constitutive model is morphed into the determination of the velocity slip in the second-order slip and jump model. Also, the smaller velocity slip and shear stress are shown to be caused by the shear-thinning property of the second-order constitutive model, that is, vanishing effective viscosity at high Knudsen number.


### 2.3. Knudsen Diffusion

Diffusive transport in micro- and nanoscale mediums can occur through three different mechanisms: (i) molecular bulk diffusion, (ii) surface diffusion (sorption), and (iii) Knudsen diffusion [67]. The proportion of each of this mechanism depends on the magnitude of the Knudsen number Kn, and Knudsen diffusion is an important mode of diffusion for large Kn such that the molecules interact more often with the walls of the medium than with each other.The molar flow rate of gas molecules in one direction of an axis (n˙±) is defined from ensemble averaging Einstein’s equation over a large enough number of molecules [68], and it is given as follows [69]:(17)n˙±=±16λAv¯∂c∂x.The entire molar flow rate (n˙) from the two opposite directions across an axis would be(18)n˙=16λAv¯∂c∂x−−16λAv¯∂c∂x=13λAv¯∂c∂x,where λ is the mean free path, A is the cross-sectional area, v¯ is the average velocity of gas molecules, and c is the gas concentration in mole per unit volume. Now, the diffusive flux of the gas would be(19)J=n˙A=13λv¯∂c∂x=−D∂c∂x,where D=1/3λv¯ is the coefficient of diffusion that varies with the kind of diffusion mechanism; for Knudsen diffusion, D would be equivalent to diffusion coefficient for Knudsen diffusion if λ is replaced by characteristic length of the medium s and v¯ follows Maxwellian distribution equal to 8RT/πM:(20)DKn=13s8RT/πM,where M is the molar mass of gas, R is the universal gas constant, and T is the temperature.The resistances of these three diffusive fluxes exist in series, so the effective diffusive flux as a result of these three mechanisms is the harmonic average [70, 71]:(21)1Jeff=1Jm+1JKn+1Jsorp,where Jeff is the effective flux that is the harmonic average of flux due to molecular diffusion Jm, flux due to Knudsen diffusion (JKn), and the flux due to sorption (Jsorp).


## 2.1. Fundamentals

### 2.1.1. Rarefaction and Slip Effect

Rarefaction effect in microsystems is a function of mean free path (MFP) of the gas. Usually the degree of rarefaction effect is characterized by the Knudsen number, which is defined as the ratio of the molecular mean free path (MFP) of a gas molecule to characteristic length of the medium:(1)Kn=λs,λ=MFP=kBT2πσ2p,where s is characteristic length of the medium, σ is diameter of a molecule, kB is the Boltzmann constant, T is absolute temperature, p is pressure, and λ is MFP.TheKn provides a direct indication whether the continuum approach for a flow of gas in a medium is applicable or not. Beyond the continuum approach, the flow is no longer near equilibrium; the assumptions of the no-slip boundary condition, thermodynamic equilibrium, and linear relation between the stress and the rate of strain fail. As the Kn increases, the rarefaction effect becomes more prominent which results in the following:(i)
Flow slippage at the wall(ii)
Flow is no long in thermodynamic equilibrium(iii)
Transfer of momentum from the bulk fluid region to the walls(iv)
Increase in mass flow rate(v)
The appropriate flow and heat transfer models depend on the range ofKn


### 2.1.2. Compressibility Effect

A change in pressure causes the gas properties to change in order to satisfy mass balance. For example, if the density of the gas decreases due to decrease in pressure, then it causes the gas to accelerate. Colin [9] recommended that compressibility effect is important for gas microflows when Mach number (Ma) is greater than 0.2, while Li and Guo [10] demonstrated experimentally that for an average Ma<0.3 the effect of compressibility is negligible. The compressibility and rarefaction phenomena in gas microflows are generally coupled together, and they can have an opposite effect on each other in the slip flow and early transition regimes based on the magnitude of Reynolds number (Re):(2)Kn=cπγ2MaRe,where γ is the ratio of specific heats of a gas. In the above expression, c=1 as given by Colin [9], whereas c=λ/s in Zhang et al. [11].Tang et al. [12] demonstrated experimentally that the compressibility effect causes the gas flow characteristics to deviate from conventional laminar flow at high Re regime (about Re=2000). In general, it has been demonstrated [13, 14] that the compressibility effect results in negative curvature of the pressure distribution, while the rarefaction effect weakens the negative curvature.


### 2.1.3. Flow Regimes

On the basis of the Knudsen number, different flow regimes with their physics and applicable models can be defined as shown in Table1. In the case of rarefied gas flow, continuum models are applicable for Kn≤0.01, and for Kn≥10 free-molecular models should be used. However, neither continuum models nor free-molecular models can be used for Kn lying in the intermediate range between the continuum and free-molecular flow regimes.Table 1
Different flow regimes as a function of Knudsen number.
Kn
Flow regime
Physics
Applicable models
Lower bound
Upper bound
0
10–2
Continuum (no slip flow)
Thermodynamic equilibrium and no-slip at the boundary
(1) Navier–Stokes equations, (2) Euler equations, (3) LBM
10–2
10–1
Slip flow
Nonequilibrium effects dominate near the walls. Assumptions of no-slip boundary condition, thermodynamic equilibrium, and linear stress-strain relationship fail
(1) Navier–Stokes equations with slip velocity and temperature jump boundary conditions, (2) Boltzmann gas-kinetic equation, (3) DSMC, (4) Burnett equations, (5) LBM, (6) gas-kinetic scheme, (7) method of moments
10–1
10
Transition flow
Rarefaction effects dominate and slip models become more complex. Assumptions of no-slip boundary condition, thermodynamic equilibrium, and linear stress-strain relationship fail
(1) Navier–Stokes equations with slip velocity and temperature jump boundary conditions, (2) Boltzmann gas-kinetic equation, (3) DSMC, (4) Burnett equations, (5) LBM, (6) gas-kinetic scheme, (7) method of moments
10
∞
Free-molecule flow
Collisions between gas molecules and boundary surface become dominant compared to intermolecular collisions. Assumptions of no-slip boundary condition, thermodynamic equilibrium, and linear stress-strain relationship fail
(1) Boltzmann gas-kinetic equation, (2) DSMC method
Flow regimes as a function of medium and flow parameters.Using the general definition of the Knudsen number that is defined as the ratio of the MFP of a gas molecule to characteristic length of the medium, and the mathematical formulation of MFP, we can writeKn as a function of the characteristic length of the flowing medium (s) and the average flowing pressure in that medium (p) as follows:(3)Kn=kBTs2πσ2p=fs,p.TheKn is estimated from the above equation for various characteristic lengths and average flowing pressures, and the corresponding flow regimes using the values of Kn are plotted as a function of characteristic lengths and flowing pressures as shown in Figure 1.Figure 1
(a) Variation ofKn values as a function of channel height and pressure, and (b) graphical representation of corresponding gas flow regimes. Here, Cont., Slip, Trans., and Mol. represent continuum flow regime, slip flow regime, transition flow regime, and molecular flow regime, respectively.
(a)
(b)Another measure ofKn was proposed by Tsien [15] that was based on thickness (δ) of the boundary layer on the flowing medium as Kn=λ/δ, where δ∝s/Re. Here Re is Reynold’s number for the characteristic length s and is defined as Re=ρus/μ, where ρ,u,andμ are fluid’s density, velocity, and viscosity, respectively. Using these relationships, Kn for the boundary layer thickness (δ) can be defined in terms of the freestream Mach number (M) and Reynold’s number as follows:(4)Knδ=λδ=MReδ∼MRe.TheKn is estimated from the above equation for various M and Re, and the corresponding flow regimes using the values of Kn are plotted as a function of M and Re as shown in Figure 2.Figure 2
(a) Variation ofMa/Re values as a function of Re and Ma, and (b) graphical representation of corresponding gas flow regimes. Here, Cont., Slip, Trans., and Mol. represent continuum flow regime, slip flow regime, transition flow regime, and molecular flow regime, respectively.
(a)
(b)It must be noted that Tsien’s parameter (Knδ) reduces to the ordinary Kn defined for the characteristic length in the free-molecular flow regime where Re≪1, such that δ≈s that results in Knδ=Kn.For flows with large amounts of nonequilibrium, modifying boundary conditions of continuum equations to account for slip will not eliminate all sources of errors. For example, the pressure tensor in more rarefied flows turns anisotropic, whereas the derivation of continuum equations assumes an isotropic pressure field [16, 17]. In those cases, it is recommended to use higher-order equations such as Burnett [18] equations and super-Burnett equations [19]. However, in some specific cases N-S equations with slight correction can be used to match data for a high Kn number regime (0.04–8.3) as shown by some authors [20, 21].


### 2.1.4. Accommodation Coefficients

The accommodation coefficients describe the macroscopic estimate of the fraction of momentum and energy transfer between the impinging/incident gas molecules and the boundary wall surface; therefore, their values vary between unity (complete accommodation of gas molecules with no reflection) and zero (complete reflection of gas molecules with no accommodation on the solid surface). These two accommodation coefficients for momentum and energy transfer [22] are key parameters in the estimation of slip velocity and temperature jump, respectively. A smaller value of accommodation coefficient would lead to higher velocity slip and temperature jump, because in that case the gas molecules are less affected by the solid surface.(1) Tangential Momentum Accommodation Coefficient. This coefficient defines the fraction of momentum transfer between the incident gas molecules and the boundary wall surface:(5)σu=τi−τrτi=1−τrτi,where σu is tangential momentum accommodation coefficient (TMAC), τi is the tangential momentum of the incident gas molecules, and τr is the tangential momentum of the reflected gas molecules. It can be seen that σu=1 if the reflected component of the tangential momentum of the gas molecules is zero (or fully diffusive), and σu=0 if the tangential momentum of the reflected molecules is equal to the tangential momentum of the incident molecules (or fully reflective/specular).Agrawal and Prabhu [23] presented the values of TMAC as a function of gas rarefaction magnitude (Kn) by using the data reported by various researchers in the literature for both monoatomic and diatomic gases. By fitting a curve to this data, they obtained and proposed the following expression for the TMAC:(6)σu=1−log1+Kn0.7.(2) Thermal Accommodation Coefficient. This coefficient defines the fraction of kinetic energy transfer between the incident gas molecules and the boundary wall surface. For a molecular beam type measurement, thermal accommodation coefficient is defined as follows:(7)σT=Ei−ErEi−Es,where σT is thermal accommodation coefficient (TAC), Ei and Er are the kinetic energies (for the normal component of velocity) of incoming and reflected gas molecules, respectively, and Es is the kinetic energy that would be produced by a diffusive reflection at the temperature of the solid surface Ts. Another definition of TAC that is directly related to measurement of temperatures is defined as follows:(8)σT=Ti−TrTi−Ts,where Ti, Tr, and Ts are the temperatures of the incident gas molecules, reflected gas molecules, and solid surface, respectively.Song and Yovanovich [24] presented the values of TAC as a function of gas molecular weight (Mg) and Ts by using the data reported by various researchers in the literature for both monoatomic, diatomic, and polyatomic gases. By fitting a curve to this data, they obtained and proposed the following expression for the TAC:(9)σT=expC0Ts−T0T0Mg∗C1+Mg∗+1−expC0Ts−T0T02.4μ1+μ2,where C0=−0.57 (dimensionless), C1=6.8g/mole, Ms is molecular weight of solid, μ=Mg/Ms, T0=273K, and Mg∗ is defined as follows:(10)Mg∗=Mg,for  monoatomic  gas,1.4Mg,for  dia/polyatomic  gases.


### 2.1.5. Knudsen Layer

As the flow deviates from the continuum regime, a sublayer with thickness of few mean free paths (MFPs) starts to develop from the wall of the medium such that it starts to become dominant between the bulk of the fluid and the wall of the medium (Figure3). The molecules within the KL collide more with the boundary wall than with other molecules [25], and this results in fluid (i) having nonzero velocity at the boundary surface and (ii) exhibiting non-Newtonian properties, thus introducing discontinuity in fluid properties between the bulk fluid and the KL.Figure 3
Schematic of the Knudsen layer near a boundary wall. The property in KL is discontinuous across which certain slip and jump conditions expressing conservation of momentum, mass, and energy are applied. The KL region is in nonequilibrium, whereas the bulk fluid flow region is in equilibrium.The flow in theKnudsen layer (KL) cannot be analyzed with conventional N-S equations, which does not account for certain slip and jump conditions expressing conservation of momentum, mass, and energy. The flow within the KL is typically regarded as being governed by the linearized Boltzmann gas-kinetic equation (LBE) [26–29] or by direct simulation Monte Carlo (DSMC) method [30–32]. Velocity slip results in increasing the mass flow rate by approximately 70%, and about one-third of this increase in flow rate is contributed by the non-Newtonian structure of the KL [33]. Therefore, accurate characterization of the KL is important in modeling of gas flow through microsystems.(1) Description of Flow in Knudsen Layer. Accurate simulation of rarefied flow effects requires accurate description and modeling of the KL. For this, it is important to know the performance of parameters at the nonequilibrium state within the KL. The flow within the KL can be described by three different approaches: (i) wall-function method, (ii) higher-order continuum method, and (iii) power-law method.The wall-function [33–35] is independent of the accommodation coefficient and is valid for planar surfaces up to Kn of 0.1. This function predicts the velocity gradient at the wall to be 1.7.Higher-order continuum models [33, 36] derived from the linear Boltzmann equation (LBE) [18, 19, 33] and truncated at disparate orders have been used to describe the flow in KL. One of the motivations for the development of these higher-order hydrodynamic models to simulate rarefied flows was to reduce the prohibitive computational demand of the DSMC method [37, 38]. However, there are questions regarding validity of these models [39] and in accurately capturing the KL [33, 40]. Additionally, these models cannot treat the boundary conditions at the wall [41, 42].Power-law model proposed by Lilley and Sader [40, 43] can be used to accurately describe the flow within the KL, and they suggest that this method describes a general physical phenomenon. The power-law model was derived by using the solution of LBE and DSMC, and it depicts the physics of velocity gradient not captured by the wall-function method [44] and the higher-order continuum method [33].(2) Thickness of Knudsen Layer. Lockerby et al. [33] presented the thicknesses of KL calculated from eight different models that used three different solution approaches of kinetic theory, higher-order continuum models, and molecular dynamics (MD) simulation. The thickness of KL obtained from these eight models varied between 0.9 λ and 4.9 λ. Gusarov and Smurov [45] presented the following expression to estimate an approximate value of KL thickness lc:(11)lc=kBTπd2p,where kB is the Boltzmann constant, T is the absolute temperature, and p is the gas pressure.


## 2.1.1. Rarefaction and Slip Effect

Rarefaction effect in microsystems is a function of mean free path (MFP) of the gas. Usually the degree of rarefaction effect is characterized by the Knudsen number, which is defined as the ratio of the molecular mean free path (MFP) of a gas molecule to characteristic length of the medium:(1)Kn=λs,λ=MFP=kBT2πσ2p,where s is characteristic length of the medium, σ is diameter of a molecule, kB is the Boltzmann constant, T is absolute temperature, p is pressure, and λ is MFP.TheKn provides a direct indication whether the continuum approach for a flow of gas in a medium is applicable or not. Beyond the continuum approach, the flow is no longer near equilibrium; the assumptions of the no-slip boundary condition, thermodynamic equilibrium, and linear relation between the stress and the rate of strain fail. As the Kn increases, the rarefaction effect becomes more prominent which results in the following:(i)
Flow slippage at the wall(ii)
Flow is no long in thermodynamic equilibrium(iii)
Transfer of momentum from the bulk fluid region to the walls(iv)
Increase in mass flow rate(v)
The appropriate flow and heat transfer models depend on the range ofKn


## 2.1.2. Compressibility Effect

A change in pressure causes the gas properties to change in order to satisfy mass balance. For example, if the density of the gas decreases due to decrease in pressure, then it causes the gas to accelerate. Colin [9] recommended that compressibility effect is important for gas microflows when Mach number (Ma) is greater than 0.2, while Li and Guo [10] demonstrated experimentally that for an average Ma<0.3 the effect of compressibility is negligible. The compressibility and rarefaction phenomena in gas microflows are generally coupled together, and they can have an opposite effect on each other in the slip flow and early transition regimes based on the magnitude of Reynolds number (Re):(2)Kn=cπγ2MaRe,where γ is the ratio of specific heats of a gas. In the above expression, c=1 as given by Colin [9], whereas c=λ/s in Zhang et al. [11].Tang et al. [12] demonstrated experimentally that the compressibility effect causes the gas flow characteristics to deviate from conventional laminar flow at high Re regime (about Re=2000). In general, it has been demonstrated [13, 14] that the compressibility effect results in negative curvature of the pressure distribution, while the rarefaction effect weakens the negative curvature.


## 2.1.3. Flow Regimes

On the basis of the Knudsen number, different flow regimes with their physics and applicable models can be defined as shown in Table1. In the case of rarefied gas flow, continuum models are applicable for Kn≤0.01, and for Kn≥10 free-molecular models should be used. However, neither continuum models nor free-molecular models can be used for Kn lying in the intermediate range between the continuum and free-molecular flow regimes.Table 1
Different flow regimes as a function of Knudsen number.
Kn
Flow regime
Physics
Applicable models
Lower bound
Upper bound
0
10–2
Continuum (no slip flow)
Thermodynamic equilibrium and no-slip at the boundary
(1) Navier–Stokes equations, (2) Euler equations, (3) LBM
10–2
10–1
Slip flow
Nonequilibrium effects dominate near the walls. Assumptions of no-slip boundary condition, thermodynamic equilibrium, and linear stress-strain relationship fail
(1) Navier–Stokes equations with slip velocity and temperature jump boundary conditions, (2) Boltzmann gas-kinetic equation, (3) DSMC, (4) Burnett equations, (5) LBM, (6) gas-kinetic scheme, (7) method of moments
10–1
10
Transition flow
Rarefaction effects dominate and slip models become more complex. Assumptions of no-slip boundary condition, thermodynamic equilibrium, and linear stress-strain relationship fail
(1) Navier–Stokes equations with slip velocity and temperature jump boundary conditions, (2) Boltzmann gas-kinetic equation, (3) DSMC, (4) Burnett equations, (5) LBM, (6) gas-kinetic scheme, (7) method of moments
10
∞
Free-molecule flow
Collisions between gas molecules and boundary surface become dominant compared to intermolecular collisions. Assumptions of no-slip boundary condition, thermodynamic equilibrium, and linear stress-strain relationship fail
(1) Boltzmann gas-kinetic equation, (2) DSMC method
Flow regimes as a function of medium and flow parameters.Using the general definition of the Knudsen number that is defined as the ratio of the MFP of a gas molecule to characteristic length of the medium, and the mathematical formulation of MFP, we can writeKn as a function of the characteristic length of the flowing medium (s) and the average flowing pressure in that medium (p) as follows:(3)Kn=kBTs2πσ2p=fs,p.TheKn is estimated from the above equation for various characteristic lengths and average flowing pressures, and the corresponding flow regimes using the values of Kn are plotted as a function of characteristic lengths and flowing pressures as shown in Figure 1.Figure 1
(a) Variation ofKn values as a function of channel height and pressure, and (b) graphical representation of corresponding gas flow regimes. Here, Cont., Slip, Trans., and Mol. represent continuum flow regime, slip flow regime, transition flow regime, and molecular flow regime, respectively.
(a)
(b)Another measure ofKn was proposed by Tsien [15] that was based on thickness (δ) of the boundary layer on the flowing medium as Kn=λ/δ, where δ∝s/Re. Here Re is Reynold’s number for the characteristic length s and is defined as Re=ρus/μ, where ρ,u,andμ are fluid’s density, velocity, and viscosity, respectively. Using these relationships, Kn for the boundary layer thickness (δ) can be defined in terms of the freestream Mach number (M) and Reynold’s number as follows:(4)Knδ=λδ=MReδ∼MRe.TheKn is estimated from the above equation for various M and Re, and the corresponding flow regimes using the values of Kn are plotted as a function of M and Re as shown in Figure 2.Figure 2
(a) Variation ofMa/Re values as a function of Re and Ma, and (b) graphical representation of corresponding gas flow regimes. Here, Cont., Slip, Trans., and Mol. represent continuum flow regime, slip flow regime, transition flow regime, and molecular flow regime, respectively.
(a)
(b)It must be noted that Tsien’s parameter (Knδ) reduces to the ordinary Kn defined for the characteristic length in the free-molecular flow regime where Re≪1, such that δ≈s that results in Knδ=Kn.For flows with large amounts of nonequilibrium, modifying boundary conditions of continuum equations to account for slip will not eliminate all sources of errors. For example, the pressure tensor in more rarefied flows turns anisotropic, whereas the derivation of continuum equations assumes an isotropic pressure field [16, 17]. In those cases, it is recommended to use higher-order equations such as Burnett [18] equations and super-Burnett equations [19]. However, in some specific cases N-S equations with slight correction can be used to match data for a high Kn number regime (0.04–8.3) as shown by some authors [20, 21].


## 2.1.4. Accommodation Coefficients

The accommodation coefficients describe the macroscopic estimate of the fraction of momentum and energy transfer between the impinging/incident gas molecules and the boundary wall surface; therefore, their values vary between unity (complete accommodation of gas molecules with no reflection) and zero (complete reflection of gas molecules with no accommodation on the solid surface). These two accommodation coefficients for momentum and energy transfer [22] are key parameters in the estimation of slip velocity and temperature jump, respectively. A smaller value of accommodation coefficient would lead to higher velocity slip and temperature jump, because in that case the gas molecules are less affected by the solid surface.(1) Tangential Momentum Accommodation Coefficient. This coefficient defines the fraction of momentum transfer between the incident gas molecules and the boundary wall surface:(5)σu=τi−τrτi=1−τrτi,where σu is tangential momentum accommodation coefficient (TMAC), τi is the tangential momentum of the incident gas molecules, and τr is the tangential momentum of the reflected gas molecules. It can be seen that σu=1 if the reflected component of the tangential momentum of the gas molecules is zero (or fully diffusive), and σu=0 if the tangential momentum of the reflected molecules is equal to the tangential momentum of the incident molecules (or fully reflective/specular).Agrawal and Prabhu [23] presented the values of TMAC as a function of gas rarefaction magnitude (Kn) by using the data reported by various researchers in the literature for both monoatomic and diatomic gases. By fitting a curve to this data, they obtained and proposed the following expression for the TMAC:(6)σu=1−log1+Kn0.7.(2) Thermal Accommodation Coefficient. This coefficient defines the fraction of kinetic energy transfer between the incident gas molecules and the boundary wall surface. For a molecular beam type measurement, thermal accommodation coefficient is defined as follows:(7)σT=Ei−ErEi−Es,where σT is thermal accommodation coefficient (TAC), Ei and Er are the kinetic energies (for the normal component of velocity) of incoming and reflected gas molecules, respectively, and Es is the kinetic energy that would be produced by a diffusive reflection at the temperature of the solid surface Ts. Another definition of TAC that is directly related to measurement of temperatures is defined as follows:(8)σT=Ti−TrTi−Ts,where Ti, Tr, and Ts are the temperatures of the incident gas molecules, reflected gas molecules, and solid surface, respectively.Song and Yovanovich [24] presented the values of TAC as a function of gas molecular weight (Mg) and Ts by using the data reported by various researchers in the literature for both monoatomic, diatomic, and polyatomic gases. By fitting a curve to this data, they obtained and proposed the following expression for the TAC:(9)σT=expC0Ts−T0T0Mg∗C1+Mg∗+1−expC0Ts−T0T02.4μ1+μ2,where C0=−0.57 (dimensionless), C1=6.8g/mole, Ms is molecular weight of solid, μ=Mg/Ms, T0=273K, and Mg∗ is defined as follows:(10)Mg∗=Mg,for  monoatomic  gas,1.4Mg,for  dia/polyatomic  gases.


## 2.1.5. Knudsen Layer

As the flow deviates from the continuum regime, a sublayer with thickness of few mean free paths (MFPs) starts to develop from the wall of the medium such that it starts to become dominant between the bulk of the fluid and the wall of the medium (Figure3). The molecules within the KL collide more with the boundary wall than with other molecules [25], and this results in fluid (i) having nonzero velocity at the boundary surface and (ii) exhibiting non-Newtonian properties, thus introducing discontinuity in fluid properties between the bulk fluid and the KL.Figure 3
Schematic of the Knudsen layer near a boundary wall. The property in KL is discontinuous across which certain slip and jump conditions expressing conservation of momentum, mass, and energy are applied. The KL region is in nonequilibrium, whereas the bulk fluid flow region is in equilibrium.The flow in theKnudsen layer (KL) cannot be analyzed with conventional N-S equations, which does not account for certain slip and jump conditions expressing conservation of momentum, mass, and energy. The flow within the KL is typically regarded as being governed by the linearized Boltzmann gas-kinetic equation (LBE) [26–29] or by direct simulation Monte Carlo (DSMC) method [30–32]. Velocity slip results in increasing the mass flow rate by approximately 70%, and about one-third of this increase in flow rate is contributed by the non-Newtonian structure of the KL [33]. Therefore, accurate characterization of the KL is important in modeling of gas flow through microsystems.(1) Description of Flow in Knudsen Layer. Accurate simulation of rarefied flow effects requires accurate description and modeling of the KL. For this, it is important to know the performance of parameters at the nonequilibrium state within the KL. The flow within the KL can be described by three different approaches: (i) wall-function method, (ii) higher-order continuum method, and (iii) power-law method.The wall-function [33–35] is independent of the accommodation coefficient and is valid for planar surfaces up to Kn of 0.1. This function predicts the velocity gradient at the wall to be 1.7.Higher-order continuum models [33, 36] derived from the linear Boltzmann equation (LBE) [18, 19, 33] and truncated at disparate orders have been used to describe the flow in KL. One of the motivations for the development of these higher-order hydrodynamic models to simulate rarefied flows was to reduce the prohibitive computational demand of the DSMC method [37, 38]. However, there are questions regarding validity of these models [39] and in accurately capturing the KL [33, 40]. Additionally, these models cannot treat the boundary conditions at the wall [41, 42].Power-law model proposed by Lilley and Sader [40, 43] can be used to accurately describe the flow within the KL, and they suggest that this method describes a general physical phenomenon. The power-law model was derived by using the solution of LBE and DSMC, and it depicts the physics of velocity gradient not captured by the wall-function method [44] and the higher-order continuum method [33].(2) Thickness of Knudsen Layer. Lockerby et al. [33] presented the thicknesses of KL calculated from eight different models that used three different solution approaches of kinetic theory, higher-order continuum models, and molecular dynamics (MD) simulation. The thickness of KL obtained from these eight models varied between 0.9 λ and 4.9 λ. Gusarov and Smurov [45] presented the following expression to estimate an approximate value of KL thickness lc:(11)lc=kBTπd2p,where kB is the Boltzmann constant, T is the absolute temperature, and p is the gas pressure.


## 2.2. Slip and Jump

Knudsen’s pioneering experimental and theoretical work [4, 5] on rarefied gas flows showed that the compressible form of the Navier–Stokes–Fourier model [6–8] is inadequate in describing rarefied fluid flows subject to the no-slip boundary condition. A common way to model the flow of rarefied gas is to modify the no-slip and no-jump boundary conditions in continuum models by revising the basic equations of momentum, mass, and heat transfer. There have been a number of approaches to account for slip boundary conditions through the Maxwell slip (1878), second-order slip [46], and Langmuir slip [47–49]. Maxwell boundary slip condition is the first-order approximation from kinetic theory of gases that has been widely used along with its derived forms in the past 135 years. One challenge in most gas flows at micro- and nanometer scale is the wide range of Knudsen numbers and the lack of a general slip model that can be applicable from the slip regime to the transition regime. Under this scenario, a different approach using molecular models that simulate the gas flow at the molecular level can be adopted. For example, one of most commonly used molecular simulation models is the Boltzmann gas-kinetic model [26–28], as well as many approximate methods based on particle and moments [50] like direct simulation Monte Carlo (DSMC) method and Burnett equations set [18, 26], respectively.Several studies on theory [16, 46, 51–53], Zhang et al. [11], and experimental data [54–57] of gaseous slip flow have cumulatively added to the knowledge base of gaseous slip. Below, we briefly review the fundamentals of velocity slip and temperature jump in terms of their governing equations and provide references for further discussions on these topics.


### 2.2.1. Velocity Slip

Although the concept of slip was first proposed by Navier in his model in which the magnitude of the fluid velocity is proportional to the magnitude of the flow shear rate, Maxwell [58] was the first person to quantify the slip length of a gas flowing over a solid boundary surface. According to Maxwell [58], for a dilute, ideal monoatomic gas, the first-order slip condition can be described as follows:(12)uslip=ut−uw=2−σuσuλ∂ut∂zw+34μρT∂T∂xw,where σu is TMAC, uslip is the gas slip velocity at the wall, ut is the gas velocity in the tangential direction, λ is the mean free path of the gas molecules, uw is the velocity of the solid wall, ρ and T are the density and temperature of the gas, respectively, and subscript w represents wall surface. The directional axes are represented as x,y,andz with respect to the direction of the flow. The second term on the right hand side represents the effect of thermal creep due to an axial temperature gradient, and it provides slip to the velocity at the wall in the presence of varying temperature.In case of a medium with curvature, the gas slip velocity at the wall can be expressed by adding an additional derivative of the gas velocity inz-direction as follows [59, 60]:(13)uslip=ut−uw=2−σuσuλ∂ut∂z+∂uz∂xw+34μρT∂T∂xw.Colin [61] suggested that for significantly rough surfaces, expression for slip velocity on curved surfaces is equally applicable. Other authors who studied the effect of roughness on rarified gas flow through microchannels [62, 63] made three important conclusions regarding the impact of surface roughness on the flow behavior of gas: (i) the flow behavior in the transition regime is more sensitive to the roughness height than that in the slip flow regime, (ii) influence of roughness height on the flow behavior is more significant than that of fractal dimension, and (iii) presence of surface roughness reduces the gas slip at the boundary wall irrespective of the flow regime. Pelević and van der Meer [64] studied the effect of surface roughness on heat transfer in microflows and reported that the roughness has a marginal effect in increasing heat transfer in microflows. For further readings on various approaches to model gaseous slip flow, readers can refer to the reviews by Cao et al. [16] and Zhang et al. [11].(1) Measurement of Gas Slip Velocity. Despite the presence of different correction coefficients [9] to account for slip flow at the boundary, there is no consensus on the most appropriate equation that can account for gas slip [65]. Additionally, accurate values of TMAC are generally not well known in the literature as it is a function of type of gas, material, and degree of rarefaction. These limitations are usually accounted by experimental data on slip length of rarefied gas flows that is defined by Maxwell’s equation in its simplest form for an isothermal surface as follows:(14)uslip=2−σuσuλ∂ut∂zw=b∂ut∂zw,where b=2−σu/σuλ, which is called as the slip length of rarefied gas flow. Maali et al. [65] reviewed the experimental techniques used to measure the slip length and TMAC for gas flow on isothermal surface. These techniques are Millikan experiment, rotating body techniques, mass flow rate measurements, and atomic force microscopy method, respectively. The gas slip is measured indirectly from these experimental techniques by using slip length, which is deduced from the measurement of the velocity of the drop in the Millikan method, the deflection of the cylinder in rotating body technique, the mass flow rate in rarefied gas flow through microchannels, and the force in atomic force microscopy method, respectively. For detailed discussion on these experimental techniques and their use in measuring gas slip, readers can refer to the review by Maali et al. [65].


### 2.2.2. Temperature Jump

Similar to the velocity slip expressions presented above, temperature jump boundary conditions can be written as follows:(15)Tjump=Tg−Tw=2−σTσT2γγ+1λPr∂T∂zw,where σT is thermal accommodation coefficient, Tjump is the temperature jump at the wall, Tg is the gas temperature, Tw is the temperature of the solid wall, γ is specific heat ratio, and Pr is the Prandtl number. For further discussions on temperature jump, readers can refer to the reviews by Sharipov [56] and Colin [61].


### 2.2.3. Second-Order Velocity Slip and Temperature Jump

The linear velocity slip equation (12) and temperature jump equation (15) were obtained from the original (nonlinear) models in which the velocity slip and temperature jump are assumed directly proportional to the degree of nonequilibrium near the wall surface [60]:(16)uslip=ut−uw=2−σuσuλΠμw+34γ−1γ/PrQxpw,Tjump=Tg−Tw=2−σTσT2γγ+1λPrQzkw.Only after the degree of nonequilibrium is taken as linear with the first-order accuracy (that is, shear stressΠ∼∂ut/∂z, tangential heat flux Qx∼∂T/∂x, and normal heat flux Qz∼∂T/∂z), the models are reduced to the well-known first-order models (Equations (12) and (15)). Therefore, if the linear constitutive relations are replaced by higher-order relations like the second-order constitutive relations, the resulting models become a second-order velocity slip and temperature jump. Such models were recently developed by Myong [66] and were applied to the theoretical analysis of the gaseous Knudsen layer in Couette flow. It was shown that the velocity gradient singularity in the Knudsen layer can be explained within the continuum framework, when the nonlinearity of the constitutive model is morphed into the determination of the velocity slip in the second-order slip and jump model. Also, the smaller velocity slip and shear stress are shown to be caused by the shear-thinning property of the second-order constitutive model, that is, vanishing effective viscosity at high Knudsen number.


## 2.2.1. Velocity Slip

Although the concept of slip was first proposed by Navier in his model in which the magnitude of the fluid velocity is proportional to the magnitude of the flow shear rate, Maxwell [58] was the first person to quantify the slip length of a gas flowing over a solid boundary surface. According to Maxwell [58], for a dilute, ideal monoatomic gas, the first-order slip condition can be described as follows:(12)uslip=ut−uw=2−σuσuλ∂ut∂zw+34μρT∂T∂xw,where σu is TMAC, uslip is the gas slip velocity at the wall, ut is the gas velocity in the tangential direction, λ is the mean free path of the gas molecules, uw is the velocity of the solid wall, ρ and T are the density and temperature of the gas, respectively, and subscript w represents wall surface. The directional axes are represented as x,y,andz with respect to the direction of the flow. The second term on the right hand side represents the effect of thermal creep due to an axial temperature gradient, and it provides slip to the velocity at the wall in the presence of varying temperature.In case of a medium with curvature, the gas slip velocity at the wall can be expressed by adding an additional derivative of the gas velocity inz-direction as follows [59, 60]:(13)uslip=ut−uw=2−σuσuλ∂ut∂z+∂uz∂xw+34μρT∂T∂xw.Colin [61] suggested that for significantly rough surfaces, expression for slip velocity on curved surfaces is equally applicable. Other authors who studied the effect of roughness on rarified gas flow through microchannels [62, 63] made three important conclusions regarding the impact of surface roughness on the flow behavior of gas: (i) the flow behavior in the transition regime is more sensitive to the roughness height than that in the slip flow regime, (ii) influence of roughness height on the flow behavior is more significant than that of fractal dimension, and (iii) presence of surface roughness reduces the gas slip at the boundary wall irrespective of the flow regime. Pelević and van der Meer [64] studied the effect of surface roughness on heat transfer in microflows and reported that the roughness has a marginal effect in increasing heat transfer in microflows. For further readings on various approaches to model gaseous slip flow, readers can refer to the reviews by Cao et al. [16] and Zhang et al. [11].(1) Measurement of Gas Slip Velocity. Despite the presence of different correction coefficients [9] to account for slip flow at the boundary, there is no consensus on the most appropriate equation that can account for gas slip [65]. Additionally, accurate values of TMAC are generally not well known in the literature as it is a function of type of gas, material, and degree of rarefaction. These limitations are usually accounted by experimental data on slip length of rarefied gas flows that is defined by Maxwell’s equation in its simplest form for an isothermal surface as follows:(14)uslip=2−σuσuλ∂ut∂zw=b∂ut∂zw,where b=2−σu/σuλ, which is called as the slip length of rarefied gas flow. Maali et al. [65] reviewed the experimental techniques used to measure the slip length and TMAC for gas flow on isothermal surface. These techniques are Millikan experiment, rotating body techniques, mass flow rate measurements, and atomic force microscopy method, respectively. The gas slip is measured indirectly from these experimental techniques by using slip length, which is deduced from the measurement of the velocity of the drop in the Millikan method, the deflection of the cylinder in rotating body technique, the mass flow rate in rarefied gas flow through microchannels, and the force in atomic force microscopy method, respectively. For detailed discussion on these experimental techniques and their use in measuring gas slip, readers can refer to the review by Maali et al. [65].


## 2.2.2. Temperature Jump

Similar to the velocity slip expressions presented above, temperature jump boundary conditions can be written as follows:(15)Tjump=Tg−Tw=2−σTσT2γγ+1λPr∂T∂zw,where σT is thermal accommodation coefficient, Tjump is the temperature jump at the wall, Tg is the gas temperature, Tw is the temperature of the solid wall, γ is specific heat ratio, and Pr is the Prandtl number. For further discussions on temperature jump, readers can refer to the reviews by Sharipov [56] and Colin [61].


## 2.2.3. Second-Order Velocity Slip and Temperature Jump

The linear velocity slip equation (12) and temperature jump equation (15) were obtained from the original (nonlinear) models in which the velocity slip and temperature jump are assumed directly proportional to the degree of nonequilibrium near the wall surface [60]:(16)uslip=ut−uw=2−σuσuλΠμw+34γ−1γ/PrQxpw,Tjump=Tg−Tw=2−σTσT2γγ+1λPrQzkw.Only after the degree of nonequilibrium is taken as linear with the first-order accuracy (that is, shear stressΠ∼∂ut/∂z, tangential heat flux Qx∼∂T/∂x, and normal heat flux Qz∼∂T/∂z), the models are reduced to the well-known first-order models (Equations (12) and (15)). Therefore, if the linear constitutive relations are replaced by higher-order relations like the second-order constitutive relations, the resulting models become a second-order velocity slip and temperature jump. Such models were recently developed by Myong [66] and were applied to the theoretical analysis of the gaseous Knudsen layer in Couette flow. It was shown that the velocity gradient singularity in the Knudsen layer can be explained within the continuum framework, when the nonlinearity of the constitutive model is morphed into the determination of the velocity slip in the second-order slip and jump model. Also, the smaller velocity slip and shear stress are shown to be caused by the shear-thinning property of the second-order constitutive model, that is, vanishing effective viscosity at high Knudsen number.


## 2.3. Knudsen Diffusion

Diffusive transport in micro- and nanoscale mediums can occur through three different mechanisms: (i) molecular bulk diffusion, (ii) surface diffusion (sorption), and (iii) Knudsen diffusion [67]. The proportion of each of this mechanism depends on the magnitude of the Knudsen number Kn, and Knudsen diffusion is an important mode of diffusion for large Kn such that the molecules interact more often with the walls of the medium than with each other.The molar flow rate of gas molecules in one direction of an axis (n˙±) is defined from ensemble averaging Einstein’s equation over a large enough number of molecules [68], and it is given as follows [69]:(17)n˙±=±16λAv¯∂c∂x.The entire molar flow rate (n˙) from the two opposite directions across an axis would be(18)n˙=16λAv¯∂c∂x−−16λAv¯∂c∂x=13λAv¯∂c∂x,where λ is the mean free path, A is the cross-sectional area, v¯ is the average velocity of gas molecules, and c is the gas concentration in mole per unit volume. Now, the diffusive flux of the gas would be(19)J=n˙A=13λv¯∂c∂x=−D∂c∂x,where D=1/3λv¯ is the coefficient of diffusion that varies with the kind of diffusion mechanism; for Knudsen diffusion, D would be equivalent to diffusion coefficient for Knudsen diffusion if λ is replaced by characteristic length of the medium s and v¯ follows Maxwellian distribution equal to 8RT/πM:(20)DKn=13s8RT/πM,where M is the molar mass of gas, R is the universal gas constant, and T is the temperature.The resistances of these three diffusive fluxes exist in series, so the effective diffusive flux as a result of these three mechanisms is the harmonic average [70, 71]:(21)1Jeff=1Jm+1JKn+1Jsorp,where Jeff is the effective flux that is the harmonic average of flux due to molecular diffusion Jm, flux due to Knudsen diffusion (JKn), and the flux due to sorption (Jsorp).


## 3. Liquid Flow

### 3.1. Fundamentals

Unlike the well-developed kinetic theory for gases, no such molecular level transport theory exists for liquids because of which it is difficult to predict molecular effects in liquids. Even though a detailed molecular theory for the thermal conductivity of monatomic liquids was developed more than half a century ago [72], the means to implement it for practical calculations have not been developed. As a result of that gap, some rough theories and empirical estimation methods are used instead [73], and one that is generally used quite often was proposed by Bridgman [74] to predict the energy transport in pure liquids. Bridgman [74] assumed the molecules of a liquid are arranged in cubic lattice, with a center-to-center spacing δ given by(22)δ=V¯NA1/3,where V¯ is the molar volume and NA is Avogadro’s number. He further hypothesized that the energy transfers from one lattice plane to the next plane with a sonic velocity vs for a given liquid, and by rearranging the thermal conductivity from the rigid-sphere gas theory he derived the following equation that is known as Bridgman’s equation:(23)k=3NAV¯2/3κvs.The above equation shows good agreement with experimental data even for polyatomic liquids, but later it was suggested that a slightly lower coefficient value of 2.80 instead of 3 gives a better agreement with the data. Here, the velocity of low-frequency sound is given as follows:(24)vs=CpCV∂p∂ρT,where Cp/CV is the ratio of heat capacities that is nearly unity for liquids (except near the critical point) and ∂p/∂ρT can be determined from isothermal compressibility measurements of change in density ρ with the corresponding change in pressure (p) of the liquid.


#### 3.1.1. Flow Regimes

The flow regimes in gases are based on the concept ofKn, which is the ratio of the mean free path of a molecule to the characteristic length of the medium. The spacing between molecules in liquids is much smaller than the spacing between molecules in gases, so the concept of mean free path and the flow regimes based on Kn are not valid in liquids. Because of the lack of a well-developed molecular theory for liquids, no flow regimes have been proposed for liquids in the literature. However, there is one study by Liu and Li [75] that talks about the liquid flow regimes in nanochannels. Liu and Li [75] investigated the liquid motion as a function of four parameters (fluid-fluid binding energy εff, fluid-solid wall binding energy εfw, temperature of the system T, and external driving force Fe corresponding to the pressure drop along the channel) by using nonequilibrium molecular dynamics simulations (NEMDs). They define two dimensionless numbers to understand the competing effect between the liquid-liquid interaction and liquid-surface interaction as εff/kBT and εfw/kBT×εfw/Feσfw, respectively, where kB is the Boltzmann constant and σfw is the collision diameter (separation) at which the potential vanishes. They compute the flux using liquid argon and liquid helium flowing through a channel made of silver walls at temperatures varying from 100 K to 300 K and pressure difference varying from 33 MPa to 132 MPa. They find liquid flux values as a function of these two dimensionless numbers and illustrate their results by dividing them into different regimes, each associated with a distinct mechanism. The authors noted that the relation of the four parameters was not clear if the slip length was used instead of flux. The reason why slip length may not directly relate to the flow is because the density and viscosity of the confined fluid are nonuniform.Using the data from the study of Liu and Li [75], we graphically illustrate the flow regimes in liquid in Figure 4 and summarize them theoretically in Table 2.Figure 4
(a) Variation of flux values as a function of fluid-wall interaction and fluid-fluid interaction, and (b) graphical representation of corresponding liquid flow regimes.
(a)
(b)Table 2
Different flow regimes as a function of liquid/solid binding force.
Interaction between
Flux
Flow regime
Physics
Liquid-liquid
Liquid-solid
Low
Low
High
Flux independent of solid-fluid binding force and increases with decreasing temperature
Surface effects are weak with no liquid adsorption. Flux is dominated by the temperature-driven friction between the solid and the liquid
Low
High
Intermediate
Flux independent of temperature and dependent on solid-fluid binding force
Surface effects become strong and liquid adsorption starts to take place that results in increasing fluid density close to the solid surface
High
Low
Intermediate
Flux independent of solid-fluid interaction and increases with decreasing temperature
Surface effects are weak with no liquid adsorption. Flux is dominated by the temperature-driven friction between the solid and the liquid
High
High
Low
Flux independent of solid-fluid interaction and increases with increasing temperature
Surface effects become strong and liquid adsorption starts to take place that results in increasing fluid density close to the solid surface. Flux is enhanced with increasing temperature as it results in desorption of the fluid molecules that merge together with the bulk liquid


### 3.2. Liquid Slip

Fundamental understanding of the mechanisms that control liquid slip is still elusive. Unlike gases, the concept of mean free path is not useful for liquids because the average distance between liquid molecules is comparable to the molecular diameter. Even though it is believed unanimously that the flow enhancement in nanoscale systems is caused by liquid slip velocity at the solid surface [76, 77], the debate on understanding its mechanisms mostly revolves around four key parameters. Many theories have been developed to associate liquid slip length (Ls) as a function of liquid properties [78, 79], interaction strength between the fluid and wall material [80–84], wettability of the medium [77, 85–89], surface roughness [90–95], and liquid shear rate [81, 86, 96–100].


#### 3.2.1. Physical Mechanisms of Liquid Slip

The physical mechanisms behind the slip are still elusive, though there are few studies [101–107] that conjecture the reasons behind the liquid slip as due to the presence of either (i) a depleted water layer or (ii) an effective air gap at the wall formed by the nanobubbles. Another approach to hypothesize liquid slip is based on intermolecular forces affecting the momentum transport from the liquid to the wall [16]. Barring two studies [81, 108], there are no molecular-based theories to characterize this process [109, 110]. Lichter et al. [81] and Martini et al. [108] attribute the slip mechanism to hopping of liquid atoms from one equilibrium site to another, which follows Arrhenius dynamics based on the rate of hopping.


#### 3.2.2. Slip Length

The slip velocity of the liquid at the surface boundary is characterized by an averaged property of all the liquid molecules called slip lengthLs, which is defined as the linearly extrapolated perpendicular distance from the surface boundary where fluid velocity becomes zero. In other words, slip length is the perpendicular distance from the solid boundary, outside the region of liquid flow, where the no-flow boundary condition is achieved. Figure 5 shows a schematic to illustrate the concept of slip for liquid flow; this figure shows three types of flow from left to right that indicate flow with no slip, flow with partial slip, and flow with full slip, respectively.Figure 5
Schematic to illustrate liquid slip on a stationary solid surface.Navier first introduced the linear boundary condition that is considered to be the standard and most basic characterization of liquid slip. Theoretically, liquid slip velocity (vs) is defined as follows:(25)vs=Lsn⋅dvxzdz+dvxzdzT⋅1−nn.Pressure gradient-driven fluid velocity (vx) in the presence of slip can then be theoretically calculated as follows [111]:(26)vx=h212η−dPdx1+6Lsh,where Ls is the slip length, n denotes the normal to the surface along the z direction, vs is the slip velocity at the surface, vx is the liquid velocity in a channel, η is liquid shear viscosity, dP/dx is the pressure gradient along x-direction, and h is the height of the channel.Research on liquid slip started with a focus on slip over smooth and hydrophobic surfaces [79, 105, 112, 113], but there has been an increasing evidence of liquid slip at different types of surface, including rough surface [90, 114] and hydrophilic surface [88, 115, 116]. Recent progress in liquid slip has broadened to various other parameters that have shown to influence liquid slip, and as a result, there is a continuous effort in developing appropriate models to characterize the slip on different types of surfaces and under different conditions. Below, we review the progress of liquid slip models as a function of (i) surface wettability, (ii) surface property, (iii) shear rate, and (iv) surface roughness.(1) Effect of Surface Wettability. The liquid slip at nanoscale is greatly affected by the surface wettability of the medium. Macroscopically, the wettability of a surface is classified as wetting, nonwetting, or neutral depending on the contact angle formed by a liquid droplet with the surface. A surface forming contact angles larger than 90° with a liquid is termed as nonwetting, whereas a surface forming contact angles smaller than 90° with a liquid is termed as wetting fluid. A surface forming contact angle equal to 90° with a liquid is termed as neutral. Nonwetting surfaces are known to have a high disposition towards introducing liquid slip [77, 86, 87, 117]. For example, hydrodynamic slip length of simple liquids flowing on smooth nonwetting surfaces (also referred as hydrophobic surfaces in the context of water as a liquid) is typically 20 nm in mediums with size varying from one to several hundreds of nanometers [112]. On the other hand, disposition of wetting surfaces in introducing slip is not as straight forward as it is for the nonwetting surfaces. There are reports [118] that suggest wetting surfaces are neutral in their disposition to introduce liquid slip, whereas others [88] have demonstrated that the wettability factor alone is not necessary to attain liquid slip and that it should be combined with the strength of the solid-fluid interactions. For example, if there is favorable sorption sites that are separated from each other by well-defined sub-nanometer distances, no slip takes place. However, liquid slip can occur if those favorable sorption sites are close to each other, provided liquid-solid attractions are not too strong. That is because nonwetting surfaces such as graphite are usually characterized by uniform distributions of fluid molecules, while the wetting surfaces such as crystalline silica are characterized by discontinuous sorption sites with high affinity for adsorbing fluid molecules. Voronov et al. [89] found that the slip length grows if the surface’s wettability is turned into more wetting by modifying the fluid-wall interfacial energy through affinity strength parameters. Gao and McCarthy [119] demonstrated experimentally that static contact angle of a fluid on a surface should not be the sole criteria to characterize the nature of wettability in order to determine the tendency of a liquid to slip. Instead, they suggested that contact angle hysteresis should be used to determine the nonwetting vs. wetting behavior of a surface.It is believed that Tolstoi [87] was the first person to propose the liquid slip by linking the macroscopic thermodynamics concept at the molecular scale. He quantified the mobility of liquid molecules with the surface energies (=σ2γ) through spreading coefficient (=γ1−cosθ) as follows:(27)Lsσ∼expασ2γ1−cosθkBT−1,where α is a dimensionless geometrical parameter of first order that represents the microcavity area within the solid, σ is molecular diameter, γ is the surface tension between the liquid and the surface, θ is the equilibrium/static contact angle, T is the absolute temperature, kB is the Boltzmann constant. Another theory [85, 120, 121] that describes liquid slip at the molecular level with thermodynamic equilibrium uses the Green-Kubo relation and fluctuation-dissipation theorem (linear regression of fluctuations) and is given by the following equation:(28)Lsσ∼D∗StcLS2ρcσ3,where D∗=Dp/D0, Dp is the collective molecular diffusion coefficient, D0 is the bulk diffusion coefficient, St is a factor that represents first molecular layer structure, ρc is the fluid density at the first molecular layer, σ is the molecular diameter, and cLS is the dimensionless coefficient for solid-liquid Lennard-Jones potential.Huang et al. [122] and Sendner et al. [116] studied the slip length of shear-driven water on various hydrophobic surfaces using nonequilibrium molecular dynamic simulations (NEMDs). They show that the slip length as a function of different surfaces (both organic-like silane monolayers and inorganic-like Lennard-Jones models) can be characterized using the static contact angle (θ) of the surface as shown below, which they describe as the crucial parameter in controlling the water slippage:(29)Ls∝1+cosθ−2.The authors demonstrate that their formula can also be obtained from the linear response theory by using the LJ fluid-solid energy parameter and the Young-Laplace equation. Their estimated slippage length varies from a few nanometers to tens of nanometers. They also proposed the following scaling expression for slip length as a function of liquid depletion length, which they define as the layer between the water and the solid surface occupied by the gas:(30)Ls∝δ4,where δ is the thickness of liquid depletion and given by(31)δ=∫0∞1−ρszρsb−ρlzρlbdz,where ρsz and ρlz are the densities of the solid and liquid along the height of the channel, whereas ρsb and ρlb are the bulk densities of the solid and liquid, respectively.Bakli and Chakraborty [123] presented a constitutive model to describe the nanopore imbibition characteristics of water as a function of the dynamic slip length using NEMDs. Their derived slip length depends on the contact angle and an acceleration parameter that comes from the pressure gradients due to capillary forces. They show that for low acceleration values, their model approaches the model given by Huang et al. [122]:(32)Ls=A1a1+cosθ2−A2a,a=γcosθρRL,a∗=γcosθηLdL/dt,arough∗=asmooth∗+a¯,a¯=0.35exp−0.6N,where L is the capillary imbibition length at time t, N is the parameter representing the surface roughness (lower value represents higher roughness), a is the imbibition driving acceleration parameter, a∗ is the dimensionless form of a, a¯ is a deceleration parameter to mimic flow hindrance due to surface roughness, and A1a and A2a are functions of a.(2) Effect of Shear Rate. The effect of shear rate on slip length has been investigated experimentally by many researchers [90, 96, 97, 99, 100] and observed to depend on the shear rate. In such a scenario where slip length is found to depend on the shear rate, the slip boundary condition proposed by Navier, and presented earlier, shows a nonlinear behavior with shear rate. Although several slip length relationships have been proposed as a function of shear rate, the relationship proposed by Priezjev [98] is an important one as it accounts for the solid-fluid interaction strength that is crucial in micro- or nanoporous media. Priezjev [98] proposed slip length for smooth surfaces as a function of shear rate that predicts that the slip length increases nonlinearly with the shear rate for weak solid-fluid interactions provided that the solid-liquid interface forms a discrete structure, whereas if the solid-fluid interaction increases, the slip length becomes linearly proportional to the shear rate:(33)Qslip=23ρfxh3μ+2hVs.The above equation can be normalized in terms of the no-slip flow rate as follows:(34)Qslip=Qno−slip1+6Ls2h,where fx is an externally applied force in the x direction, ρ is the fluid density, h is the distance between the confining solid walls, μ is the fluid viscosity, Vs is the fluid velocity at the confining parallel walls that can also be written as Vs=v−h=vh, and Qslip and Qno−slip are the flow rates due to slip and no-slip boundary conditions, respectively. The term 2hVs in the above equation represents the correction to the flow rate due to slip boundary conditions. Priezjev [98] reported that the fluid viscosity, shown below, was found to be shear rate (γ˙) independent in the range of γ˙≤0.16τ−1, where τ is the characteristic Lennard-Jones (LJ) time:(35)μ=2.2±0.2ετσ−3,where ε and σ represent the energy and length scales of the fluid phase.(3) Effect of Surface Roughness. Even though most of the solids are naturally rough at micro- to nanoscale, slip length of liquid as a function of surface roughness is not well studied [16, 124]. The problem of calculating slip as a function of surface roughness is challenging [16, 125, 126] because of four reasons: (i) controlling surface roughness in a natural and a realistic manner at nanometer scale is not easy, (ii) surface roughness can introduce changes to some interfacial properties such as interfacial forces, wettability, and trapped gases, which may be undesirable, (iii) rough surface will introduce undulating boundaries that will introduce complexities in making correct interpretations with respect to varying boundary positions, and (iv) there is no analytical characterization for the nonperiodic and natural randomness of surface roughness.Although several slip length relationships have been proposed as a function of surface roughness, two relationships by Lund et al. [127] and Misra and Bakli [128] are discussed here because of their simplicity and effectiveness.Lund et al. [127] derived an expression of effective slip length on a surface with periodic roughness in terms of intrinsic slip length and contact area of the surface, using steady-state Stokes flow and homogenization method. Their derived effective slip length represents the harmonic mean weighted by the area of contact between surface and fluid as follows:(36)Leff=1+s2Ls−1,where Ls is the intrinsic slip length, s is the slope, and 1+s2=1/L∫0L1+s2dx is the arc length of the surface function over one period, normalized by the length of the period (L). For a flat surface (s=0), the above expression reduces to the intrinsic slip length. Their result is applicable for rough surfaces where L≪Lsmin.Misra and Bakli [128] recently proposed an expression of effective slip length for slippage of water considering coupled effect of the surface chemistry and surface roughness by performing molecular dynamics simulations. They found that the slip length without surface roughness match closely with the result of 1+cosθ2 proposed by Huang et al. [122] for smooth surface, while the presence of surface roughness always reduces the effective slip length as follows:(37)Leff=k1+cosθ21+α∗2,α∗=∑fx,y−fx,y¯N,where k is the constant of proportionality and α∗ is a parameter that characterizes the roughness of a surface through its amplitude N and roughness function fx,y that modifies the attractive part of the Lennard-Jones (LJ) potential.(4) Electrokinetic Effects. The amount of liquid slip on a surface is found to be affected by the presence of a charged surface and/or liquids with charged ions such as electrolytes. The presence of charge leads to attraction of oppositely charged ions and repulsion of similarly charged ions, and in the presence of a charged surface, this leads to the formation of a charged interfacial layer near the surface. The fluid flow in the presence of charge (also referred to as electrokinetic phenomenon) is centrally connected to the key concept of electrical double layer (EDL) that is made up of a stern layer and a diffuse layer as illustrated in Figure 6. Stern layer is the width of a layer close to the surface where the underlying solid surface charge (cations or anions) attracts the opposite charge of the solvent/liquid to form a layer of highly concentrated ions, whereas the diffuse layer is the width of a layer starting from the stern layer and ending where the electrokinetic potential goes to zero.Figure 6
Schematic of a flow in a channel with electrokinetic effects.The electrokinetic effects result in enhanced solid-fluid interaction close to the surface, which may also influence the liquid slip flow and consequently the overall fluid transport. Two examples of electrokinetic flows relevant here are electro-osmosis (liquid flow induced by application of an electric field) and streaming currents (electric current induced by application of a pressure gradient).The surface-fluid interaction due to electrokinetic effects is more complex than the surface-fluid interaction without any charge. Research on this area is sparse due to limited applications that have surface charge present. Out of the few models that exist in the literature, models proposed by Joly et al. [129] and Choudhary et al. [130] are discussed here because of parameters that are directly related to the electrokinetic effects and solid-fluid interaction compared to other models where the parameters are not that apparent.Joly et al. [129] explored the effect of surface charge to solid-liquid interfacial friction and its effect on the slip length boundary condition. They study the hydrodynamics of liquid coupled with electrostatics within the EDL using molecular dynamics simulations. Their study found that the presence of surface charge on the solid-liquid interface enhances the frictional force between the solid-liquid and as a result decreases the slip length.(38)Leff=Ls1+1/ασd2/e2lB/d2Ls,where Ls is the slip length without considering surface charge, α≈1 is a numerical factor, σ is the surface charge density, d is the equilibrium distance of Lennard-Jones potential, e is the elementary charge, and lB (∼0.7 nm in water at room temperature) is the Bjerrum length that represents a length scale below which electrostatic interactions dominate thermal effects.Choudhary et al. [130] studied the liquid slip for electro-osmotic flow through a nanochannel with hydrophobic walls of sinusoidally varying slippage using an asymptotic theory that uses the ratio of pattern amplitude (perpendicular to the applied electric field) to the average slip as a parameter. They proposed the following analytical expression for effective slip that can be used to design slip through engineered variations in topography and/or chemistry:(39)L¯eff=L¯s−α2L¯s2k¯s2sinhk¯s1+cosθ2coshk¯s+L¯s.k¯ssinhk¯s+coshk¯s1−cosθ2sinhk¯s+L¯s.k¯scoshk¯s,where L¯eff is the dimensionless effective slip length due to electrokinetic effects and L¯s is the dimensionless average slip length (=Lsx/h) that is equal to Ls01+αcosk¯sx¯+θ and Ls01+αcosk¯sx¯ on the top and bottom walls of the channel with height 2h, respectively. Here, Ls0 is the average slip length without electrokinetic effects, α is the dimensionless amplitude (=a/Ls0) of the pattern on either walls charged to facilitate electro-osmosis flow, k¯s is the wavenumber, and θ is the phase angle between the slip waves on the two walls.


## 3.1. Fundamentals

Unlike the well-developed kinetic theory for gases, no such molecular level transport theory exists for liquids because of which it is difficult to predict molecular effects in liquids. Even though a detailed molecular theory for the thermal conductivity of monatomic liquids was developed more than half a century ago [72], the means to implement it for practical calculations have not been developed. As a result of that gap, some rough theories and empirical estimation methods are used instead [73], and one that is generally used quite often was proposed by Bridgman [74] to predict the energy transport in pure liquids. Bridgman [74] assumed the molecules of a liquid are arranged in cubic lattice, with a center-to-center spacing δ given by(22)δ=V¯NA1/3,where V¯ is the molar volume and NA is Avogadro’s number. He further hypothesized that the energy transfers from one lattice plane to the next plane with a sonic velocity vs for a given liquid, and by rearranging the thermal conductivity from the rigid-sphere gas theory he derived the following equation that is known as Bridgman’s equation:(23)k=3NAV¯2/3κvs.The above equation shows good agreement with experimental data even for polyatomic liquids, but later it was suggested that a slightly lower coefficient value of 2.80 instead of 3 gives a better agreement with the data. Here, the velocity of low-frequency sound is given as follows:(24)vs=CpCV∂p∂ρT,where Cp/CV is the ratio of heat capacities that is nearly unity for liquids (except near the critical point) and ∂p/∂ρT can be determined from isothermal compressibility measurements of change in density ρ with the corresponding change in pressure (p) of the liquid.


### 3.1.1. Flow Regimes

The flow regimes in gases are based on the concept ofKn, which is the ratio of the mean free path of a molecule to the characteristic length of the medium. The spacing between molecules in liquids is much smaller than the spacing between molecules in gases, so the concept of mean free path and the flow regimes based on Kn are not valid in liquids. Because of the lack of a well-developed molecular theory for liquids, no flow regimes have been proposed for liquids in the literature. However, there is one study by Liu and Li [75] that talks about the liquid flow regimes in nanochannels. Liu and Li [75] investigated the liquid motion as a function of four parameters (fluid-fluid binding energy εff, fluid-solid wall binding energy εfw, temperature of the system T, and external driving force Fe corresponding to the pressure drop along the channel) by using nonequilibrium molecular dynamics simulations (NEMDs). They define two dimensionless numbers to understand the competing effect between the liquid-liquid interaction and liquid-surface interaction as εff/kBT and εfw/kBT×εfw/Feσfw, respectively, where kB is the Boltzmann constant and σfw is the collision diameter (separation) at which the potential vanishes. They compute the flux using liquid argon and liquid helium flowing through a channel made of silver walls at temperatures varying from 100 K to 300 K and pressure difference varying from 33 MPa to 132 MPa. They find liquid flux values as a function of these two dimensionless numbers and illustrate their results by dividing them into different regimes, each associated with a distinct mechanism. The authors noted that the relation of the four parameters was not clear if the slip length was used instead of flux. The reason why slip length may not directly relate to the flow is because the density and viscosity of the confined fluid are nonuniform.Using the data from the study of Liu and Li [75], we graphically illustrate the flow regimes in liquid in Figure 4 and summarize them theoretically in Table 2.Figure 4
(a) Variation of flux values as a function of fluid-wall interaction and fluid-fluid interaction, and (b) graphical representation of corresponding liquid flow regimes.
(a)
(b)Table 2
Different flow regimes as a function of liquid/solid binding force.
Interaction between
Flux
Flow regime
Physics
Liquid-liquid
Liquid-solid
Low
Low
High
Flux independent of solid-fluid binding force and increases with decreasing temperature
Surface effects are weak with no liquid adsorption. Flux is dominated by the temperature-driven friction between the solid and the liquid
Low
High
Intermediate
Flux independent of temperature and dependent on solid-fluid binding force
Surface effects become strong and liquid adsorption starts to take place that results in increasing fluid density close to the solid surface
High
Low
Intermediate
Flux independent of solid-fluid interaction and increases with decreasing temperature
Surface effects are weak with no liquid adsorption. Flux is dominated by the temperature-driven friction between the solid and the liquid
High
High
Low
Flux independent of solid-fluid interaction and increases with increasing temperature
Surface effects become strong and liquid adsorption starts to take place that results in increasing fluid density close to the solid surface. Flux is enhanced with increasing temperature as it results in desorption of the fluid molecules that merge together with the bulk liquid


## 3.1.1. Flow Regimes

The flow regimes in gases are based on the concept ofKn, which is the ratio of the mean free path of a molecule to the characteristic length of the medium. The spacing between molecules in liquids is much smaller than the spacing between molecules in gases, so the concept of mean free path and the flow regimes based on Kn are not valid in liquids. Because of the lack of a well-developed molecular theory for liquids, no flow regimes have been proposed for liquids in the literature. However, there is one study by Liu and Li [75] that talks about the liquid flow regimes in nanochannels. Liu and Li [75] investigated the liquid motion as a function of four parameters (fluid-fluid binding energy εff, fluid-solid wall binding energy εfw, temperature of the system T, and external driving force Fe corresponding to the pressure drop along the channel) by using nonequilibrium molecular dynamics simulations (NEMDs). They define two dimensionless numbers to understand the competing effect between the liquid-liquid interaction and liquid-surface interaction as εff/kBT and εfw/kBT×εfw/Feσfw, respectively, where kB is the Boltzmann constant and σfw is the collision diameter (separation) at which the potential vanishes. They compute the flux using liquid argon and liquid helium flowing through a channel made of silver walls at temperatures varying from 100 K to 300 K and pressure difference varying from 33 MPa to 132 MPa. They find liquid flux values as a function of these two dimensionless numbers and illustrate their results by dividing them into different regimes, each associated with a distinct mechanism. The authors noted that the relation of the four parameters was not clear if the slip length was used instead of flux. The reason why slip length may not directly relate to the flow is because the density and viscosity of the confined fluid are nonuniform.Using the data from the study of Liu and Li [75], we graphically illustrate the flow regimes in liquid in Figure 4 and summarize them theoretically in Table 2.Figure 4
(a) Variation of flux values as a function of fluid-wall interaction and fluid-fluid interaction, and (b) graphical representation of corresponding liquid flow regimes.
(a)
(b)Table 2
Different flow regimes as a function of liquid/solid binding force.
Interaction between
Flux
Flow regime
Physics
Liquid-liquid
Liquid-solid
Low
Low
High
Flux independent of solid-fluid binding force and increases with decreasing temperature
Surface effects are weak with no liquid adsorption. Flux is dominated by the temperature-driven friction between the solid and the liquid
Low
High
Intermediate
Flux independent of temperature and dependent on solid-fluid binding force
Surface effects become strong and liquid adsorption starts to take place that results in increasing fluid density close to the solid surface
High
Low
Intermediate
Flux independent of solid-fluid interaction and increases with decreasing temperature
Surface effects are weak with no liquid adsorption. Flux is dominated by the temperature-driven friction between the solid and the liquid
High
High
Low
Flux independent of solid-fluid interaction and increases with increasing temperature
Surface effects become strong and liquid adsorption starts to take place that results in increasing fluid density close to the solid surface. Flux is enhanced with increasing temperature as it results in desorption of the fluid molecules that merge together with the bulk liquid


## 3.2. Liquid Slip

Fundamental understanding of the mechanisms that control liquid slip is still elusive. Unlike gases, the concept of mean free path is not useful for liquids because the average distance between liquid molecules is comparable to the molecular diameter. Even though it is believed unanimously that the flow enhancement in nanoscale systems is caused by liquid slip velocity at the solid surface [76, 77], the debate on understanding its mechanisms mostly revolves around four key parameters. Many theories have been developed to associate liquid slip length (Ls) as a function of liquid properties [78, 79], interaction strength between the fluid and wall material [80–84], wettability of the medium [77, 85–89], surface roughness [90–95], and liquid shear rate [81, 86, 96–100].


### 3.2.1. Physical Mechanisms of Liquid Slip

The physical mechanisms behind the slip are still elusive, though there are few studies [101–107] that conjecture the reasons behind the liquid slip as due to the presence of either (i) a depleted water layer or (ii) an effective air gap at the wall formed by the nanobubbles. Another approach to hypothesize liquid slip is based on intermolecular forces affecting the momentum transport from the liquid to the wall [16]. Barring two studies [81, 108], there are no molecular-based theories to characterize this process [109, 110]. Lichter et al. [81] and Martini et al. [108] attribute the slip mechanism to hopping of liquid atoms from one equilibrium site to another, which follows Arrhenius dynamics based on the rate of hopping.


### 3.2.2. Slip Length

The slip velocity of the liquid at the surface boundary is characterized by an averaged property of all the liquid molecules called slip lengthLs, which is defined as the linearly extrapolated perpendicular distance from the surface boundary where fluid velocity becomes zero. In other words, slip length is the perpendicular distance from the solid boundary, outside the region of liquid flow, where the no-flow boundary condition is achieved. Figure 5 shows a schematic to illustrate the concept of slip for liquid flow; this figure shows three types of flow from left to right that indicate flow with no slip, flow with partial slip, and flow with full slip, respectively.Figure 5
Schematic to illustrate liquid slip on a stationary solid surface.Navier first introduced the linear boundary condition that is considered to be the standard and most basic characterization of liquid slip. Theoretically, liquid slip velocity (vs) is defined as follows:(25)vs=Lsn⋅dvxzdz+dvxzdzT⋅1−nn.Pressure gradient-driven fluid velocity (vx) in the presence of slip can then be theoretically calculated as follows [111]:(26)vx=h212η−dPdx1+6Lsh,where Ls is the slip length, n denotes the normal to the surface along the z direction, vs is the slip velocity at the surface, vx is the liquid velocity in a channel, η is liquid shear viscosity, dP/dx is the pressure gradient along x-direction, and h is the height of the channel.Research on liquid slip started with a focus on slip over smooth and hydrophobic surfaces [79, 105, 112, 113], but there has been an increasing evidence of liquid slip at different types of surface, including rough surface [90, 114] and hydrophilic surface [88, 115, 116]. Recent progress in liquid slip has broadened to various other parameters that have shown to influence liquid slip, and as a result, there is a continuous effort in developing appropriate models to characterize the slip on different types of surfaces and under different conditions. Below, we review the progress of liquid slip models as a function of (i) surface wettability, (ii) surface property, (iii) shear rate, and (iv) surface roughness.(1) Effect of Surface Wettability. The liquid slip at nanoscale is greatly affected by the surface wettability of the medium. Macroscopically, the wettability of a surface is classified as wetting, nonwetting, or neutral depending on the contact angle formed by a liquid droplet with the surface. A surface forming contact angles larger than 90° with a liquid is termed as nonwetting, whereas a surface forming contact angles smaller than 90° with a liquid is termed as wetting fluid. A surface forming contact angle equal to 90° with a liquid is termed as neutral. Nonwetting surfaces are known to have a high disposition towards introducing liquid slip [77, 86, 87, 117]. For example, hydrodynamic slip length of simple liquids flowing on smooth nonwetting surfaces (also referred as hydrophobic surfaces in the context of water as a liquid) is typically 20 nm in mediums with size varying from one to several hundreds of nanometers [112]. On the other hand, disposition of wetting surfaces in introducing slip is not as straight forward as it is for the nonwetting surfaces. There are reports [118] that suggest wetting surfaces are neutral in their disposition to introduce liquid slip, whereas others [88] have demonstrated that the wettability factor alone is not necessary to attain liquid slip and that it should be combined with the strength of the solid-fluid interactions. For example, if there is favorable sorption sites that are separated from each other by well-defined sub-nanometer distances, no slip takes place. However, liquid slip can occur if those favorable sorption sites are close to each other, provided liquid-solid attractions are not too strong. That is because nonwetting surfaces such as graphite are usually characterized by uniform distributions of fluid molecules, while the wetting surfaces such as crystalline silica are characterized by discontinuous sorption sites with high affinity for adsorbing fluid molecules. Voronov et al. [89] found that the slip length grows if the surface’s wettability is turned into more wetting by modifying the fluid-wall interfacial energy through affinity strength parameters. Gao and McCarthy [119] demonstrated experimentally that static contact angle of a fluid on a surface should not be the sole criteria to characterize the nature of wettability in order to determine the tendency of a liquid to slip. Instead, they suggested that contact angle hysteresis should be used to determine the nonwetting vs. wetting behavior of a surface.It is believed that Tolstoi [87] was the first person to propose the liquid slip by linking the macroscopic thermodynamics concept at the molecular scale. He quantified the mobility of liquid molecules with the surface energies (=σ2γ) through spreading coefficient (=γ1−cosθ) as follows:(27)Lsσ∼expασ2γ1−cosθkBT−1,where α is a dimensionless geometrical parameter of first order that represents the microcavity area within the solid, σ is molecular diameter, γ is the surface tension between the liquid and the surface, θ is the equilibrium/static contact angle, T is the absolute temperature, kB is the Boltzmann constant. Another theory [85, 120, 121] that describes liquid slip at the molecular level with thermodynamic equilibrium uses the Green-Kubo relation and fluctuation-dissipation theorem (linear regression of fluctuations) and is given by the following equation:(28)Lsσ∼D∗StcLS2ρcσ3,where D∗=Dp/D0, Dp is the collective molecular diffusion coefficient, D0 is the bulk diffusion coefficient, St is a factor that represents first molecular layer structure, ρc is the fluid density at the first molecular layer, σ is the molecular diameter, and cLS is the dimensionless coefficient for solid-liquid Lennard-Jones potential.Huang et al. [122] and Sendner et al. [116] studied the slip length of shear-driven water on various hydrophobic surfaces using nonequilibrium molecular dynamic simulations (NEMDs). They show that the slip length as a function of different surfaces (both organic-like silane monolayers and inorganic-like Lennard-Jones models) can be characterized using the static contact angle (θ) of the surface as shown below, which they describe as the crucial parameter in controlling the water slippage:(29)Ls∝1+cosθ−2.The authors demonstrate that their formula can also be obtained from the linear response theory by using the LJ fluid-solid energy parameter and the Young-Laplace equation. Their estimated slippage length varies from a few nanometers to tens of nanometers. They also proposed the following scaling expression for slip length as a function of liquid depletion length, which they define as the layer between the water and the solid surface occupied by the gas:(30)Ls∝δ4,where δ is the thickness of liquid depletion and given by(31)δ=∫0∞1−ρszρsb−ρlzρlbdz,where ρsz and ρlz are the densities of the solid and liquid along the height of the channel, whereas ρsb and ρlb are the bulk densities of the solid and liquid, respectively.Bakli and Chakraborty [123] presented a constitutive model to describe the nanopore imbibition characteristics of water as a function of the dynamic slip length using NEMDs. Their derived slip length depends on the contact angle and an acceleration parameter that comes from the pressure gradients due to capillary forces. They show that for low acceleration values, their model approaches the model given by Huang et al. [122]:(32)Ls=A1a1+cosθ2−A2a,a=γcosθρRL,a∗=γcosθηLdL/dt,arough∗=asmooth∗+a¯,a¯=0.35exp−0.6N,where L is the capillary imbibition length at time t, N is the parameter representing the surface roughness (lower value represents higher roughness), a is the imbibition driving acceleration parameter, a∗ is the dimensionless form of a, a¯ is a deceleration parameter to mimic flow hindrance due to surface roughness, and A1a and A2a are functions of a.(2) Effect of Shear Rate. The effect of shear rate on slip length has been investigated experimentally by many researchers [90, 96, 97, 99, 100] and observed to depend on the shear rate. In such a scenario where slip length is found to depend on the shear rate, the slip boundary condition proposed by Navier, and presented earlier, shows a nonlinear behavior with shear rate. Although several slip length relationships have been proposed as a function of shear rate, the relationship proposed by Priezjev [98] is an important one as it accounts for the solid-fluid interaction strength that is crucial in micro- or nanoporous media. Priezjev [98] proposed slip length for smooth surfaces as a function of shear rate that predicts that the slip length increases nonlinearly with the shear rate for weak solid-fluid interactions provided that the solid-liquid interface forms a discrete structure, whereas if the solid-fluid interaction increases, the slip length becomes linearly proportional to the shear rate:(33)Qslip=23ρfxh3μ+2hVs.The above equation can be normalized in terms of the no-slip flow rate as follows:(34)Qslip=Qno−slip1+6Ls2h,where fx is an externally applied force in the x direction, ρ is the fluid density, h is the distance between the confining solid walls, μ is the fluid viscosity, Vs is the fluid velocity at the confining parallel walls that can also be written as Vs=v−h=vh, and Qslip and Qno−slip are the flow rates due to slip and no-slip boundary conditions, respectively. The term 2hVs in the above equation represents the correction to the flow rate due to slip boundary conditions. Priezjev [98] reported that the fluid viscosity, shown below, was found to be shear rate (γ˙) independent in the range of γ˙≤0.16τ−1, where τ is the characteristic Lennard-Jones (LJ) time:(35)μ=2.2±0.2ετσ−3,where ε and σ represent the energy and length scales of the fluid phase.(3) Effect of Surface Roughness. Even though most of the solids are naturally rough at micro- to nanoscale, slip length of liquid as a function of surface roughness is not well studied [16, 124]. The problem of calculating slip as a function of surface roughness is challenging [16, 125, 126] because of four reasons: (i) controlling surface roughness in a natural and a realistic manner at nanometer scale is not easy, (ii) surface roughness can introduce changes to some interfacial properties such as interfacial forces, wettability, and trapped gases, which may be undesirable, (iii) rough surface will introduce undulating boundaries that will introduce complexities in making correct interpretations with respect to varying boundary positions, and (iv) there is no analytical characterization for the nonperiodic and natural randomness of surface roughness.Although several slip length relationships have been proposed as a function of surface roughness, two relationships by Lund et al. [127] and Misra and Bakli [128] are discussed here because of their simplicity and effectiveness.Lund et al. [127] derived an expression of effective slip length on a surface with periodic roughness in terms of intrinsic slip length and contact area of the surface, using steady-state Stokes flow and homogenization method. Their derived effective slip length represents the harmonic mean weighted by the area of contact between surface and fluid as follows:(36)Leff=1+s2Ls−1,where Ls is the intrinsic slip length, s is the slope, and 1+s2=1/L∫0L1+s2dx is the arc length of the surface function over one period, normalized by the length of the period (L). For a flat surface (s=0), the above expression reduces to the intrinsic slip length. Their result is applicable for rough surfaces where L≪Lsmin.Misra and Bakli [128] recently proposed an expression of effective slip length for slippage of water considering coupled effect of the surface chemistry and surface roughness by performing molecular dynamics simulations. They found that the slip length without surface roughness match closely with the result of 1+cosθ2 proposed by Huang et al. [122] for smooth surface, while the presence of surface roughness always reduces the effective slip length as follows:(37)Leff=k1+cosθ21+α∗2,α∗=∑fx,y−fx,y¯N,where k is the constant of proportionality and α∗ is a parameter that characterizes the roughness of a surface through its amplitude N and roughness function fx,y that modifies the attractive part of the Lennard-Jones (LJ) potential.(4) Electrokinetic Effects. The amount of liquid slip on a surface is found to be affected by the presence of a charged surface and/or liquids with charged ions such as electrolytes. The presence of charge leads to attraction of oppositely charged ions and repulsion of similarly charged ions, and in the presence of a charged surface, this leads to the formation of a charged interfacial layer near the surface. The fluid flow in the presence of charge (also referred to as electrokinetic phenomenon) is centrally connected to the key concept of electrical double layer (EDL) that is made up of a stern layer and a diffuse layer as illustrated in Figure 6. Stern layer is the width of a layer close to the surface where the underlying solid surface charge (cations or anions) attracts the opposite charge of the solvent/liquid to form a layer of highly concentrated ions, whereas the diffuse layer is the width of a layer starting from the stern layer and ending where the electrokinetic potential goes to zero.Figure 6
Schematic of a flow in a channel with electrokinetic effects.The electrokinetic effects result in enhanced solid-fluid interaction close to the surface, which may also influence the liquid slip flow and consequently the overall fluid transport. Two examples of electrokinetic flows relevant here are electro-osmosis (liquid flow induced by application of an electric field) and streaming currents (electric current induced by application of a pressure gradient).The surface-fluid interaction due to electrokinetic effects is more complex than the surface-fluid interaction without any charge. Research on this area is sparse due to limited applications that have surface charge present. Out of the few models that exist in the literature, models proposed by Joly et al. [129] and Choudhary et al. [130] are discussed here because of parameters that are directly related to the electrokinetic effects and solid-fluid interaction compared to other models where the parameters are not that apparent.Joly et al. [129] explored the effect of surface charge to solid-liquid interfacial friction and its effect on the slip length boundary condition. They study the hydrodynamics of liquid coupled with electrostatics within the EDL using molecular dynamics simulations. Their study found that the presence of surface charge on the solid-liquid interface enhances the frictional force between the solid-liquid and as a result decreases the slip length.(38)Leff=Ls1+1/ασd2/e2lB/d2Ls,where Ls is the slip length without considering surface charge, α≈1 is a numerical factor, σ is the surface charge density, d is the equilibrium distance of Lennard-Jones potential, e is the elementary charge, and lB (∼0.7 nm in water at room temperature) is the Bjerrum length that represents a length scale below which electrostatic interactions dominate thermal effects.Choudhary et al. [130] studied the liquid slip for electro-osmotic flow through a nanochannel with hydrophobic walls of sinusoidally varying slippage using an asymptotic theory that uses the ratio of pattern amplitude (perpendicular to the applied electric field) to the average slip as a parameter. They proposed the following analytical expression for effective slip that can be used to design slip through engineered variations in topography and/or chemistry:(39)L¯eff=L¯s−α2L¯s2k¯s2sinhk¯s1+cosθ2coshk¯s+L¯s.k¯ssinhk¯s+coshk¯s1−cosθ2sinhk¯s+L¯s.k¯scoshk¯s,where L¯eff is the dimensionless effective slip length due to electrokinetic effects and L¯s is the dimensionless average slip length (=Lsx/h) that is equal to Ls01+αcosk¯sx¯+θ and Ls01+αcosk¯sx¯ on the top and bottom walls of the channel with height 2h, respectively. Here, Ls0 is the average slip length without electrokinetic effects, α is the dimensionless amplitude (=a/Ls0) of the pattern on either walls charged to facilitate electro-osmosis flow, k¯s is the wavenumber, and θ is the phase angle between the slip waves on the two walls.


## 3.2.1. Physical Mechanisms of Liquid Slip

The physical mechanisms behind the slip are still elusive, though there are few studies [101–107] that conjecture the reasons behind the liquid slip as due to the presence of either (i) a depleted water layer or (ii) an effective air gap at the wall formed by the nanobubbles. Another approach to hypothesize liquid slip is based on intermolecular forces affecting the momentum transport from the liquid to the wall [16]. Barring two studies [81, 108], there are no molecular-based theories to characterize this process [109, 110]. Lichter et al. [81] and Martini et al. [108] attribute the slip mechanism to hopping of liquid atoms from one equilibrium site to another, which follows Arrhenius dynamics based on the rate of hopping.


## 3.2.2. Slip Length

The slip velocity of the liquid at the surface boundary is characterized by an averaged property of all the liquid molecules called slip lengthLs, which is defined as the linearly extrapolated perpendicular distance from the surface boundary where fluid velocity becomes zero. In other words, slip length is the perpendicular distance from the solid boundary, outside the region of liquid flow, where the no-flow boundary condition is achieved. Figure 5 shows a schematic to illustrate the concept of slip for liquid flow; this figure shows three types of flow from left to right that indicate flow with no slip, flow with partial slip, and flow with full slip, respectively.Figure 5
Schematic to illustrate liquid slip on a stationary solid surface.Navier first introduced the linear boundary condition that is considered to be the standard and most basic characterization of liquid slip. Theoretically, liquid slip velocity (vs) is defined as follows:(25)vs=Lsn⋅dvxzdz+dvxzdzT⋅1−nn.Pressure gradient-driven fluid velocity (vx) in the presence of slip can then be theoretically calculated as follows [111]:(26)vx=h212η−dPdx1+6Lsh,where Ls is the slip length, n denotes the normal to the surface along the z direction, vs is the slip velocity at the surface, vx is the liquid velocity in a channel, η is liquid shear viscosity, dP/dx is the pressure gradient along x-direction, and h is the height of the channel.Research on liquid slip started with a focus on slip over smooth and hydrophobic surfaces [79, 105, 112, 113], but there has been an increasing evidence of liquid slip at different types of surface, including rough surface [90, 114] and hydrophilic surface [88, 115, 116]. Recent progress in liquid slip has broadened to various other parameters that have shown to influence liquid slip, and as a result, there is a continuous effort in developing appropriate models to characterize the slip on different types of surfaces and under different conditions. Below, we review the progress of liquid slip models as a function of (i) surface wettability, (ii) surface property, (iii) shear rate, and (iv) surface roughness.(1) Effect of Surface Wettability. The liquid slip at nanoscale is greatly affected by the surface wettability of the medium. Macroscopically, the wettability of a surface is classified as wetting, nonwetting, or neutral depending on the contact angle formed by a liquid droplet with the surface. A surface forming contact angles larger than 90° with a liquid is termed as nonwetting, whereas a surface forming contact angles smaller than 90° with a liquid is termed as wetting fluid. A surface forming contact angle equal to 90° with a liquid is termed as neutral. Nonwetting surfaces are known to have a high disposition towards introducing liquid slip [77, 86, 87, 117]. For example, hydrodynamic slip length of simple liquids flowing on smooth nonwetting surfaces (also referred as hydrophobic surfaces in the context of water as a liquid) is typically 20 nm in mediums with size varying from one to several hundreds of nanometers [112]. On the other hand, disposition of wetting surfaces in introducing slip is not as straight forward as it is for the nonwetting surfaces. There are reports [118] that suggest wetting surfaces are neutral in their disposition to introduce liquid slip, whereas others [88] have demonstrated that the wettability factor alone is not necessary to attain liquid slip and that it should be combined with the strength of the solid-fluid interactions. For example, if there is favorable sorption sites that are separated from each other by well-defined sub-nanometer distances, no slip takes place. However, liquid slip can occur if those favorable sorption sites are close to each other, provided liquid-solid attractions are not too strong. That is because nonwetting surfaces such as graphite are usually characterized by uniform distributions of fluid molecules, while the wetting surfaces such as crystalline silica are characterized by discontinuous sorption sites with high affinity for adsorbing fluid molecules. Voronov et al. [89] found that the slip length grows if the surface’s wettability is turned into more wetting by modifying the fluid-wall interfacial energy through affinity strength parameters. Gao and McCarthy [119] demonstrated experimentally that static contact angle of a fluid on a surface should not be the sole criteria to characterize the nature of wettability in order to determine the tendency of a liquid to slip. Instead, they suggested that contact angle hysteresis should be used to determine the nonwetting vs. wetting behavior of a surface.It is believed that Tolstoi [87] was the first person to propose the liquid slip by linking the macroscopic thermodynamics concept at the molecular scale. He quantified the mobility of liquid molecules with the surface energies (=σ2γ) through spreading coefficient (=γ1−cosθ) as follows:(27)Lsσ∼expασ2γ1−cosθkBT−1,where α is a dimensionless geometrical parameter of first order that represents the microcavity area within the solid, σ is molecular diameter, γ is the surface tension between the liquid and the surface, θ is the equilibrium/static contact angle, T is the absolute temperature, kB is the Boltzmann constant. Another theory [85, 120, 121] that describes liquid slip at the molecular level with thermodynamic equilibrium uses the Green-Kubo relation and fluctuation-dissipation theorem (linear regression of fluctuations) and is given by the following equation:(28)Lsσ∼D∗StcLS2ρcσ3,where D∗=Dp/D0, Dp is the collective molecular diffusion coefficient, D0 is the bulk diffusion coefficient, St is a factor that represents first molecular layer structure, ρc is the fluid density at the first molecular layer, σ is the molecular diameter, and cLS is the dimensionless coefficient for solid-liquid Lennard-Jones potential.Huang et al. [122] and Sendner et al. [116] studied the slip length of shear-driven water on various hydrophobic surfaces using nonequilibrium molecular dynamic simulations (NEMDs). They show that the slip length as a function of different surfaces (both organic-like silane monolayers and inorganic-like Lennard-Jones models) can be characterized using the static contact angle (θ) of the surface as shown below, which they describe as the crucial parameter in controlling the water slippage:(29)Ls∝1+cosθ−2.The authors demonstrate that their formula can also be obtained from the linear response theory by using the LJ fluid-solid energy parameter and the Young-Laplace equation. Their estimated slippage length varies from a few nanometers to tens of nanometers. They also proposed the following scaling expression for slip length as a function of liquid depletion length, which they define as the layer between the water and the solid surface occupied by the gas:(30)Ls∝δ4,where δ is the thickness of liquid depletion and given by(31)δ=∫0∞1−ρszρsb−ρlzρlbdz,where ρsz and ρlz are the densities of the solid and liquid along the height of the channel, whereas ρsb and ρlb are the bulk densities of the solid and liquid, respectively.Bakli and Chakraborty [123] presented a constitutive model to describe the nanopore imbibition characteristics of water as a function of the dynamic slip length using NEMDs. Their derived slip length depends on the contact angle and an acceleration parameter that comes from the pressure gradients due to capillary forces. They show that for low acceleration values, their model approaches the model given by Huang et al. [122]:(32)Ls=A1a1+cosθ2−A2a,a=γcosθρRL,a∗=γcosθηLdL/dt,arough∗=asmooth∗+a¯,a¯=0.35exp−0.6N,where L is the capillary imbibition length at time t, N is the parameter representing the surface roughness (lower value represents higher roughness), a is the imbibition driving acceleration parameter, a∗ is the dimensionless form of a, a¯ is a deceleration parameter to mimic flow hindrance due to surface roughness, and A1a and A2a are functions of a.(2) Effect of Shear Rate. The effect of shear rate on slip length has been investigated experimentally by many researchers [90, 96, 97, 99, 100] and observed to depend on the shear rate. In such a scenario where slip length is found to depend on the shear rate, the slip boundary condition proposed by Navier, and presented earlier, shows a nonlinear behavior with shear rate. Although several slip length relationships have been proposed as a function of shear rate, the relationship proposed by Priezjev [98] is an important one as it accounts for the solid-fluid interaction strength that is crucial in micro- or nanoporous media. Priezjev [98] proposed slip length for smooth surfaces as a function of shear rate that predicts that the slip length increases nonlinearly with the shear rate for weak solid-fluid interactions provided that the solid-liquid interface forms a discrete structure, whereas if the solid-fluid interaction increases, the slip length becomes linearly proportional to the shear rate:(33)Qslip=23ρfxh3μ+2hVs.The above equation can be normalized in terms of the no-slip flow rate as follows:(34)Qslip=Qno−slip1+6Ls2h,where fx is an externally applied force in the x direction, ρ is the fluid density, h is the distance between the confining solid walls, μ is the fluid viscosity, Vs is the fluid velocity at the confining parallel walls that can also be written as Vs=v−h=vh, and Qslip and Qno−slip are the flow rates due to slip and no-slip boundary conditions, respectively. The term 2hVs in the above equation represents the correction to the flow rate due to slip boundary conditions. Priezjev [98] reported that the fluid viscosity, shown below, was found to be shear rate (γ˙) independent in the range of γ˙≤0.16τ−1, where τ is the characteristic Lennard-Jones (LJ) time:(35)μ=2.2±0.2ετσ−3,where ε and σ represent the energy and length scales of the fluid phase.(3) Effect of Surface Roughness. Even though most of the solids are naturally rough at micro- to nanoscale, slip length of liquid as a function of surface roughness is not well studied [16, 124]. The problem of calculating slip as a function of surface roughness is challenging [16, 125, 126] because of four reasons: (i) controlling surface roughness in a natural and a realistic manner at nanometer scale is not easy, (ii) surface roughness can introduce changes to some interfacial properties such as interfacial forces, wettability, and trapped gases, which may be undesirable, (iii) rough surface will introduce undulating boundaries that will introduce complexities in making correct interpretations with respect to varying boundary positions, and (iv) there is no analytical characterization for the nonperiodic and natural randomness of surface roughness.Although several slip length relationships have been proposed as a function of surface roughness, two relationships by Lund et al. [127] and Misra and Bakli [128] are discussed here because of their simplicity and effectiveness.Lund et al. [127] derived an expression of effective slip length on a surface with periodic roughness in terms of intrinsic slip length and contact area of the surface, using steady-state Stokes flow and homogenization method. Their derived effective slip length represents the harmonic mean weighted by the area of contact between surface and fluid as follows:(36)Leff=1+s2Ls−1,where Ls is the intrinsic slip length, s is the slope, and 1+s2=1/L∫0L1+s2dx is the arc length of the surface function over one period, normalized by the length of the period (L). For a flat surface (s=0), the above expression reduces to the intrinsic slip length. Their result is applicable for rough surfaces where L≪Lsmin.Misra and Bakli [128] recently proposed an expression of effective slip length for slippage of water considering coupled effect of the surface chemistry and surface roughness by performing molecular dynamics simulations. They found that the slip length without surface roughness match closely with the result of 1+cosθ2 proposed by Huang et al. [122] for smooth surface, while the presence of surface roughness always reduces the effective slip length as follows:(37)Leff=k1+cosθ21+α∗2,α∗=∑fx,y−fx,y¯N,where k is the constant of proportionality and α∗ is a parameter that characterizes the roughness of a surface through its amplitude N and roughness function fx,y that modifies the attractive part of the Lennard-Jones (LJ) potential.(4) Electrokinetic Effects. The amount of liquid slip on a surface is found to be affected by the presence of a charged surface and/or liquids with charged ions such as electrolytes. The presence of charge leads to attraction of oppositely charged ions and repulsion of similarly charged ions, and in the presence of a charged surface, this leads to the formation of a charged interfacial layer near the surface. The fluid flow in the presence of charge (also referred to as electrokinetic phenomenon) is centrally connected to the key concept of electrical double layer (EDL) that is made up of a stern layer and a diffuse layer as illustrated in Figure 6. Stern layer is the width of a layer close to the surface where the underlying solid surface charge (cations or anions) attracts the opposite charge of the solvent/liquid to form a layer of highly concentrated ions, whereas the diffuse layer is the width of a layer starting from the stern layer and ending where the electrokinetic potential goes to zero.Figure 6
Schematic of a flow in a channel with electrokinetic effects.The electrokinetic effects result in enhanced solid-fluid interaction close to the surface, which may also influence the liquid slip flow and consequently the overall fluid transport. Two examples of electrokinetic flows relevant here are electro-osmosis (liquid flow induced by application of an electric field) and streaming currents (electric current induced by application of a pressure gradient).The surface-fluid interaction due to electrokinetic effects is more complex than the surface-fluid interaction without any charge. Research on this area is sparse due to limited applications that have surface charge present. Out of the few models that exist in the literature, models proposed by Joly et al. [129] and Choudhary et al. [130] are discussed here because of parameters that are directly related to the electrokinetic effects and solid-fluid interaction compared to other models where the parameters are not that apparent.Joly et al. [129] explored the effect of surface charge to solid-liquid interfacial friction and its effect on the slip length boundary condition. They study the hydrodynamics of liquid coupled with electrostatics within the EDL using molecular dynamics simulations. Their study found that the presence of surface charge on the solid-liquid interface enhances the frictional force between the solid-liquid and as a result decreases the slip length.(38)Leff=Ls1+1/ασd2/e2lB/d2Ls,where Ls is the slip length without considering surface charge, α≈1 is a numerical factor, σ is the surface charge density, d is the equilibrium distance of Lennard-Jones potential, e is the elementary charge, and lB (∼0.7 nm in water at room temperature) is the Bjerrum length that represents a length scale below which electrostatic interactions dominate thermal effects.Choudhary et al. [130] studied the liquid slip for electro-osmotic flow through a nanochannel with hydrophobic walls of sinusoidally varying slippage using an asymptotic theory that uses the ratio of pattern amplitude (perpendicular to the applied electric field) to the average slip as a parameter. They proposed the following analytical expression for effective slip that can be used to design slip through engineered variations in topography and/or chemistry:(39)L¯eff=L¯s−α2L¯s2k¯s2sinhk¯s1+cosθ2coshk¯s+L¯s.k¯ssinhk¯s+coshk¯s1−cosθ2sinhk¯s+L¯s.k¯scoshk¯s,where L¯eff is the dimensionless effective slip length due to electrokinetic effects and L¯s is the dimensionless average slip length (=Lsx/h) that is equal to Ls01+αcosk¯sx¯+θ and Ls01+αcosk¯sx¯ on the top and bottom walls of the channel with height 2h, respectively. Here, Ls0 is the average slip length without electrokinetic effects, α is the dimensionless amplitude (=a/Ls0) of the pattern on either walls charged to facilitate electro-osmosis flow, k¯s is the wavenumber, and θ is the phase angle between the slip waves on the two walls.


## 4. Porous Media Applications in the Energy Sector

In this section, the above-reviewed knowledge on gas and liquid flow at micro- to nanoscale is used to critically examine the modeling gaps in four different porous media applications from the energy sector and how that can be improved. Specifically, for each application we provide a brief introduction in context of its importance, a brief literature review on that application with the present state of progress, followed by identification of critical gaps and proposed approach to improve modeling at micro- and nanoscale in each application. The four applications from energy sector considered here are (i) electrokinetic energy conversion devices, (ii) membrane-based water desalination through reverse osmosis, (iii) shale reservoirs, and (iv) hydrogen storage, respectively.


### 4.1. Energy Conversion through Electrokinetics

#### 4.1.1. Introduction

Increasing demand for energy has led to harnessing of energy from various forms such as mechanical, chemical, internal, electrical, thermal, etc. and various other combinations. To utilize these various sources of energy, micro- and nanofluidics provide the capability of converting energy from one form to another form to cater to different applications. Energy conversion is possible through electrokinetics, which is a term used to describe different transport mechanisms in the presence of an electric field and liquid flow through a medium. Electrokinetics can be used for energy conversion in four different ways: (i) drive liquid by applying an electric potential (electro-osmosis), (ii) transport and separate particles by applying an electric potential (electrophoresis), (iii) generate electric potential by pressure gradient-driven liquid flow (streaming potential), and (iv) generate electric potential by movement of charged particles under gravitational or centrifugal force (sedimentation potential). In particular, micro- and nanofluidics have found application in converting energy from heat, electrical, chemical/biochemical, pressure, etc. to electrical energy using micro- and nanoscale technologies with a wide range of efficiencies varying from 3% for streaming potential devices [131] to 60% for micro-fuel cells [132].


#### 4.1.2. Literature Review

Electric double layer (∼1 nm thickness) is a common phenomenon in small energy conversion devices that are increasingly based on nanostructured materials. As an example to illustrate the scale at which fluids flow in electrokinetic devices, Figure7 shows transmission electron microscope images of carbon composites [133, 134] with application in electrochemical energy conversion.Figure 7
Transmission electron microscope images of carbon composites for application in electrochemical energy conversion [133, 134].The electrokinetic effects at confined scales result in enhanced solid-fluid interaction close to the surface, which may also influence the liquid slip flow and consequently the overall fluid transport. The amount of liquid slip on a surface is found to be affected by the presence of a charged surface and/or liquids with charged ions such as electrolytes. Two examples of electrokinetic flows relevant here are electro-osmosis (liquid flow induced by application of an electric field) and streaming currents (electric current induced by application of a pressure gradient).(1) Key Challenge. One of the key challenges in the design of electrokinetic energy conversion devices is their low conversion efficiencies. The maximum efficiency of energy conversion ηmax, which is defined as the ratio of the output power to the input power, is given as follows:(40)ηmax=PoweroutputPowerinput.(2) Governing Equations. The mathematical theory used to model the transport in the presence of chemical potential (due to surface charge and charged ions in the solution), externally applied electric potential, and hydrostatic potential is briefly reviewed below. The chemical potential and ionic concentrations in the diffuse part of the EDL are governed by the Poisson-Boltzmann equation as follows:(41)∇2ψy=−1ερy=−eε∑i=1Nziciexp−eziψykBT,where ψy is the potential in the EDL due to surface charge and the charged ions in the solution, y is the radial distance from the surface of the channel, ρy is the volume charge density of all the ions present in the neighborhood of the surface, e is the unit charge, ε is the permittivity of the liquid, ci and zi are the concentration and the ionic valence of the ith ionic species, respectively, kB is the Boltzmann constant, and T is the absolute temperature.The Poisson-Boltzmann equation does not have a general analytical solution, but its solutions are available for specific cases [135]. For our interest, we consider the solution for a flat surface with a low potential (ψsurface) at its surface for which the above equation can be linearized based on the Debye-Hückel approximation to obtain the following solution:(42)ψy=ψsurfaceexp−κy,λD=κ−1=εkBTNAe2∑i=1Nzi2ci,where κ−1 is called the Debye length λD that gives a measure of the EDL and NA is the Avogadro number. The Debye length is independent of the charge on the solid surface and increases as the solute ion concentration in the solution becomes dilute such that λD for a solution without electrolytes can be considered as infinitely thick. Typical values of λD can vary from ∼3nm to ∼300nm for electrolytes with ionic concentration of 10−2M and 10−6M, respectively [135].The above-defined definitions are finally used to describe the solvent and solute transport in the presence of hydrostatic, chemical, and electric potential within a micro- or a nanochannel of width2h as follows.(i) Solvent Transport. The solvent flux across the channel due to hydrostatic potential and electric potential is the sum of the average of the fluxes due to these two potentials and given as follows [136, 137]:(43)vs=1h∫0hvsydy,vs=1h∫0h12ηy2−2hydPdx︸flow  due  tohydrostatic  potentialdy+1h∫0hεηζ−ψydϕdx︸flow  due  toelectric  potentialdy,where vs is the average solvent velocity, ζ is the electric potential at the slipping plane (also called zeta potential), ψy is the potential in the EDL, η is the solution viscosity, ε is the permittivity of the solution, and ϕ is an externally applied electric potential.(ii) Solute Transport. The solute flux across the channel is given as the averaged sum of the flux due to hydrostatic potential, flux due to diffusion, and flux due to electric potential. The solute transport in the presence of these three mechanisms is also called the Nernst-Planck equation and given as follows [137–139]:(44)J=1h∫0hcyvsy︸flux  due  tohydrostatic  potentialdy−1h∫0hD∞dcydx︸flux  due  todiffusiondy−1h∫0hzeD∞kBTcydϕdx︸flux  due  toelectric  potentialdy,where J is the area-averaged solute flux, D∞ is the diffusion coefficient of the solute in the free solution, cy is the solute concentration at distance y from the surface, and ze is the ionic charge.(iii) Solute Partitioning. Solute ions experience electric, chemical, and surface forces that vary in magnitude based on the size and charge on the ions; therefore, this leads to nonuniform distribution of ions along the width of the flowing medium. This partitioning of the solute ions can be characterized through the ratio of the average solute concentration c to the bulk concentration c0 in the solution, which is also referred to as the partitioning coefficient βpart [140]:(45)βpart=cc0=1h∫∼0h−rexp−EykBTdy,where βpart is the solute partition coefficient that is the measure of ion permeability-selectivity of the pore, r is the radius of the spherical solute, and E is the electrostatic interaction potential. It should be noted that when the confinement width becomes comparable to the Debye length, the overlapped EDL gets deformed significantly due to steric effects and in that case the modified Poisson-Boltzmann equation must be employed [141]. When the transport of solute through the medium is strongly dependent on solute size and independent of frictional and hydrostatic (convective) forces, in that case [137, 142],(46)βpart=1−rh.(iv) Conversion Efficiency. The maximum efficiency of energy conversion ηmax is defined as the ratio of the output power to the input power; however, some researchers [143–145] also use an alternate measure for the conversion efficiency using a measure called electrokinetic figure of merit [146], which is based on phenomenological transport equations. The following equation is an alternate form for the maximum efficiency in terms of figure of merit β:(47)ηmax=1+β−11+β+1,where β=ν2σ/kH. Here, ν, σ, and kH are the streaming potential coefficient, the ion conductivity, and the hydraulic permeability, respectively.


#### 4.1.3. Critical Gaps and Proposed Approach

About a decade ago it was hypothesized [147, 148] that the efficiencies of the energy conversion devices based on electrokinetic effects can be optimized by enhancing the liquid slip. A theoretical study by Davidson and Xuan [143] showed that the efficiency of electrokinetic energy conversion devices is almost a linear function of slip length for small values of zeta potentials. Even though there is a consensus that liquid slip plays a strong role in the device efficiency among the electrokinetic energy research community, some issues remain in regard to understanding of slip length as a function of key device parameters and its relationship with flux.(1) Slip Length as a Function of Device Parameters. In the electrokinetic energy research community, performance of the device is predicted by either completely neglecting the slip [149, 150] or an assumed value of slip length used by others [143, 148, 151, 152] that is independent of the surface and flow conditions. Further, some researchers [153, 154] suggested that the slip length increases as the surface becomes slippery, or in other words slip length increases on surfaces with low friction coefficient. However, this suggestion is contrary to the observations reported by researchers [155–158] in the field of liquid slip who reported that the surface friction results in increase in the effective viscosity that consequently increases the liquid slip length [156]. This result, even though counterintuitive in nature, is supported by experimental observations [90, 96, 99, 159]. Therefore, as noted by others [160], there is a lack of understanding about the liquid slip length when it comes to predicting the conversion efficiencies. Developing the fundamental understanding of liquid slip as a function of key device parameters presents an opportunity of increasing the efficiency of energy conversion devices by exploiting their potential of developing large slip lengths. For instance, we discussed many slip length models earlier as a function of surface wettability, surface friction, shear rate, surface roughness, and electrokinetic effects that can be used to either estimate the slip lengths for a given device or design the device parameters for a given slip length.(2) Slip Length vs. Flux. Following theoretical considerations initiated by some studies [143, 147, 148], it is now believed by most researchers that the efficiency of energy conversion devices can be increased by increasing the liquid slip length. In this regard, some comments were made on how to optimally exploit this potential. Having said that, the arguments made in favor of slip are theoretical in nature, and there are no experimental studies specific to energy conversion devices that verify the amplified values of conversion efficiencies predicted theoretically and usually cited by many authors [137, 153]. The objective of this note is to caution that a theoretical definition of liquid slip where flux is proportional to the slip length (vx∝1+6Ls/h) may not be universally true. To support this hypothesis, we refer to the experimental study of Liu and Li [75] who measured the flux and slip length for liquid argon and liquid helium in nanochannels of various widths and noted that the relation of slip length to flux was not clear from the data. They attributed this reason to nonuniform distribution of density and viscosity of the confined fluid.


### 4.2. Membrane-Based Water Desalination

#### 4.2.1. Introduction

Currently, there are 4 major sources of water that can be used for water desalination and treatment, i.e., produced water from oil and gas wells, brackish groundwater, municipal wastewater, and seawater, respectively. Out of these, seawater is a source with continuous and large amounts of saline water that can be potentially desalinated and supplied to various applications on large scale as illustrated in Figure8.Figure 8
Role of water in the energy sector.Desalination using the membrane-based reverse osmosis process is a potentially attractive technology that can be used to desalinate enormous amount of seawater. Reverse osmosis (RO) is considered to be a leading technology for desalination primarily because of its relatively low cost (one-half to one-third of the cost of distillation) due to relatively low energy consumption (1.8 kW·h/m3) that contributes to about 44% of production capacity in the world [161]. RO uses a semipermeable membrane to stop the contaminants at the level of 0.1 nm to 1 nm by applying hydrostatic pressure to overcome osmotic pressure, which is a chemical potential developed due to salinity difference of the solvents that drives the solvent from an aqueous solution of low salt concentration to an aqueous solution of high salt concentration. In the RO, the movement of solvent from low concentration to high concentration due to chemical potential is reversed by applying an external pressure.


#### 4.2.2. Literature Review

The study by Nair et al. [162] was the first to discover that graphene oxide has the membrane properties suitable for desalination, i.e., unique water permeation pathway and retention for other gases and ions/solutes. The water flows through the multiple layers of graphene oxide through two-dimensional capillaries formed by closely spaced graphene sheets where the salt ions are retained on the sheets as shown in Figure 9.Figure 9
Graphene oxide membrane with nanopores [162].Conventional understanding of the flow fluid through pores suggests that this dual-task is paradoxical in nature, but development of nanomaterials such as carbon nanotubes (CNTs), boron nitride nanotubes (BNNTs), and graphene that show unconventional fluid transport behavior due to their smallest possible pore volumes combined with ionic dehydration [163], electrokinetic effects [164], and size exclusion [165] properties provide opportunity for their application in optimizing two coupled tasks of desalination mentioned above. Specifically, high selective-solute retention and high permeability to water is achieved by drilling sub-nanometer holes through the solid-state nanomaterials such as graphene by using focus ion beam [166, 167] or oxygen plasma etching process [168]. The water molecules are too big to pass through graphene’s fine mesh without making holes/pores in the graphene sheet and poking these sub-nanometer holes uniformly is another challenge. Additionally, in order to drive satisfactory magnitude of water flux through these nanopores, a large pressure difference is required that cannot be solely achieved by osmotic pressure. However, compared to the fine mesh of nanoporous graphene, Nair et al. [162] proposed graphene oxide membrane as an alternate to the graphene because graphene oxide relies on its naturally occurring tortuous path to flow out water (Figure 9), unlike the challenging task of artificially poking holes through graphene for creating a path to flow out water. Other than this major advantage of graphene oxide, there are several other reasons as presented by You et al. [169] that make graphene oxide a preferred choice over nanoporous graphene.(1) Key Challenge. One of the key design challenges in the desalination technology using RO membranes is the binding of the following two conditions: (i) increasing retention for selective solutes through sieving using small nanometer-sized apertures, and (ii) allowing large flux of pure water to pass through the membrane by increasing membrane permeability. These two challenges can be visualized through a schematic [170] in Figure 10 that shows the nanosized pores of the graphene membrane and water molecules with salt are comparable to each other.Figure 10
Schematic of RO membrane with nanopores [170].Finally, a recent study by Werber et al. [171] questioned and critiqued the coupled-importance of a high selective-solute retention, and high water permeability in the desalination performance. They suggest that increased water permeability seems to have little impact on the efficiency and that increased retention of selective solutes should be the most important measure for desalination efficiency. The water/salt selectivity tradeoff relationship is based on empirical relationship [172] that is still continued to be used [173] because of incomplete theoretical understanding of many physical/chemical processes that occur during desalination through membrane-based RO.(2) Governing Equations. Most of the studies use the simple convective flux and diffusive flux to model the transport of water and solute, respectively, through the nanoscale membranes as reviewed by Werber et al. [174]:(48)Jw=AΔP−Δπm,Js=BΔcm,where A=PwVw/δmRgT and B=Ps/δm. Here, AandB are the permeability coefficients for water and solute, respectively, PwandPs are the diffusive permeability for water and solute, respectively, Vw is the molar volume of water, Rg is the gas constant, T is the absolute temperature, ΔP is the applied hydraulic pressure difference across the membrane, Δπm is the osmotic pressure difference across the active layer, and δm is the active layer thickness of the membrane. Werber et al. [174] suggest that the coefficients A and B together define the selective layer performance of nonporous membranes.


#### 4.2.3. Critical Gaps and Proposed Approach

Despite the rapid progress made in developing nanomaterials for membranes, very little research has been performed in modeling the flow of water and solute retention through graphene-based membranes. Based on above discussions of liquid flow at nanoscale, water and solute transport in graphene-based membranes should be a strong function of wettability and slip length. The strong effect of membrane wettability and slip length on water permeability and volumetric flux was shown recently by Shahbabaei et al. [175] and Chen et al. [176], respectively, through molecular dynamics simulations. Werber et al. [174] reviewed important polymeric membrane materials as a measure of controlling the wettability from an experimental perspective. Therefore, in this regard it is important to understand how wettability may affect the slip length and a method to model the effect of their coupled behavior on flux.(1) Coupled Effect of Slip Length and Wettability. Bakli and Chakraborty (2012) presented a constitutive model to describe the liquid filling of water in nanopores as a function of the dynamic slip length that changes with the contact angle and the fluid acceleration due to pressure gradient as follows:(49)Ls=A1a1+cosθ2−A2a,a=γcosθρRL,a∗=γcosθηLdL/dt,where Ls is the slip length, L is the water filling length at time t, N is the parameter representing the surface roughness (lower value represents higher roughness), a is the imbibition driving acceleration parameter, a∗ is the dimensionless form of a, and A1a and A2a are functions of a.Parobek and Liu [177] reviewed the wettability of graphene and found that the intrinsic wettability of graphene based on advancing contact angle measurements suggests that graphene is a nonwetting (hydrophobic) material and its wettability is independent of the substrate underneath graphene and the number of graphene layers. Based on our earlier discussions, it is known that slip lengths on the order of tens of nanometers are observed on hydrophobic (nonwetting) surfaces that can go up to hundreds of nanometers [178, 179] for super-hydrophobic surfaces. Therefore, the hydrophobic nature of graphene further emphasizes the need to incorporate the slip length in modeling the flux of water, which can be obtained using the Hagen-Poiseuille (HP) equation for cylindrical-shaped capillary or using the HP equation modified for geometric corrections and slip flow through capillaries of different shapes:(50)q=cπr4+4r3Ls8μ.1L0,where c is a geometry correction factor [180] to account for different shapes and eccentricity of capillaries (k=1,1.43,and1.98 for capillaries with circular, square, and equilateral triangle cross-sections, respectively), r is the pore radius, μ is the water viscosity, Ls is the slip length, and L0 is the capillary length within the membrane.


### 4.3. Shale Reservoirs

#### 4.3.1. Introduction

Shales are unconventional rocks that are both the source and the reservoir for oil and gas resources. Economical production from shale formations requires creating hydraulic fractures by pumping large volume of water (∼2–6 million gallons) at high pressure. Shale rocks are predominantly composed of consolidated clay-sized particles with a high organic content [181, 182] that result in extremely small pore sizes with 80% of pores in the range of 2 nm to 15 nm [183, 184]. Because of their extremely small pore sizes, shale rocks present an opportunity to apply the fundamentals of nanoscale fluid flow physics over existing conventional theories of porous media fluid flow.


#### 4.3.2. Literature Review

The typical characteristics of a shale reservoir include low porosity, low permeability [185, 186], complex network of matrix-fracture systems [187–191], and heterogeneous mineralogy [181, 182, 192–195], and the combination of these characteristics presents a unique challenge in understanding the fluid flow and fluid-rock and fluid-fluid interactions in these reservoirs. Following theoretical considerations, some studies explored the fluid flow in these reservoirs as a function of convection, slip flow, Knudsen diffusion [196–201], surface diffusion/sorption [202–210], rock-mechanics [189, 211–215], and osmosis [216–218]. Unlike conventional reservoir rocks, shales are typically oil-wet or mixed-wet as discussed in a review by Singh [218] for various shale rocks across the world [219–222]. The wettability of shales cannot be characterized by the macroscopic method of measuring contact angle because of the nontrivial effect of line tension in nanopores, and to account for this effect, a modified contact angle method was discussed in a review by Singh [218].Other areas that have recently caught interest due to their importance in understanding the fluid behavior in shale or increasing the productivity are (i) spontaneous imbibition in shales [218, 223–225], (ii) thermodynamic phase behavior of fluids inside shales as a function of pressure drawdown during production [226–231], (iii) alternate fracturing fluids [232–237] for better reservoir stimulation and to avoid water usage, and (iv) enhancing oil recovery by injection of fluids [238–240].(1) Key Challenge. The knowledge about fluid flow in shale reservoirs is still at its infancy compared to our knowledge of fluid flow in conventional oil and gas reservoirs that has been rigorously tested over several decades. The key challenge in production from shale reservoirs is developing verifiable theories that describe flow and mechanics from these extremely tigh pore spaces, which are heterogeneous in terms of the mineralogical composition and pore structure.


#### 4.3.3. Critical Gaps and Proposed Approach

Even though the importance of shale reservoirs as a source of hydrocarbon production continues to increase with time, the production of fluids from these reservoirs remain poorly understood at the pore level and that drives many researchers to look at these problems from continuum perspective at the reservoir scale [241]. From a theoretical perspective, there are many transport mechanisms that seem to be contributing to the hydrocarbon production from the shale rocks. A better understanding to characterize those mechanisms in terms of their contribution to the flow, the period for which they remain activated, and their likely place of origin in the reservoir can be useful to optimize the production from these reservoirs. In this regard, we can elaborate this into three specific questions: (i) What is the contribution of each transport mechanism in the total production? (ii) At what point of time and for how long during a reservoir’s life do these mechanisms remain activated? (iii) What is the contribution of different sections of reservoir, such as matrix, natural fractures, and hydraulic fractures, in the transport of hydrocarbons?(1) Reservoir-on-Chip Approach. It is clear from the discussions in the literature that ultra-tight nanopore structure of the shales plays a major role in controlling the transport of fluids and their production, but despite this consensus most of the experimental studies focus on investigating the fluid transport from a continuum perspective using core-scale experiments [219, 242–246]. The core-scale investigations can be employed to explore the flow due to interaction between matrix and fractures; however, the appropriate way to investigate nanoscale transport mechanisms would be through the nanofluidics (e.g., chip-based models) framework. Recently, there has been some effort [159, 226, 247–252] in investigating the transport from nanopores using some innovative ideas that mimic the heterogeneous nanopore and fracture network through the reservoir-on-chip approach as illustrated in Figure 11. A more realistic network of matrix and fracture mimicking actual core images can be developed synthetically using micro-CT images, characterizing the fracture and matrix system in those images using an image processing software and fabricating the processed images using standard lithography techniques [247, 253]. The features of fracture network from the CT images can be patterned and etched using deep reactive ion etching. An innovative and sophisticated reservoir-on-chip model was developed by Gerami et al. [247] to study coal reservoirs, and a similar approach can be adopted to develop a reservoir-on-chip model for shale rocks.Figure 11
Schematic illustrating (a) a reservoir-on-chip framework and (b) a micromodel with a network of fractures and heterogeneous pores (inset) (adapted after [252] and [248], respectively).(2) Contribution of Liquid Slip to Total Flow. Discovery of slip flow’s contribution in enhancing performance in other micro- and nanoscale applications, for instance, energy conversion devices [143, 148], has been acknowledged with great enthusiasm by others [153] in their field. However, similar interest has not been visible in the research community studying flow from shale reservoirs, despite some theoretical studies that illustrated the importance of gas slip in shale [196, 206, 254]. Until now, there have been no experimental studies to investigate the theoretical contribution of gas slip to the total flux. Theoretical studies considering slip flow of oil and its contribution to total flux are still evading. Perhaps, the major reason why micro- and nanoscale transport mechanisms such as slip flow have not been experimentally investigated by shale research community is because of the general inclination to perform experiments using continuum scale cores that are not suitable to investigate the fluid slip. Therefore, it is suggested that to study the micro- and nanoscale transport mechanisms in shale rocks, an emphasis must be given to employ the reservoir-on-chip approach discussed earlier. Following that, a systematic approach can be adopted to study the effect of gas and liquid slip on the total flow. As a starting point, well-established theories on liquid slip can be used to verify their applicability in shales; for instance, it can be investigated whether slip length (Ls) for oil varies as Ls∝1+cosθ−2 in synthetically fabricated chips mimicking shale rock structure by varying its wettability (θ). If the surface of the fabricated chip has been characterized with some roughness function fx,y, then a modified expression for slip length accounting for roughness [128] can be used as follows:(51)Ls=k1+cosθ21+α∗2,α∗=∑fx,y−fx,y¯N,where k is the constant of proportionality and α∗ is a parameter that characterizes the roughness of a surface through its amplitude N and roughness function fx,y. Now, the flow rate as a function of slip length would be given as follows:(52)Qslip=Qno−slip1+6Ls2h,where h is the distance between the confining pore walls and Qslip and Qno−slip are the flow rates due to slip and no-slip boundary conditions, respectively.(3) Contribution of Evaporation to Water Loss. Hydraulic fracturing requires injecting large volume of water (∼2–6 million gallons); however, on average, only 6–10% of the injected water is recovered in the US across all shale plays [255]. This abundant retention of fracturing fluid inside the shale formations is a cause of major concern because it keeps the hydrocarbons from flowing out of the reservoir by reducing the relative permeability of the hydrocarbons inside the formation. Singh [218] reviewed the process of water uptake by shales and suggested evaporation as one of the possible mechanisms that may be contributing to loss of water in shales. To support their hypothesis, they provided three characteristics associated with shale formations that could facilitate evaporation of water: (1) gas expansion, (2) water droplets formed at asperities or cavities in fractures, and (3) pore confinement. Since the hypothesis of water loss in shales due to evaporation has not been investigated before, it can be tested using the reservoir-on-chip approach with the network of small pores and fractures. Among the three characteristics noted above that could be facilitating evaporation, gas expansion seems to be the most intuitive source and can be mimicked by injecting gas into the water saturated chip with a network of matrix pores and fractures under reservoir conditions. The process of fabricating such a chip and other experimental devices needed to test this process can be set up based on the reservoir-on-chip model developed by Gerami et al. [247] to study coal reservoirs and shown in Figure 12.Figure 12
Reservoir-on-chip model used to depict the fracture structure in a coal reservoir [247].Keeping track of water saturation and its phase change with time can indicate whether expansion of gas as it moves from pores to fractures results in evaporating some water. The phase change can be studied using direct visual observation of water phase behavior due to light scattering, a method that has been reliably tested for vapor detection by Yang et al. [256] and more recently in the context of shales by Ally et al. [227] to study the fluid phase behavior in nanopores.


### 4.4. Hydrogen Storage

#### 4.4.1. Introduction

Hydrogen-based energy is one of the promising sources of energy [257] because hydrogen is a zero-emission fuel and can provide solutions to many versatile problems such as predictable supply that is not dependent on weather conditions, its capability to integrate within existing infrastructure developed for other forms of energy, and its applications in many different industries [258, 259]. As such, hydrogen-based energy solutions have versatile applications, such as sources of power [260] (fuel cells [261]) in transportation [262, 263] and as a commodity in many industrial processes [264]. Figure 13 shows a general illustration of a hydrogen energy system with sources, production, storage, and its applications. Among all sources of hydrogen production worldwide, 95% of the hydrogen production comes from nonrenewable sources [265] such as natural gas, coal, gasoline, diesel fuels, etc.Figure 13
Schematic representing a hydrogen energy system. About 95% of the hydrogen currently produced across the world comes from nonrenewable sources.Among various steps in the hydrogen energy system, hydrogen storage is a major challenge and a key enabling technology for the advancement of the hydrogen energy-based applications.


#### 4.4.2. Literature Review

Hydrogen can be stored either physically (as compressed gas or liquid) or chemically (on the surfaces of solid materials by adsorption or within the solid materials by absorption). Exploiting hydrogen energy on commercial scale as planned by the US Department of Energy through portable power sources, fuel cells, and transportation requires a compact system that can store hydrogen on these applications. The fundamental difficulty in achieving this task is low energy density per unit volume of hydrogen, which means hydrogen needs relatively larger volume to pack same amount of energy compared to other fuels. The storage of hydrogen can be divided in three sectors [266–269], namely: (i) high-pressured gaseous state, (ii) cryogenic liquid state, and (iii) solid state. Because of the low storage efficiency of liquid and gaseous hydrogen, research efforts to develop storage technologies have focused towards increasing the density of hydrogen either through chemical storage technology using new materials [270–273] or through physical storage using compression [274–276].Even though some researchers [271] suggest that the chemical method is more efficient in storing hydrogen than the physical method, they do not provide a comparison or a reference to support such a claim. However, Petitpas et al. [275] provided an objective comparison between the two storage methods, i.e., physical storage and chemical storage, which suggested that at low pressures the density of stored hydrogen through the chemical method is significantly higher than the density of hydrogen stored through the physical method. As the pressure is increased, the density of hydrogen through the chemical method starts to flatten, whereas the density of hydrogen through the physical method continues to increase, eventually overtaking the density of hydrogen stored through the chemical method at a breakeven pressure beyond which the density of hydrogen stored through the physical method remains larger than the hydrogen stored through the chemical method (Figure 14). They reported that this breakeven pressure strongly depends on temperature.Figure 14
Comparison of physical and chemical storage capacities [275].Since the choice of hydrogen storage technology would depend on the end application, Petitpas et al. [275] suggested that the conditions relevant to a specific application may supersede the chemistry and physics required to enhance the hydrogen density. For instance, the added mass of the adsorbent for chemical storage may discourage the use of this storage technology in the case of weight-sensitive applications, whereas in the case of storage by compression the use of higher pressures may discourage its use for some applications where high pressure is a safety concern. Additionally, the cryogenic temperature needed for both techniques could be discouraging in the absence of liquid hydrogen (or liquid nitrogen for 80 K fueling). Petitpas et al. [275] also compared the use of a “hybrid” system that can store hydrogen by both compression and adsorbent and found that they can store 10% to 20% (with metal-organic framework as the sorbent) more hydrogen compared to storage by compression alone.Despite the variation in density of stored hydrogen by the two approaches, the research on hydrogen storage in past few years has mostly focused on chemical storage by adsorbents because of its expected use in many commercial and portable applications that may have safety concerns related to high pressure. Therefore, below we review the modeling scheme to store hydrogen by adsorption on porous and deformable materials.(1) Key Challenge. Increasing the efficiency of hydrogen storage is one of the major challenges in the application of hydrogen for consumer applications such as transportation. Figure 15 shows a schematic comparing the volumetric scale of storage space required for 4 kg of hydrogen storage in different materials and storage methods.Figure 15
Schematic representation of the volumetric scale of storage space required for 4 kg of hydrogen with different modes of storage [259, 268]. Car image is a courtesy of Toyota press information, 33rd Tokyo Motor Show.(2) Governing Equations. The governing equations for the storage of hydrogen follow the conservation of mass and energy.Mass balance:(53)∂ϕρg∂t+∇.ϕρgu=−∂sm∂t,∂sm∂t=∂1−ϕρs∂t=∂ρs∂t−∂ϕρs∂t,where ρg is the hydrogen gas density, ϕ is the porosity of the adsorbent, u is the hydrogen gas velocity that is computed using either Navier–Stokes equation or Darcy’s law, sm is the mass of hydrogen adsorbed or desorbed per unit volume of the sorbent, and ρs is the density of the sorbent.Energy balance:(54)∂ρcp¯T∂t+∇.ρgcp,guT=∇.keff∇T+∂sT∂t,where ρcp¯ is the effective heat capacity, keff is the effective thermal conductivity [277], T is the temperature, and sT is the product of the enthalpy change (as a result of sorption).


#### 4.4.3. Critical Gaps and Proposed Approach

Even though the experimental research on developing superior nanostructured adsorbents [271, 273, 278, 279] for hydrogen storage continues to increase with time, the same effort has not been devoted to model the transport of hydrogen in these nanostructured porous materials. For instance, transport of hydrogen is modeled using only the continuum scale transport [277, 280–283] with no contribution of nanoscale transport mechanisms despite the role of nanoconfinement [273, 284] as an effective strategy in modifying the hydrogen storage performance. From a theoretical perspective, transport of hydrogen in nanoporous mediums is controlled by two more mechanisms (Knudsen diffusion and slip flow) other than viscous/convective flow, heat transfer, and sorption (surface diffusion). A better understanding to characterize the transport mechanisms at nanoscale in terms of their contribution to the storage and the operating conditions for which they remain activated can be useful to optimize the storage of hydrogen in nanostructured materials. In this regard, we discuss two transport mechanisms of Knudsen diffusion and slip flow at nanoscale and how they can be added in the governing equations to model hydrogen storage.(1) Accounting for Slip Flow and Knudsen Diffusion in Total Flux. The flow behavior at the continuum regime is governed by the Navier–Stokes equations; however, as Kn increases beyond the range of the continuum regime, the no-slip boundary conditions start to break such that the N-S equations should be modified to account for slip boundary conditions or molecular models should be employed. Another transport mechanism that plays an important role when the mean free path of the gas molecules is larger than or equal to the pore size is the form of diffusive transport called Knudsen diffusion. Diffusive transport in micro- and nanopores can occur through three different mechanisms: (i) bulk diffusion, (ii) surface diffusion (sorption), and (iii) Knudsen diffusion [67]. The proportion of each mechanism depends on the magnitude of the Knudsen number Kn, and for large Kn, the molecules interact more often with the pore walls than with each other. Typical nanoporous materials for storing hydrogen have pore size that can vary from less than 1 nm to about 10 nm [270, 271, 273, 285, 286] for varied materials ranging from carbons, zeolites, metal-organic framework (MOF), complex hydrides, etc., whereas mean free path of hydrogen molecules at ambient conditions estimated from kinetic theory is about 160 nm. Therefore, it is evident from this comparison that Knudsen diffusion may have an important role to play in the transport of hydrogen.(i) Slip Flow and Knudsen Diffusion: The convective mass flow in the presence of slip boundary conditions and diffusive mass flow in the presence of Knudsen diffusion is given as follows [69]:(55)m˙C=MuARTp0l,m˙D=DKnMARTΔpl=13s8RTπMMARTΔpl,s=4hπ,channel  shaped  pores,2r,tube  shaped  pores.The solution of the Navier–Stokes equation and slip flow boundary conditions are used to solve for hydrogen velocity due to convective flux for two pore shapes of channels and tubes, respectively, as follows:(56)u=2h∫−h/2h/212μ∂p∂xy2+c1y+c2dy,channel‐shaped  pores,1πr∗2∫0r∗14μ∂p∂xr2+c3lnr+c42πrdr,tube‐shaped  pores,uslip=2−σuσuλ∂u∂yy=h,channel‐shaped  pores,2−σuσuλ∂u∂rr=r∗,tube‐shaped  pores,where m˙D is the diffusive mass flow, DKn is the Knudsen diffusion coefficient, s is the characteristic length, p is the pressure, u is the hydrogen velocity due to convective flux, σu is TMAC, λ is the mean free path, p0 is the outlet pressure, M is the molar mass, A is the cross-sectional area of the medium, h is the height of the channel-shaped pore, r∗ is the radius of a cylindrical-shaped pore, l is the length of the medium, R is the universal gas constant, T is the temperature, and c1,c2,c3,andc4 are constants from integration of the Navier–Stokes equation.The total flux is the sum of the mass flow due to convective flux (including slip flow) and diffusive flux (from Knudsen diffusion), eventually leading to an effective velocity of hydrogen molecules (ueff) that can be substituted in place of u in the mass balance equation for hydrogen:(57)m˙=m˙C+m˙D,⇒m˙=MuARTp0l+DKnMARTΔpl,ueff=m˙ρgA,ueff=MρgRTlup0+DKnΔp.


## 4.1. Energy Conversion through Electrokinetics

### 4.1.1. Introduction

Increasing demand for energy has led to harnessing of energy from various forms such as mechanical, chemical, internal, electrical, thermal, etc. and various other combinations. To utilize these various sources of energy, micro- and nanofluidics provide the capability of converting energy from one form to another form to cater to different applications. Energy conversion is possible through electrokinetics, which is a term used to describe different transport mechanisms in the presence of an electric field and liquid flow through a medium. Electrokinetics can be used for energy conversion in four different ways: (i) drive liquid by applying an electric potential (electro-osmosis), (ii) transport and separate particles by applying an electric potential (electrophoresis), (iii) generate electric potential by pressure gradient-driven liquid flow (streaming potential), and (iv) generate electric potential by movement of charged particles under gravitational or centrifugal force (sedimentation potential). In particular, micro- and nanofluidics have found application in converting energy from heat, electrical, chemical/biochemical, pressure, etc. to electrical energy using micro- and nanoscale technologies with a wide range of efficiencies varying from 3% for streaming potential devices [131] to 60% for micro-fuel cells [132].


### 4.1.2. Literature Review

Electric double layer (∼1 nm thickness) is a common phenomenon in small energy conversion devices that are increasingly based on nanostructured materials. As an example to illustrate the scale at which fluids flow in electrokinetic devices, Figure7 shows transmission electron microscope images of carbon composites [133, 134] with application in electrochemical energy conversion.Figure 7
Transmission electron microscope images of carbon composites for application in electrochemical energy conversion [133, 134].The electrokinetic effects at confined scales result in enhanced solid-fluid interaction close to the surface, which may also influence the liquid slip flow and consequently the overall fluid transport. The amount of liquid slip on a surface is found to be affected by the presence of a charged surface and/or liquids with charged ions such as electrolytes. Two examples of electrokinetic flows relevant here are electro-osmosis (liquid flow induced by application of an electric field) and streaming currents (electric current induced by application of a pressure gradient).(1) Key Challenge. One of the key challenges in the design of electrokinetic energy conversion devices is their low conversion efficiencies. The maximum efficiency of energy conversion ηmax, which is defined as the ratio of the output power to the input power, is given as follows:(40)ηmax=PoweroutputPowerinput.(2) Governing Equations. The mathematical theory used to model the transport in the presence of chemical potential (due to surface charge and charged ions in the solution), externally applied electric potential, and hydrostatic potential is briefly reviewed below. The chemical potential and ionic concentrations in the diffuse part of the EDL are governed by the Poisson-Boltzmann equation as follows:(41)∇2ψy=−1ερy=−eε∑i=1Nziciexp−eziψykBT,where ψy is the potential in the EDL due to surface charge and the charged ions in the solution, y is the radial distance from the surface of the channel, ρy is the volume charge density of all the ions present in the neighborhood of the surface, e is the unit charge, ε is the permittivity of the liquid, ci and zi are the concentration and the ionic valence of the ith ionic species, respectively, kB is the Boltzmann constant, and T is the absolute temperature.The Poisson-Boltzmann equation does not have a general analytical solution, but its solutions are available for specific cases [135]. For our interest, we consider the solution for a flat surface with a low potential (ψsurface) at its surface for which the above equation can be linearized based on the Debye-Hückel approximation to obtain the following solution:(42)ψy=ψsurfaceexp−κy,λD=κ−1=εkBTNAe2∑i=1Nzi2ci,where κ−1 is called the Debye length λD that gives a measure of the EDL and NA is the Avogadro number. The Debye length is independent of the charge on the solid surface and increases as the solute ion concentration in the solution becomes dilute such that λD for a solution without electrolytes can be considered as infinitely thick. Typical values of λD can vary from ∼3nm to ∼300nm for electrolytes with ionic concentration of 10−2M and 10−6M, respectively [135].The above-defined definitions are finally used to describe the solvent and solute transport in the presence of hydrostatic, chemical, and electric potential within a micro- or a nanochannel of width2h as follows.(i) Solvent Transport. The solvent flux across the channel due to hydrostatic potential and electric potential is the sum of the average of the fluxes due to these two potentials and given as follows [136, 137]:(43)vs=1h∫0hvsydy,vs=1h∫0h12ηy2−2hydPdx︸flow  due  tohydrostatic  potentialdy+1h∫0hεηζ−ψydϕdx︸flow  due  toelectric  potentialdy,where vs is the average solvent velocity, ζ is the electric potential at the slipping plane (also called zeta potential), ψy is the potential in the EDL, η is the solution viscosity, ε is the permittivity of the solution, and ϕ is an externally applied electric potential.(ii) Solute Transport. The solute flux across the channel is given as the averaged sum of the flux due to hydrostatic potential, flux due to diffusion, and flux due to electric potential. The solute transport in the presence of these three mechanisms is also called the Nernst-Planck equation and given as follows [137–139]:(44)J=1h∫0hcyvsy︸flux  due  tohydrostatic  potentialdy−1h∫0hD∞dcydx︸flux  due  todiffusiondy−1h∫0hzeD∞kBTcydϕdx︸flux  due  toelectric  potentialdy,where J is the area-averaged solute flux, D∞ is the diffusion coefficient of the solute in the free solution, cy is the solute concentration at distance y from the surface, and ze is the ionic charge.(iii) Solute Partitioning. Solute ions experience electric, chemical, and surface forces that vary in magnitude based on the size and charge on the ions; therefore, this leads to nonuniform distribution of ions along the width of the flowing medium. This partitioning of the solute ions can be characterized through the ratio of the average solute concentration c to the bulk concentration c0 in the solution, which is also referred to as the partitioning coefficient βpart [140]:(45)βpart=cc0=1h∫∼0h−rexp−EykBTdy,where βpart is the solute partition coefficient that is the measure of ion permeability-selectivity of the pore, r is the radius of the spherical solute, and E is the electrostatic interaction potential. It should be noted that when the confinement width becomes comparable to the Debye length, the overlapped EDL gets deformed significantly due to steric effects and in that case the modified Poisson-Boltzmann equation must be employed [141]. When the transport of solute through the medium is strongly dependent on solute size and independent of frictional and hydrostatic (convective) forces, in that case [137, 142],(46)βpart=1−rh.(iv) Conversion Efficiency. The maximum efficiency of energy conversion ηmax is defined as the ratio of the output power to the input power; however, some researchers [143–145] also use an alternate measure for the conversion efficiency using a measure called electrokinetic figure of merit [146], which is based on phenomenological transport equations. The following equation is an alternate form for the maximum efficiency in terms of figure of merit β:(47)ηmax=1+β−11+β+1,where β=ν2σ/kH. Here, ν, σ, and kH are the streaming potential coefficient, the ion conductivity, and the hydraulic permeability, respectively.


### 4.1.3. Critical Gaps and Proposed Approach

About a decade ago it was hypothesized [147, 148] that the efficiencies of the energy conversion devices based on electrokinetic effects can be optimized by enhancing the liquid slip. A theoretical study by Davidson and Xuan [143] showed that the efficiency of electrokinetic energy conversion devices is almost a linear function of slip length for small values of zeta potentials. Even though there is a consensus that liquid slip plays a strong role in the device efficiency among the electrokinetic energy research community, some issues remain in regard to understanding of slip length as a function of key device parameters and its relationship with flux.(1) Slip Length as a Function of Device Parameters. In the electrokinetic energy research community, performance of the device is predicted by either completely neglecting the slip [149, 150] or an assumed value of slip length used by others [143, 148, 151, 152] that is independent of the surface and flow conditions. Further, some researchers [153, 154] suggested that the slip length increases as the surface becomes slippery, or in other words slip length increases on surfaces with low friction coefficient. However, this suggestion is contrary to the observations reported by researchers [155–158] in the field of liquid slip who reported that the surface friction results in increase in the effective viscosity that consequently increases the liquid slip length [156]. This result, even though counterintuitive in nature, is supported by experimental observations [90, 96, 99, 159]. Therefore, as noted by others [160], there is a lack of understanding about the liquid slip length when it comes to predicting the conversion efficiencies. Developing the fundamental understanding of liquid slip as a function of key device parameters presents an opportunity of increasing the efficiency of energy conversion devices by exploiting their potential of developing large slip lengths. For instance, we discussed many slip length models earlier as a function of surface wettability, surface friction, shear rate, surface roughness, and electrokinetic effects that can be used to either estimate the slip lengths for a given device or design the device parameters for a given slip length.(2) Slip Length vs. Flux. Following theoretical considerations initiated by some studies [143, 147, 148], it is now believed by most researchers that the efficiency of energy conversion devices can be increased by increasing the liquid slip length. In this regard, some comments were made on how to optimally exploit this potential. Having said that, the arguments made in favor of slip are theoretical in nature, and there are no experimental studies specific to energy conversion devices that verify the amplified values of conversion efficiencies predicted theoretically and usually cited by many authors [137, 153]. The objective of this note is to caution that a theoretical definition of liquid slip where flux is proportional to the slip length (vx∝1+6Ls/h) may not be universally true. To support this hypothesis, we refer to the experimental study of Liu and Li [75] who measured the flux and slip length for liquid argon and liquid helium in nanochannels of various widths and noted that the relation of slip length to flux was not clear from the data. They attributed this reason to nonuniform distribution of density and viscosity of the confined fluid.


## 4.1.1. Introduction

Increasing demand for energy has led to harnessing of energy from various forms such as mechanical, chemical, internal, electrical, thermal, etc. and various other combinations. To utilize these various sources of energy, micro- and nanofluidics provide the capability of converting energy from one form to another form to cater to different applications. Energy conversion is possible through electrokinetics, which is a term used to describe different transport mechanisms in the presence of an electric field and liquid flow through a medium. Electrokinetics can be used for energy conversion in four different ways: (i) drive liquid by applying an electric potential (electro-osmosis), (ii) transport and separate particles by applying an electric potential (electrophoresis), (iii) generate electric potential by pressure gradient-driven liquid flow (streaming potential), and (iv) generate electric potential by movement of charged particles under gravitational or centrifugal force (sedimentation potential). In particular, micro- and nanofluidics have found application in converting energy from heat, electrical, chemical/biochemical, pressure, etc. to electrical energy using micro- and nanoscale technologies with a wide range of efficiencies varying from 3% for streaming potential devices [131] to 60% for micro-fuel cells [132].


## 4.1.2. Literature Review

Electric double layer (∼1 nm thickness) is a common phenomenon in small energy conversion devices that are increasingly based on nanostructured materials. As an example to illustrate the scale at which fluids flow in electrokinetic devices, Figure7 shows transmission electron microscope images of carbon composites [133, 134] with application in electrochemical energy conversion.Figure 7
Transmission electron microscope images of carbon composites for application in electrochemical energy conversion [133, 134].The electrokinetic effects at confined scales result in enhanced solid-fluid interaction close to the surface, which may also influence the liquid slip flow and consequently the overall fluid transport. The amount of liquid slip on a surface is found to be affected by the presence of a charged surface and/or liquids with charged ions such as electrolytes. Two examples of electrokinetic flows relevant here are electro-osmosis (liquid flow induced by application of an electric field) and streaming currents (electric current induced by application of a pressure gradient).(1) Key Challenge. One of the key challenges in the design of electrokinetic energy conversion devices is their low conversion efficiencies. The maximum efficiency of energy conversion ηmax, which is defined as the ratio of the output power to the input power, is given as follows:(40)ηmax=PoweroutputPowerinput.(2) Governing Equations. The mathematical theory used to model the transport in the presence of chemical potential (due to surface charge and charged ions in the solution), externally applied electric potential, and hydrostatic potential is briefly reviewed below. The chemical potential and ionic concentrations in the diffuse part of the EDL are governed by the Poisson-Boltzmann equation as follows:(41)∇2ψy=−1ερy=−eε∑i=1Nziciexp−eziψykBT,where ψy is the potential in the EDL due to surface charge and the charged ions in the solution, y is the radial distance from the surface of the channel, ρy is the volume charge density of all the ions present in the neighborhood of the surface, e is the unit charge, ε is the permittivity of the liquid, ci and zi are the concentration and the ionic valence of the ith ionic species, respectively, kB is the Boltzmann constant, and T is the absolute temperature.The Poisson-Boltzmann equation does not have a general analytical solution, but its solutions are available for specific cases [135]. For our interest, we consider the solution for a flat surface with a low potential (ψsurface) at its surface for which the above equation can be linearized based on the Debye-Hückel approximation to obtain the following solution:(42)ψy=ψsurfaceexp−κy,λD=κ−1=εkBTNAe2∑i=1Nzi2ci,where κ−1 is called the Debye length λD that gives a measure of the EDL and NA is the Avogadro number. The Debye length is independent of the charge on the solid surface and increases as the solute ion concentration in the solution becomes dilute such that λD for a solution without electrolytes can be considered as infinitely thick. Typical values of λD can vary from ∼3nm to ∼300nm for electrolytes with ionic concentration of 10−2M and 10−6M, respectively [135].The above-defined definitions are finally used to describe the solvent and solute transport in the presence of hydrostatic, chemical, and electric potential within a micro- or a nanochannel of width2h as follows.(i) Solvent Transport. The solvent flux across the channel due to hydrostatic potential and electric potential is the sum of the average of the fluxes due to these two potentials and given as follows [136, 137]:(43)vs=1h∫0hvsydy,vs=1h∫0h12ηy2−2hydPdx︸flow  due  tohydrostatic  potentialdy+1h∫0hεηζ−ψydϕdx︸flow  due  toelectric  potentialdy,where vs is the average solvent velocity, ζ is the electric potential at the slipping plane (also called zeta potential), ψy is the potential in the EDL, η is the solution viscosity, ε is the permittivity of the solution, and ϕ is an externally applied electric potential.(ii) Solute Transport. The solute flux across the channel is given as the averaged sum of the flux due to hydrostatic potential, flux due to diffusion, and flux due to electric potential. The solute transport in the presence of these three mechanisms is also called the Nernst-Planck equation and given as follows [137–139]:(44)J=1h∫0hcyvsy︸flux  due  tohydrostatic  potentialdy−1h∫0hD∞dcydx︸flux  due  todiffusiondy−1h∫0hzeD∞kBTcydϕdx︸flux  due  toelectric  potentialdy,where J is the area-averaged solute flux, D∞ is the diffusion coefficient of the solute in the free solution, cy is the solute concentration at distance y from the surface, and ze is the ionic charge.(iii) Solute Partitioning. Solute ions experience electric, chemical, and surface forces that vary in magnitude based on the size and charge on the ions; therefore, this leads to nonuniform distribution of ions along the width of the flowing medium. This partitioning of the solute ions can be characterized through the ratio of the average solute concentration c to the bulk concentration c0 in the solution, which is also referred to as the partitioning coefficient βpart [140]:(45)βpart=cc0=1h∫∼0h−rexp−EykBTdy,where βpart is the solute partition coefficient that is the measure of ion permeability-selectivity of the pore, r is the radius of the spherical solute, and E is the electrostatic interaction potential. It should be noted that when the confinement width becomes comparable to the Debye length, the overlapped EDL gets deformed significantly due to steric effects and in that case the modified Poisson-Boltzmann equation must be employed [141]. When the transport of solute through the medium is strongly dependent on solute size and independent of frictional and hydrostatic (convective) forces, in that case [137, 142],(46)βpart=1−rh.(iv) Conversion Efficiency. The maximum efficiency of energy conversion ηmax is defined as the ratio of the output power to the input power; however, some researchers [143–145] also use an alternate measure for the conversion efficiency using a measure called electrokinetic figure of merit [146], which is based on phenomenological transport equations. The following equation is an alternate form for the maximum efficiency in terms of figure of merit β:(47)ηmax=1+β−11+β+1,where β=ν2σ/kH. Here, ν, σ, and kH are the streaming potential coefficient, the ion conductivity, and the hydraulic permeability, respectively.


## 4.1.3. Critical Gaps and Proposed Approach

About a decade ago it was hypothesized [147, 148] that the efficiencies of the energy conversion devices based on electrokinetic effects can be optimized by enhancing the liquid slip. A theoretical study by Davidson and Xuan [143] showed that the efficiency of electrokinetic energy conversion devices is almost a linear function of slip length for small values of zeta potentials. Even though there is a consensus that liquid slip plays a strong role in the device efficiency among the electrokinetic energy research community, some issues remain in regard to understanding of slip length as a function of key device parameters and its relationship with flux.(1) Slip Length as a Function of Device Parameters. In the electrokinetic energy research community, performance of the device is predicted by either completely neglecting the slip [149, 150] or an assumed value of slip length used by others [143, 148, 151, 152] that is independent of the surface and flow conditions. Further, some researchers [153, 154] suggested that the slip length increases as the surface becomes slippery, or in other words slip length increases on surfaces with low friction coefficient. However, this suggestion is contrary to the observations reported by researchers [155–158] in the field of liquid slip who reported that the surface friction results in increase in the effective viscosity that consequently increases the liquid slip length [156]. This result, even though counterintuitive in nature, is supported by experimental observations [90, 96, 99, 159]. Therefore, as noted by others [160], there is a lack of understanding about the liquid slip length when it comes to predicting the conversion efficiencies. Developing the fundamental understanding of liquid slip as a function of key device parameters presents an opportunity of increasing the efficiency of energy conversion devices by exploiting their potential of developing large slip lengths. For instance, we discussed many slip length models earlier as a function of surface wettability, surface friction, shear rate, surface roughness, and electrokinetic effects that can be used to either estimate the slip lengths for a given device or design the device parameters for a given slip length.(2) Slip Length vs. Flux. Following theoretical considerations initiated by some studies [143, 147, 148], it is now believed by most researchers that the efficiency of energy conversion devices can be increased by increasing the liquid slip length. In this regard, some comments were made on how to optimally exploit this potential. Having said that, the arguments made in favor of slip are theoretical in nature, and there are no experimental studies specific to energy conversion devices that verify the amplified values of conversion efficiencies predicted theoretically and usually cited by many authors [137, 153]. The objective of this note is to caution that a theoretical definition of liquid slip where flux is proportional to the slip length (vx∝1+6Ls/h) may not be universally true. To support this hypothesis, we refer to the experimental study of Liu and Li [75] who measured the flux and slip length for liquid argon and liquid helium in nanochannels of various widths and noted that the relation of slip length to flux was not clear from the data. They attributed this reason to nonuniform distribution of density and viscosity of the confined fluid.


## 4.2. Membrane-Based Water Desalination

### 4.2.1. Introduction

Currently, there are 4 major sources of water that can be used for water desalination and treatment, i.e., produced water from oil and gas wells, brackish groundwater, municipal wastewater, and seawater, respectively. Out of these, seawater is a source with continuous and large amounts of saline water that can be potentially desalinated and supplied to various applications on large scale as illustrated in Figure8.Figure 8
Role of water in the energy sector.Desalination using the membrane-based reverse osmosis process is a potentially attractive technology that can be used to desalinate enormous amount of seawater. Reverse osmosis (RO) is considered to be a leading technology for desalination primarily because of its relatively low cost (one-half to one-third of the cost of distillation) due to relatively low energy consumption (1.8 kW·h/m3) that contributes to about 44% of production capacity in the world [161]. RO uses a semipermeable membrane to stop the contaminants at the level of 0.1 nm to 1 nm by applying hydrostatic pressure to overcome osmotic pressure, which is a chemical potential developed due to salinity difference of the solvents that drives the solvent from an aqueous solution of low salt concentration to an aqueous solution of high salt concentration. In the RO, the movement of solvent from low concentration to high concentration due to chemical potential is reversed by applying an external pressure.


### 4.2.2. Literature Review

The study by Nair et al. [162] was the first to discover that graphene oxide has the membrane properties suitable for desalination, i.e., unique water permeation pathway and retention for other gases and ions/solutes. The water flows through the multiple layers of graphene oxide through two-dimensional capillaries formed by closely spaced graphene sheets where the salt ions are retained on the sheets as shown in Figure 9.Figure 9
Graphene oxide membrane with nanopores [162].Conventional understanding of the flow fluid through pores suggests that this dual-task is paradoxical in nature, but development of nanomaterials such as carbon nanotubes (CNTs), boron nitride nanotubes (BNNTs), and graphene that show unconventional fluid transport behavior due to their smallest possible pore volumes combined with ionic dehydration [163], electrokinetic effects [164], and size exclusion [165] properties provide opportunity for their application in optimizing two coupled tasks of desalination mentioned above. Specifically, high selective-solute retention and high permeability to water is achieved by drilling sub-nanometer holes through the solid-state nanomaterials such as graphene by using focus ion beam [166, 167] or oxygen plasma etching process [168]. The water molecules are too big to pass through graphene’s fine mesh without making holes/pores in the graphene sheet and poking these sub-nanometer holes uniformly is another challenge. Additionally, in order to drive satisfactory magnitude of water flux through these nanopores, a large pressure difference is required that cannot be solely achieved by osmotic pressure. However, compared to the fine mesh of nanoporous graphene, Nair et al. [162] proposed graphene oxide membrane as an alternate to the graphene because graphene oxide relies on its naturally occurring tortuous path to flow out water (Figure 9), unlike the challenging task of artificially poking holes through graphene for creating a path to flow out water. Other than this major advantage of graphene oxide, there are several other reasons as presented by You et al. [169] that make graphene oxide a preferred choice over nanoporous graphene.(1) Key Challenge. One of the key design challenges in the desalination technology using RO membranes is the binding of the following two conditions: (i) increasing retention for selective solutes through sieving using small nanometer-sized apertures, and (ii) allowing large flux of pure water to pass through the membrane by increasing membrane permeability. These two challenges can be visualized through a schematic [170] in Figure 10 that shows the nanosized pores of the graphene membrane and water molecules with salt are comparable to each other.Figure 10
Schematic of RO membrane with nanopores [170].Finally, a recent study by Werber et al. [171] questioned and critiqued the coupled-importance of a high selective-solute retention, and high water permeability in the desalination performance. They suggest that increased water permeability seems to have little impact on the efficiency and that increased retention of selective solutes should be the most important measure for desalination efficiency. The water/salt selectivity tradeoff relationship is based on empirical relationship [172] that is still continued to be used [173] because of incomplete theoretical understanding of many physical/chemical processes that occur during desalination through membrane-based RO.(2) Governing Equations. Most of the studies use the simple convective flux and diffusive flux to model the transport of water and solute, respectively, through the nanoscale membranes as reviewed by Werber et al. [174]:(48)Jw=AΔP−Δπm,Js=BΔcm,where A=PwVw/δmRgT and B=Ps/δm. Here, AandB are the permeability coefficients for water and solute, respectively, PwandPs are the diffusive permeability for water and solute, respectively, Vw is the molar volume of water, Rg is the gas constant, T is the absolute temperature, ΔP is the applied hydraulic pressure difference across the membrane, Δπm is the osmotic pressure difference across the active layer, and δm is the active layer thickness of the membrane. Werber et al. [174] suggest that the coefficients A and B together define the selective layer performance of nonporous membranes.


### 4.2.3. Critical Gaps and Proposed Approach

Despite the rapid progress made in developing nanomaterials for membranes, very little research has been performed in modeling the flow of water and solute retention through graphene-based membranes. Based on above discussions of liquid flow at nanoscale, water and solute transport in graphene-based membranes should be a strong function of wettability and slip length. The strong effect of membrane wettability and slip length on water permeability and volumetric flux was shown recently by Shahbabaei et al. [175] and Chen et al. [176], respectively, through molecular dynamics simulations. Werber et al. [174] reviewed important polymeric membrane materials as a measure of controlling the wettability from an experimental perspective. Therefore, in this regard it is important to understand how wettability may affect the slip length and a method to model the effect of their coupled behavior on flux.(1) Coupled Effect of Slip Length and Wettability. Bakli and Chakraborty (2012) presented a constitutive model to describe the liquid filling of water in nanopores as a function of the dynamic slip length that changes with the contact angle and the fluid acceleration due to pressure gradient as follows:(49)Ls=A1a1+cosθ2−A2a,a=γcosθρRL,a∗=γcosθηLdL/dt,where Ls is the slip length, L is the water filling length at time t, N is the parameter representing the surface roughness (lower value represents higher roughness), a is the imbibition driving acceleration parameter, a∗ is the dimensionless form of a, and A1a and A2a are functions of a.Parobek and Liu [177] reviewed the wettability of graphene and found that the intrinsic wettability of graphene based on advancing contact angle measurements suggests that graphene is a nonwetting (hydrophobic) material and its wettability is independent of the substrate underneath graphene and the number of graphene layers. Based on our earlier discussions, it is known that slip lengths on the order of tens of nanometers are observed on hydrophobic (nonwetting) surfaces that can go up to hundreds of nanometers [178, 179] for super-hydrophobic surfaces. Therefore, the hydrophobic nature of graphene further emphasizes the need to incorporate the slip length in modeling the flux of water, which can be obtained using the Hagen-Poiseuille (HP) equation for cylindrical-shaped capillary or using the HP equation modified for geometric corrections and slip flow through capillaries of different shapes:(50)q=cπr4+4r3Ls8μ.1L0,where c is a geometry correction factor [180] to account for different shapes and eccentricity of capillaries (k=1,1.43,and1.98 for capillaries with circular, square, and equilateral triangle cross-sections, respectively), r is the pore radius, μ is the water viscosity, Ls is the slip length, and L0 is the capillary length within the membrane.


## 4.2.1. Introduction

Currently, there are 4 major sources of water that can be used for water desalination and treatment, i.e., produced water from oil and gas wells, brackish groundwater, municipal wastewater, and seawater, respectively. Out of these, seawater is a source with continuous and large amounts of saline water that can be potentially desalinated and supplied to various applications on large scale as illustrated in Figure8.Figure 8
Role of water in the energy sector.Desalination using the membrane-based reverse osmosis process is a potentially attractive technology that can be used to desalinate enormous amount of seawater. Reverse osmosis (RO) is considered to be a leading technology for desalination primarily because of its relatively low cost (one-half to one-third of the cost of distillation) due to relatively low energy consumption (1.8 kW·h/m3) that contributes to about 44% of production capacity in the world [161]. RO uses a semipermeable membrane to stop the contaminants at the level of 0.1 nm to 1 nm by applying hydrostatic pressure to overcome osmotic pressure, which is a chemical potential developed due to salinity difference of the solvents that drives the solvent from an aqueous solution of low salt concentration to an aqueous solution of high salt concentration. In the RO, the movement of solvent from low concentration to high concentration due to chemical potential is reversed by applying an external pressure.


## 4.2.2. Literature Review

The study by Nair et al. [162] was the first to discover that graphene oxide has the membrane properties suitable for desalination, i.e., unique water permeation pathway and retention for other gases and ions/solutes. The water flows through the multiple layers of graphene oxide through two-dimensional capillaries formed by closely spaced graphene sheets where the salt ions are retained on the sheets as shown in Figure 9.Figure 9
Graphene oxide membrane with nanopores [162].Conventional understanding of the flow fluid through pores suggests that this dual-task is paradoxical in nature, but development of nanomaterials such as carbon nanotubes (CNTs), boron nitride nanotubes (BNNTs), and graphene that show unconventional fluid transport behavior due to their smallest possible pore volumes combined with ionic dehydration [163], electrokinetic effects [164], and size exclusion [165] properties provide opportunity for their application in optimizing two coupled tasks of desalination mentioned above. Specifically, high selective-solute retention and high permeability to water is achieved by drilling sub-nanometer holes through the solid-state nanomaterials such as graphene by using focus ion beam [166, 167] or oxygen plasma etching process [168]. The water molecules are too big to pass through graphene’s fine mesh without making holes/pores in the graphene sheet and poking these sub-nanometer holes uniformly is another challenge. Additionally, in order to drive satisfactory magnitude of water flux through these nanopores, a large pressure difference is required that cannot be solely achieved by osmotic pressure. However, compared to the fine mesh of nanoporous graphene, Nair et al. [162] proposed graphene oxide membrane as an alternate to the graphene because graphene oxide relies on its naturally occurring tortuous path to flow out water (Figure 9), unlike the challenging task of artificially poking holes through graphene for creating a path to flow out water. Other than this major advantage of graphene oxide, there are several other reasons as presented by You et al. [169] that make graphene oxide a preferred choice over nanoporous graphene.(1) Key Challenge. One of the key design challenges in the desalination technology using RO membranes is the binding of the following two conditions: (i) increasing retention for selective solutes through sieving using small nanometer-sized apertures, and (ii) allowing large flux of pure water to pass through the membrane by increasing membrane permeability. These two challenges can be visualized through a schematic [170] in Figure 10 that shows the nanosized pores of the graphene membrane and water molecules with salt are comparable to each other.Figure 10
Schematic of RO membrane with nanopores [170].Finally, a recent study by Werber et al. [171] questioned and critiqued the coupled-importance of a high selective-solute retention, and high water permeability in the desalination performance. They suggest that increased water permeability seems to have little impact on the efficiency and that increased retention of selective solutes should be the most important measure for desalination efficiency. The water/salt selectivity tradeoff relationship is based on empirical relationship [172] that is still continued to be used [173] because of incomplete theoretical understanding of many physical/chemical processes that occur during desalination through membrane-based RO.(2) Governing Equations. Most of the studies use the simple convective flux and diffusive flux to model the transport of water and solute, respectively, through the nanoscale membranes as reviewed by Werber et al. [174]:(48)Jw=AΔP−Δπm,Js=BΔcm,where A=PwVw/δmRgT and B=Ps/δm. Here, AandB are the permeability coefficients for water and solute, respectively, PwandPs are the diffusive permeability for water and solute, respectively, Vw is the molar volume of water, Rg is the gas constant, T is the absolute temperature, ΔP is the applied hydraulic pressure difference across the membrane, Δπm is the osmotic pressure difference across the active layer, and δm is the active layer thickness of the membrane. Werber et al. [174] suggest that the coefficients A and B together define the selective layer performance of nonporous membranes.


## 4.2.3. Critical Gaps and Proposed Approach

Despite the rapid progress made in developing nanomaterials for membranes, very little research has been performed in modeling the flow of water and solute retention through graphene-based membranes. Based on above discussions of liquid flow at nanoscale, water and solute transport in graphene-based membranes should be a strong function of wettability and slip length. The strong effect of membrane wettability and slip length on water permeability and volumetric flux was shown recently by Shahbabaei et al. [175] and Chen et al. [176], respectively, through molecular dynamics simulations. Werber et al. [174] reviewed important polymeric membrane materials as a measure of controlling the wettability from an experimental perspective. Therefore, in this regard it is important to understand how wettability may affect the slip length and a method to model the effect of their coupled behavior on flux.(1) Coupled Effect of Slip Length and Wettability. Bakli and Chakraborty (2012) presented a constitutive model to describe the liquid filling of water in nanopores as a function of the dynamic slip length that changes with the contact angle and the fluid acceleration due to pressure gradient as follows:(49)Ls=A1a1+cosθ2−A2a,a=γcosθρRL,a∗=γcosθηLdL/dt,where Ls is the slip length, L is the water filling length at time t, N is the parameter representing the surface roughness (lower value represents higher roughness), a is the imbibition driving acceleration parameter, a∗ is the dimensionless form of a, and A1a and A2a are functions of a.Parobek and Liu [177] reviewed the wettability of graphene and found that the intrinsic wettability of graphene based on advancing contact angle measurements suggests that graphene is a nonwetting (hydrophobic) material and its wettability is independent of the substrate underneath graphene and the number of graphene layers. Based on our earlier discussions, it is known that slip lengths on the order of tens of nanometers are observed on hydrophobic (nonwetting) surfaces that can go up to hundreds of nanometers [178, 179] for super-hydrophobic surfaces. Therefore, the hydrophobic nature of graphene further emphasizes the need to incorporate the slip length in modeling the flux of water, which can be obtained using the Hagen-Poiseuille (HP) equation for cylindrical-shaped capillary or using the HP equation modified for geometric corrections and slip flow through capillaries of different shapes:(50)q=cπr4+4r3Ls8μ.1L0,where c is a geometry correction factor [180] to account for different shapes and eccentricity of capillaries (k=1,1.43,and1.98 for capillaries with circular, square, and equilateral triangle cross-sections, respectively), r is the pore radius, μ is the water viscosity, Ls is the slip length, and L0 is the capillary length within the membrane.


## 4.3. Shale Reservoirs

### 4.3.1. Introduction

Shales are unconventional rocks that are both the source and the reservoir for oil and gas resources. Economical production from shale formations requires creating hydraulic fractures by pumping large volume of water (∼2–6 million gallons) at high pressure. Shale rocks are predominantly composed of consolidated clay-sized particles with a high organic content [181, 182] that result in extremely small pore sizes with 80% of pores in the range of 2 nm to 15 nm [183, 184]. Because of their extremely small pore sizes, shale rocks present an opportunity to apply the fundamentals of nanoscale fluid flow physics over existing conventional theories of porous media fluid flow.


### 4.3.2. Literature Review

The typical characteristics of a shale reservoir include low porosity, low permeability [185, 186], complex network of matrix-fracture systems [187–191], and heterogeneous mineralogy [181, 182, 192–195], and the combination of these characteristics presents a unique challenge in understanding the fluid flow and fluid-rock and fluid-fluid interactions in these reservoirs. Following theoretical considerations, some studies explored the fluid flow in these reservoirs as a function of convection, slip flow, Knudsen diffusion [196–201], surface diffusion/sorption [202–210], rock-mechanics [189, 211–215], and osmosis [216–218]. Unlike conventional reservoir rocks, shales are typically oil-wet or mixed-wet as discussed in a review by Singh [218] for various shale rocks across the world [219–222]. The wettability of shales cannot be characterized by the macroscopic method of measuring contact angle because of the nontrivial effect of line tension in nanopores, and to account for this effect, a modified contact angle method was discussed in a review by Singh [218].Other areas that have recently caught interest due to their importance in understanding the fluid behavior in shale or increasing the productivity are (i) spontaneous imbibition in shales [218, 223–225], (ii) thermodynamic phase behavior of fluids inside shales as a function of pressure drawdown during production [226–231], (iii) alternate fracturing fluids [232–237] for better reservoir stimulation and to avoid water usage, and (iv) enhancing oil recovery by injection of fluids [238–240].(1) Key Challenge. The knowledge about fluid flow in shale reservoirs is still at its infancy compared to our knowledge of fluid flow in conventional oil and gas reservoirs that has been rigorously tested over several decades. The key challenge in production from shale reservoirs is developing verifiable theories that describe flow and mechanics from these extremely tigh pore spaces, which are heterogeneous in terms of the mineralogical composition and pore structure.


### 4.3.3. Critical Gaps and Proposed Approach

Even though the importance of shale reservoirs as a source of hydrocarbon production continues to increase with time, the production of fluids from these reservoirs remain poorly understood at the pore level and that drives many researchers to look at these problems from continuum perspective at the reservoir scale [241]. From a theoretical perspective, there are many transport mechanisms that seem to be contributing to the hydrocarbon production from the shale rocks. A better understanding to characterize those mechanisms in terms of their contribution to the flow, the period for which they remain activated, and their likely place of origin in the reservoir can be useful to optimize the production from these reservoirs. In this regard, we can elaborate this into three specific questions: (i) What is the contribution of each transport mechanism in the total production? (ii) At what point of time and for how long during a reservoir’s life do these mechanisms remain activated? (iii) What is the contribution of different sections of reservoir, such as matrix, natural fractures, and hydraulic fractures, in the transport of hydrocarbons?(1) Reservoir-on-Chip Approach. It is clear from the discussions in the literature that ultra-tight nanopore structure of the shales plays a major role in controlling the transport of fluids and their production, but despite this consensus most of the experimental studies focus on investigating the fluid transport from a continuum perspective using core-scale experiments [219, 242–246]. The core-scale investigations can be employed to explore the flow due to interaction between matrix and fractures; however, the appropriate way to investigate nanoscale transport mechanisms would be through the nanofluidics (e.g., chip-based models) framework. Recently, there has been some effort [159, 226, 247–252] in investigating the transport from nanopores using some innovative ideas that mimic the heterogeneous nanopore and fracture network through the reservoir-on-chip approach as illustrated in Figure 11. A more realistic network of matrix and fracture mimicking actual core images can be developed synthetically using micro-CT images, characterizing the fracture and matrix system in those images using an image processing software and fabricating the processed images using standard lithography techniques [247, 253]. The features of fracture network from the CT images can be patterned and etched using deep reactive ion etching. An innovative and sophisticated reservoir-on-chip model was developed by Gerami et al. [247] to study coal reservoirs, and a similar approach can be adopted to develop a reservoir-on-chip model for shale rocks.Figure 11
Schematic illustrating (a) a reservoir-on-chip framework and (b) a micromodel with a network of fractures and heterogeneous pores (inset) (adapted after [252] and [248], respectively).(2) Contribution of Liquid Slip to Total Flow. Discovery of slip flow’s contribution in enhancing performance in other micro- and nanoscale applications, for instance, energy conversion devices [143, 148], has been acknowledged with great enthusiasm by others [153] in their field. However, similar interest has not been visible in the research community studying flow from shale reservoirs, despite some theoretical studies that illustrated the importance of gas slip in shale [196, 206, 254]. Until now, there have been no experimental studies to investigate the theoretical contribution of gas slip to the total flux. Theoretical studies considering slip flow of oil and its contribution to total flux are still evading. Perhaps, the major reason why micro- and nanoscale transport mechanisms such as slip flow have not been experimentally investigated by shale research community is because of the general inclination to perform experiments using continuum scale cores that are not suitable to investigate the fluid slip. Therefore, it is suggested that to study the micro- and nanoscale transport mechanisms in shale rocks, an emphasis must be given to employ the reservoir-on-chip approach discussed earlier. Following that, a systematic approach can be adopted to study the effect of gas and liquid slip on the total flow. As a starting point, well-established theories on liquid slip can be used to verify their applicability in shales; for instance, it can be investigated whether slip length (Ls) for oil varies as Ls∝1+cosθ−2 in synthetically fabricated chips mimicking shale rock structure by varying its wettability (θ). If the surface of the fabricated chip has been characterized with some roughness function fx,y, then a modified expression for slip length accounting for roughness [128] can be used as follows:(51)Ls=k1+cosθ21+α∗2,α∗=∑fx,y−fx,y¯N,where k is the constant of proportionality and α∗ is a parameter that characterizes the roughness of a surface through its amplitude N and roughness function fx,y. Now, the flow rate as a function of slip length would be given as follows:(52)Qslip=Qno−slip1+6Ls2h,where h is the distance between the confining pore walls and Qslip and Qno−slip are the flow rates due to slip and no-slip boundary conditions, respectively.(3) Contribution of Evaporation to Water Loss. Hydraulic fracturing requires injecting large volume of water (∼2–6 million gallons); however, on average, only 6–10% of the injected water is recovered in the US across all shale plays [255]. This abundant retention of fracturing fluid inside the shale formations is a cause of major concern because it keeps the hydrocarbons from flowing out of the reservoir by reducing the relative permeability of the hydrocarbons inside the formation. Singh [218] reviewed the process of water uptake by shales and suggested evaporation as one of the possible mechanisms that may be contributing to loss of water in shales. To support their hypothesis, they provided three characteristics associated with shale formations that could facilitate evaporation of water: (1) gas expansion, (2) water droplets formed at asperities or cavities in fractures, and (3) pore confinement. Since the hypothesis of water loss in shales due to evaporation has not been investigated before, it can be tested using the reservoir-on-chip approach with the network of small pores and fractures. Among the three characteristics noted above that could be facilitating evaporation, gas expansion seems to be the most intuitive source and can be mimicked by injecting gas into the water saturated chip with a network of matrix pores and fractures under reservoir conditions. The process of fabricating such a chip and other experimental devices needed to test this process can be set up based on the reservoir-on-chip model developed by Gerami et al. [247] to study coal reservoirs and shown in Figure 12.Figure 12
Reservoir-on-chip model used to depict the fracture structure in a coal reservoir [247].Keeping track of water saturation and its phase change with time can indicate whether expansion of gas as it moves from pores to fractures results in evaporating some water. The phase change can be studied using direct visual observation of water phase behavior due to light scattering, a method that has been reliably tested for vapor detection by Yang et al. [256] and more recently in the context of shales by Ally et al. [227] to study the fluid phase behavior in nanopores.


## 4.3.1. Introduction

Shales are unconventional rocks that are both the source and the reservoir for oil and gas resources. Economical production from shale formations requires creating hydraulic fractures by pumping large volume of water (∼2–6 million gallons) at high pressure. Shale rocks are predominantly composed of consolidated clay-sized particles with a high organic content [181, 182] that result in extremely small pore sizes with 80% of pores in the range of 2 nm to 15 nm [183, 184]. Because of their extremely small pore sizes, shale rocks present an opportunity to apply the fundamentals of nanoscale fluid flow physics over existing conventional theories of porous media fluid flow.


## 4.3.2. Literature Review

The typical characteristics of a shale reservoir include low porosity, low permeability [185, 186], complex network of matrix-fracture systems [187–191], and heterogeneous mineralogy [181, 182, 192–195], and the combination of these characteristics presents a unique challenge in understanding the fluid flow and fluid-rock and fluid-fluid interactions in these reservoirs. Following theoretical considerations, some studies explored the fluid flow in these reservoirs as a function of convection, slip flow, Knudsen diffusion [196–201], surface diffusion/sorption [202–210], rock-mechanics [189, 211–215], and osmosis [216–218]. Unlike conventional reservoir rocks, shales are typically oil-wet or mixed-wet as discussed in a review by Singh [218] for various shale rocks across the world [219–222]. The wettability of shales cannot be characterized by the macroscopic method of measuring contact angle because of the nontrivial effect of line tension in nanopores, and to account for this effect, a modified contact angle method was discussed in a review by Singh [218].Other areas that have recently caught interest due to their importance in understanding the fluid behavior in shale or increasing the productivity are (i) spontaneous imbibition in shales [218, 223–225], (ii) thermodynamic phase behavior of fluids inside shales as a function of pressure drawdown during production [226–231], (iii) alternate fracturing fluids [232–237] for better reservoir stimulation and to avoid water usage, and (iv) enhancing oil recovery by injection of fluids [238–240].(1) Key Challenge. The knowledge about fluid flow in shale reservoirs is still at its infancy compared to our knowledge of fluid flow in conventional oil and gas reservoirs that has been rigorously tested over several decades. The key challenge in production from shale reservoirs is developing verifiable theories that describe flow and mechanics from these extremely tigh pore spaces, which are heterogeneous in terms of the mineralogical composition and pore structure.


## 4.3.3. Critical Gaps and Proposed Approach

Even though the importance of shale reservoirs as a source of hydrocarbon production continues to increase with time, the production of fluids from these reservoirs remain poorly understood at the pore level and that drives many researchers to look at these problems from continuum perspective at the reservoir scale [241]. From a theoretical perspective, there are many transport mechanisms that seem to be contributing to the hydrocarbon production from the shale rocks. A better understanding to characterize those mechanisms in terms of their contribution to the flow, the period for which they remain activated, and their likely place of origin in the reservoir can be useful to optimize the production from these reservoirs. In this regard, we can elaborate this into three specific questions: (i) What is the contribution of each transport mechanism in the total production? (ii) At what point of time and for how long during a reservoir’s life do these mechanisms remain activated? (iii) What is the contribution of different sections of reservoir, such as matrix, natural fractures, and hydraulic fractures, in the transport of hydrocarbons?(1) Reservoir-on-Chip Approach. It is clear from the discussions in the literature that ultra-tight nanopore structure of the shales plays a major role in controlling the transport of fluids and their production, but despite this consensus most of the experimental studies focus on investigating the fluid transport from a continuum perspective using core-scale experiments [219, 242–246]. The core-scale investigations can be employed to explore the flow due to interaction between matrix and fractures; however, the appropriate way to investigate nanoscale transport mechanisms would be through the nanofluidics (e.g., chip-based models) framework. Recently, there has been some effort [159, 226, 247–252] in investigating the transport from nanopores using some innovative ideas that mimic the heterogeneous nanopore and fracture network through the reservoir-on-chip approach as illustrated in Figure 11. A more realistic network of matrix and fracture mimicking actual core images can be developed synthetically using micro-CT images, characterizing the fracture and matrix system in those images using an image processing software and fabricating the processed images using standard lithography techniques [247, 253]. The features of fracture network from the CT images can be patterned and etched using deep reactive ion etching. An innovative and sophisticated reservoir-on-chip model was developed by Gerami et al. [247] to study coal reservoirs, and a similar approach can be adopted to develop a reservoir-on-chip model for shale rocks.Figure 11
Schematic illustrating (a) a reservoir-on-chip framework and (b) a micromodel with a network of fractures and heterogeneous pores (inset) (adapted after [252] and [248], respectively).(2) Contribution of Liquid Slip to Total Flow. Discovery of slip flow’s contribution in enhancing performance in other micro- and nanoscale applications, for instance, energy conversion devices [143, 148], has been acknowledged with great enthusiasm by others [153] in their field. However, similar interest has not been visible in the research community studying flow from shale reservoirs, despite some theoretical studies that illustrated the importance of gas slip in shale [196, 206, 254]. Until now, there have been no experimental studies to investigate the theoretical contribution of gas slip to the total flux. Theoretical studies considering slip flow of oil and its contribution to total flux are still evading. Perhaps, the major reason why micro- and nanoscale transport mechanisms such as slip flow have not been experimentally investigated by shale research community is because of the general inclination to perform experiments using continuum scale cores that are not suitable to investigate the fluid slip. Therefore, it is suggested that to study the micro- and nanoscale transport mechanisms in shale rocks, an emphasis must be given to employ the reservoir-on-chip approach discussed earlier. Following that, a systematic approach can be adopted to study the effect of gas and liquid slip on the total flow. As a starting point, well-established theories on liquid slip can be used to verify their applicability in shales; for instance, it can be investigated whether slip length (Ls) for oil varies as Ls∝1+cosθ−2 in synthetically fabricated chips mimicking shale rock structure by varying its wettability (θ). If the surface of the fabricated chip has been characterized with some roughness function fx,y, then a modified expression for slip length accounting for roughness [128] can be used as follows:(51)Ls=k1+cosθ21+α∗2,α∗=∑fx,y−fx,y¯N,where k is the constant of proportionality and α∗ is a parameter that characterizes the roughness of a surface through its amplitude N and roughness function fx,y. Now, the flow rate as a function of slip length would be given as follows:(52)Qslip=Qno−slip1+6Ls2h,where h is the distance between the confining pore walls and Qslip and Qno−slip are the flow rates due to slip and no-slip boundary conditions, respectively.(3) Contribution of Evaporation to Water Loss. Hydraulic fracturing requires injecting large volume of water (∼2–6 million gallons); however, on average, only 6–10% of the injected water is recovered in the US across all shale plays [255]. This abundant retention of fracturing fluid inside the shale formations is a cause of major concern because it keeps the hydrocarbons from flowing out of the reservoir by reducing the relative permeability of the hydrocarbons inside the formation. Singh [218] reviewed the process of water uptake by shales and suggested evaporation as one of the possible mechanisms that may be contributing to loss of water in shales. To support their hypothesis, they provided three characteristics associated with shale formations that could facilitate evaporation of water: (1) gas expansion, (2) water droplets formed at asperities or cavities in fractures, and (3) pore confinement. Since the hypothesis of water loss in shales due to evaporation has not been investigated before, it can be tested using the reservoir-on-chip approach with the network of small pores and fractures. Among the three characteristics noted above that could be facilitating evaporation, gas expansion seems to be the most intuitive source and can be mimicked by injecting gas into the water saturated chip with a network of matrix pores and fractures under reservoir conditions. The process of fabricating such a chip and other experimental devices needed to test this process can be set up based on the reservoir-on-chip model developed by Gerami et al. [247] to study coal reservoirs and shown in Figure 12.Figure 12
Reservoir-on-chip model used to depict the fracture structure in a coal reservoir [247].Keeping track of water saturation and its phase change with time can indicate whether expansion of gas as it moves from pores to fractures results in evaporating some water. The phase change can be studied using direct visual observation of water phase behavior due to light scattering, a method that has been reliably tested for vapor detection by Yang et al. [256] and more recently in the context of shales by Ally et al. [227] to study the fluid phase behavior in nanopores.


## 4.4. Hydrogen Storage

### 4.4.1. Introduction

Hydrogen-based energy is one of the promising sources of energy [257] because hydrogen is a zero-emission fuel and can provide solutions to many versatile problems such as predictable supply that is not dependent on weather conditions, its capability to integrate within existing infrastructure developed for other forms of energy, and its applications in many different industries [258, 259]. As such, hydrogen-based energy solutions have versatile applications, such as sources of power [260] (fuel cells [261]) in transportation [262, 263] and as a commodity in many industrial processes [264]. Figure 13 shows a general illustration of a hydrogen energy system with sources, production, storage, and its applications. Among all sources of hydrogen production worldwide, 95% of the hydrogen production comes from nonrenewable sources [265] such as natural gas, coal, gasoline, diesel fuels, etc.Figure 13
Schematic representing a hydrogen energy system. About 95% of the hydrogen currently produced across the world comes from nonrenewable sources.Among various steps in the hydrogen energy system, hydrogen storage is a major challenge and a key enabling technology for the advancement of the hydrogen energy-based applications.


### 4.4.2. Literature Review

Hydrogen can be stored either physically (as compressed gas or liquid) or chemically (on the surfaces of solid materials by adsorption or within the solid materials by absorption). Exploiting hydrogen energy on commercial scale as planned by the US Department of Energy through portable power sources, fuel cells, and transportation requires a compact system that can store hydrogen on these applications. The fundamental difficulty in achieving this task is low energy density per unit volume of hydrogen, which means hydrogen needs relatively larger volume to pack same amount of energy compared to other fuels. The storage of hydrogen can be divided in three sectors [266–269], namely: (i) high-pressured gaseous state, (ii) cryogenic liquid state, and (iii) solid state. Because of the low storage efficiency of liquid and gaseous hydrogen, research efforts to develop storage technologies have focused towards increasing the density of hydrogen either through chemical storage technology using new materials [270–273] or through physical storage using compression [274–276].Even though some researchers [271] suggest that the chemical method is more efficient in storing hydrogen than the physical method, they do not provide a comparison or a reference to support such a claim. However, Petitpas et al. [275] provided an objective comparison between the two storage methods, i.e., physical storage and chemical storage, which suggested that at low pressures the density of stored hydrogen through the chemical method is significantly higher than the density of hydrogen stored through the physical method. As the pressure is increased, the density of hydrogen through the chemical method starts to flatten, whereas the density of hydrogen through the physical method continues to increase, eventually overtaking the density of hydrogen stored through the chemical method at a breakeven pressure beyond which the density of hydrogen stored through the physical method remains larger than the hydrogen stored through the chemical method (Figure 14). They reported that this breakeven pressure strongly depends on temperature.Figure 14
Comparison of physical and chemical storage capacities [275].Since the choice of hydrogen storage technology would depend on the end application, Petitpas et al. [275] suggested that the conditions relevant to a specific application may supersede the chemistry and physics required to enhance the hydrogen density. For instance, the added mass of the adsorbent for chemical storage may discourage the use of this storage technology in the case of weight-sensitive applications, whereas in the case of storage by compression the use of higher pressures may discourage its use for some applications where high pressure is a safety concern. Additionally, the cryogenic temperature needed for both techniques could be discouraging in the absence of liquid hydrogen (or liquid nitrogen for 80 K fueling). Petitpas et al. [275] also compared the use of a “hybrid” system that can store hydrogen by both compression and adsorbent and found that they can store 10% to 20% (with metal-organic framework as the sorbent) more hydrogen compared to storage by compression alone.Despite the variation in density of stored hydrogen by the two approaches, the research on hydrogen storage in past few years has mostly focused on chemical storage by adsorbents because of its expected use in many commercial and portable applications that may have safety concerns related to high pressure. Therefore, below we review the modeling scheme to store hydrogen by adsorption on porous and deformable materials.(1) Key Challenge. Increasing the efficiency of hydrogen storage is one of the major challenges in the application of hydrogen for consumer applications such as transportation. Figure 15 shows a schematic comparing the volumetric scale of storage space required for 4 kg of hydrogen storage in different materials and storage methods.Figure 15
Schematic representation of the volumetric scale of storage space required for 4 kg of hydrogen with different modes of storage [259, 268]. Car image is a courtesy of Toyota press information, 33rd Tokyo Motor Show.(2) Governing Equations. The governing equations for the storage of hydrogen follow the conservation of mass and energy.Mass balance:(53)∂ϕρg∂t+∇.ϕρgu=−∂sm∂t,∂sm∂t=∂1−ϕρs∂t=∂ρs∂t−∂ϕρs∂t,where ρg is the hydrogen gas density, ϕ is the porosity of the adsorbent, u is the hydrogen gas velocity that is computed using either Navier–Stokes equation or Darcy’s law, sm is the mass of hydrogen adsorbed or desorbed per unit volume of the sorbent, and ρs is the density of the sorbent.Energy balance:(54)∂ρcp¯T∂t+∇.ρgcp,guT=∇.keff∇T+∂sT∂t,where ρcp¯ is the effective heat capacity, keff is the effective thermal conductivity [277], T is the temperature, and sT is the product of the enthalpy change (as a result of sorption).


### 4.4.3. Critical Gaps and Proposed Approach

Even though the experimental research on developing superior nanostructured adsorbents [271, 273, 278, 279] for hydrogen storage continues to increase with time, the same effort has not been devoted to model the transport of hydrogen in these nanostructured porous materials. For instance, transport of hydrogen is modeled using only the continuum scale transport [277, 280–283] with no contribution of nanoscale transport mechanisms despite the role of nanoconfinement [273, 284] as an effective strategy in modifying the hydrogen storage performance. From a theoretical perspective, transport of hydrogen in nanoporous mediums is controlled by two more mechanisms (Knudsen diffusion and slip flow) other than viscous/convective flow, heat transfer, and sorption (surface diffusion). A better understanding to characterize the transport mechanisms at nanoscale in terms of their contribution to the storage and the operating conditions for which they remain activated can be useful to optimize the storage of hydrogen in nanostructured materials. In this regard, we discuss two transport mechanisms of Knudsen diffusion and slip flow at nanoscale and how they can be added in the governing equations to model hydrogen storage.(1) Accounting for Slip Flow and Knudsen Diffusion in Total Flux. The flow behavior at the continuum regime is governed by the Navier–Stokes equations; however, as Kn increases beyond the range of the continuum regime, the no-slip boundary conditions start to break such that the N-S equations should be modified to account for slip boundary conditions or molecular models should be employed. Another transport mechanism that plays an important role when the mean free path of the gas molecules is larger than or equal to the pore size is the form of diffusive transport called Knudsen diffusion. Diffusive transport in micro- and nanopores can occur through three different mechanisms: (i) bulk diffusion, (ii) surface diffusion (sorption), and (iii) Knudsen diffusion [67]. The proportion of each mechanism depends on the magnitude of the Knudsen number Kn, and for large Kn, the molecules interact more often with the pore walls than with each other. Typical nanoporous materials for storing hydrogen have pore size that can vary from less than 1 nm to about 10 nm [270, 271, 273, 285, 286] for varied materials ranging from carbons, zeolites, metal-organic framework (MOF), complex hydrides, etc., whereas mean free path of hydrogen molecules at ambient conditions estimated from kinetic theory is about 160 nm. Therefore, it is evident from this comparison that Knudsen diffusion may have an important role to play in the transport of hydrogen.(i) Slip Flow and Knudsen Diffusion: The convective mass flow in the presence of slip boundary conditions and diffusive mass flow in the presence of Knudsen diffusion is given as follows [69]:(55)m˙C=MuARTp0l,m˙D=DKnMARTΔpl=13s8RTπMMARTΔpl,s=4hπ,channel  shaped  pores,2r,tube  shaped  pores.The solution of the Navier–Stokes equation and slip flow boundary conditions are used to solve for hydrogen velocity due to convective flux for two pore shapes of channels and tubes, respectively, as follows:(56)u=2h∫−h/2h/212μ∂p∂xy2+c1y+c2dy,channel‐shaped  pores,1πr∗2∫0r∗14μ∂p∂xr2+c3lnr+c42πrdr,tube‐shaped  pores,uslip=2−σuσuλ∂u∂yy=h,channel‐shaped  pores,2−σuσuλ∂u∂rr=r∗,tube‐shaped  pores,where m˙D is the diffusive mass flow, DKn is the Knudsen diffusion coefficient, s is the characteristic length, p is the pressure, u is the hydrogen velocity due to convective flux, σu is TMAC, λ is the mean free path, p0 is the outlet pressure, M is the molar mass, A is the cross-sectional area of the medium, h is the height of the channel-shaped pore, r∗ is the radius of a cylindrical-shaped pore, l is the length of the medium, R is the universal gas constant, T is the temperature, and c1,c2,c3,andc4 are constants from integration of the Navier–Stokes equation.The total flux is the sum of the mass flow due to convective flux (including slip flow) and diffusive flux (from Knudsen diffusion), eventually leading to an effective velocity of hydrogen molecules (ueff) that can be substituted in place of u in the mass balance equation for hydrogen:(57)m˙=m˙C+m˙D,⇒m˙=MuARTp0l+DKnMARTΔpl,ueff=m˙ρgA,ueff=MρgRTlup0+DKnΔp.


## 4.4.1. Introduction

Hydrogen-based energy is one of the promising sources of energy [257] because hydrogen is a zero-emission fuel and can provide solutions to many versatile problems such as predictable supply that is not dependent on weather conditions, its capability to integrate within existing infrastructure developed for other forms of energy, and its applications in many different industries [258, 259]. As such, hydrogen-based energy solutions have versatile applications, such as sources of power [260] (fuel cells [261]) in transportation [262, 263] and as a commodity in many industrial processes [264]. Figure 13 shows a general illustration of a hydrogen energy system with sources, production, storage, and its applications. Among all sources of hydrogen production worldwide, 95% of the hydrogen production comes from nonrenewable sources [265] such as natural gas, coal, gasoline, diesel fuels, etc.Figure 13
Schematic representing a hydrogen energy system. About 95% of the hydrogen currently produced across the world comes from nonrenewable sources.Among various steps in the hydrogen energy system, hydrogen storage is a major challenge and a key enabling technology for the advancement of the hydrogen energy-based applications.


## 4.4.2. Literature Review

Hydrogen can be stored either physically (as compressed gas or liquid) or chemically (on the surfaces of solid materials by adsorption or within the solid materials by absorption). Exploiting hydrogen energy on commercial scale as planned by the US Department of Energy through portable power sources, fuel cells, and transportation requires a compact system that can store hydrogen on these applications. The fundamental difficulty in achieving this task is low energy density per unit volume of hydrogen, which means hydrogen needs relatively larger volume to pack same amount of energy compared to other fuels. The storage of hydrogen can be divided in three sectors [266–269], namely: (i) high-pressured gaseous state, (ii) cryogenic liquid state, and (iii) solid state. Because of the low storage efficiency of liquid and gaseous hydrogen, research efforts to develop storage technologies have focused towards increasing the density of hydrogen either through chemical storage technology using new materials [270–273] or through physical storage using compression [274–276].Even though some researchers [271] suggest that the chemical method is more efficient in storing hydrogen than the physical method, they do not provide a comparison or a reference to support such a claim. However, Petitpas et al. [275] provided an objective comparison between the two storage methods, i.e., physical storage and chemical storage, which suggested that at low pressures the density of stored hydrogen through the chemical method is significantly higher than the density of hydrogen stored through the physical method. As the pressure is increased, the density of hydrogen through the chemical method starts to flatten, whereas the density of hydrogen through the physical method continues to increase, eventually overtaking the density of hydrogen stored through the chemical method at a breakeven pressure beyond which the density of hydrogen stored through the physical method remains larger than the hydrogen stored through the chemical method (Figure 14). They reported that this breakeven pressure strongly depends on temperature.Figure 14
Comparison of physical and chemical storage capacities [275].Since the choice of hydrogen storage technology would depend on the end application, Petitpas et al. [275] suggested that the conditions relevant to a specific application may supersede the chemistry and physics required to enhance the hydrogen density. For instance, the added mass of the adsorbent for chemical storage may discourage the use of this storage technology in the case of weight-sensitive applications, whereas in the case of storage by compression the use of higher pressures may discourage its use for some applications where high pressure is a safety concern. Additionally, the cryogenic temperature needed for both techniques could be discouraging in the absence of liquid hydrogen (or liquid nitrogen for 80 K fueling). Petitpas et al. [275] also compared the use of a “hybrid” system that can store hydrogen by both compression and adsorbent and found that they can store 10% to 20% (with metal-organic framework as the sorbent) more hydrogen compared to storage by compression alone.Despite the variation in density of stored hydrogen by the two approaches, the research on hydrogen storage in past few years has mostly focused on chemical storage by adsorbents because of its expected use in many commercial and portable applications that may have safety concerns related to high pressure. Therefore, below we review the modeling scheme to store hydrogen by adsorption on porous and deformable materials.(1) Key Challenge. Increasing the efficiency of hydrogen storage is one of the major challenges in the application of hydrogen for consumer applications such as transportation. Figure 15 shows a schematic comparing the volumetric scale of storage space required for 4 kg of hydrogen storage in different materials and storage methods.Figure 15
Schematic representation of the volumetric scale of storage space required for 4 kg of hydrogen with different modes of storage [259, 268]. Car image is a courtesy of Toyota press information, 33rd Tokyo Motor Show.(2) Governing Equations. The governing equations for the storage of hydrogen follow the conservation of mass and energy.Mass balance:(53)∂ϕρg∂t+∇.ϕρgu=−∂sm∂t,∂sm∂t=∂1−ϕρs∂t=∂ρs∂t−∂ϕρs∂t,where ρg is the hydrogen gas density, ϕ is the porosity of the adsorbent, u is the hydrogen gas velocity that is computed using either Navier–Stokes equation or Darcy’s law, sm is the mass of hydrogen adsorbed or desorbed per unit volume of the sorbent, and ρs is the density of the sorbent.Energy balance:(54)∂ρcp¯T∂t+∇.ρgcp,guT=∇.keff∇T+∂sT∂t,where ρcp¯ is the effective heat capacity, keff is the effective thermal conductivity [277], T is the temperature, and sT is the product of the enthalpy change (as a result of sorption).


## 4.4.3. Critical Gaps and Proposed Approach

Even though the experimental research on developing superior nanostructured adsorbents [271, 273, 278, 279] for hydrogen storage continues to increase with time, the same effort has not been devoted to model the transport of hydrogen in these nanostructured porous materials. For instance, transport of hydrogen is modeled using only the continuum scale transport [277, 280–283] with no contribution of nanoscale transport mechanisms despite the role of nanoconfinement [273, 284] as an effective strategy in modifying the hydrogen storage performance. From a theoretical perspective, transport of hydrogen in nanoporous mediums is controlled by two more mechanisms (Knudsen diffusion and slip flow) other than viscous/convective flow, heat transfer, and sorption (surface diffusion). A better understanding to characterize the transport mechanisms at nanoscale in terms of their contribution to the storage and the operating conditions for which they remain activated can be useful to optimize the storage of hydrogen in nanostructured materials. In this regard, we discuss two transport mechanisms of Knudsen diffusion and slip flow at nanoscale and how they can be added in the governing equations to model hydrogen storage.(1) Accounting for Slip Flow and Knudsen Diffusion in Total Flux. The flow behavior at the continuum regime is governed by the Navier–Stokes equations; however, as Kn increases beyond the range of the continuum regime, the no-slip boundary conditions start to break such that the N-S equations should be modified to account for slip boundary conditions or molecular models should be employed. Another transport mechanism that plays an important role when the mean free path of the gas molecules is larger than or equal to the pore size is the form of diffusive transport called Knudsen diffusion. Diffusive transport in micro- and nanopores can occur through three different mechanisms: (i) bulk diffusion, (ii) surface diffusion (sorption), and (iii) Knudsen diffusion [67]. The proportion of each mechanism depends on the magnitude of the Knudsen number Kn, and for large Kn, the molecules interact more often with the pore walls than with each other. Typical nanoporous materials for storing hydrogen have pore size that can vary from less than 1 nm to about 10 nm [270, 271, 273, 285, 286] for varied materials ranging from carbons, zeolites, metal-organic framework (MOF), complex hydrides, etc., whereas mean free path of hydrogen molecules at ambient conditions estimated from kinetic theory is about 160 nm. Therefore, it is evident from this comparison that Knudsen diffusion may have an important role to play in the transport of hydrogen.(i) Slip Flow and Knudsen Diffusion: The convective mass flow in the presence of slip boundary conditions and diffusive mass flow in the presence of Knudsen diffusion is given as follows [69]:(55)m˙C=MuARTp0l,m˙D=DKnMARTΔpl=13s8RTπMMARTΔpl,s=4hπ,channel  shaped  pores,2r,tube  shaped  pores.The solution of the Navier–Stokes equation and slip flow boundary conditions are used to solve for hydrogen velocity due to convective flux for two pore shapes of channels and tubes, respectively, as follows:(56)u=2h∫−h/2h/212μ∂p∂xy2+c1y+c2dy,channel‐shaped  pores,1πr∗2∫0r∗14μ∂p∂xr2+c3lnr+c42πrdr,tube‐shaped  pores,uslip=2−σuσuλ∂u∂yy=h,channel‐shaped  pores,2−σuσuλ∂u∂rr=r∗,tube‐shaped  pores,where m˙D is the diffusive mass flow, DKn is the Knudsen diffusion coefficient, s is the characteristic length, p is the pressure, u is the hydrogen velocity due to convective flux, σu is TMAC, λ is the mean free path, p0 is the outlet pressure, M is the molar mass, A is the cross-sectional area of the medium, h is the height of the channel-shaped pore, r∗ is the radius of a cylindrical-shaped pore, l is the length of the medium, R is the universal gas constant, T is the temperature, and c1,c2,c3,andc4 are constants from integration of the Navier–Stokes equation.The total flux is the sum of the mass flow due to convective flux (including slip flow) and diffusive flux (from Knudsen diffusion), eventually leading to an effective velocity of hydrogen molecules (ueff) that can be substituted in place of u in the mass balance equation for hydrogen:(57)m˙=m˙C+m˙D,⇒m˙=MuARTp0l+DKnMARTΔpl,ueff=m˙ρgA,ueff=MρgRTlup0+DKnΔp.


## 5. Conclusions

While there is a consensus in the literature that embracing nanodevices and nanomaterials helps in improving the efficiency and performance, the reason for the better performance is mostly subscribed to the nanosized material/structure of the system without sufficiently acknowledging the role of fluid flow mechanisms in these systems. This is evident from the literature review of fluid flow modeling in various energy-related applications, which reveals that the fundamental understanding of fluid transport at micro- and nanoscale is not adequately adapted in models. Incomplete or insufficient physics for the fluid flow can lead to untapped potential of these applications that can be used to increase their performance. This paper reviewed the current state of research for the physics of gas and liquid flow at micro- and nanoscale and identified critical gaps to improve fluid flow modeling in four different applications related to the energy sector. The review for gas flow focused on fundamentals of gas flow at rarefied conditions, the velocity slip, and temperature jump conditions. The review for liquid flow provided fundamental flow regimes of liquid flow, and liquid slip models as a function of various parameters at micro- and nanoscale. Regarding porous media applications in the energy sector, this review identified critical gaps to improve fluid flow modeling in each of the four different applications. The four applications from the energy sector considered in this review are (i) electrokinetic energy conversion devices, (ii) membrane-based water desalination through reverse osmosis, (iii) shale reservoirs, and (iv) hydrogen storage, respectively. For electrokinetic energy conversion, it was found that further improvements can be made in terms of modeling slip length as a function of key device parameters and its relationship with flux. For membrane-based water desalination, it was found that further improvements can be made if water and solute transport in graphene-based membranes are modeled by coupling the dependency of wettability and slip length in addition to typical flux from convective and diffusive flows. For shale reservoirs, it was found that researchers explore the issue at continuum scale than at pore scale, which can generally mask the phenomena contributing to flow from micro- to nanoscale; in that regard, we proposed using a reservoir-on-chip approach that can enable capturing the subcontinuum effects contributing to fluid flow in shale reservoirs. Additionally, the reservoir-on-chip approach can be used to estimate the contribution of liquid slip to total flow and the role of evaporation in the loss of fracturing water in shales. For hydrogen storage, it was found that further improvements can be made by including Knudsen diffusion and slip in the governing equations of hydrogen storage.While many experiments have been conducted to understand the contribution of some mechanisms (e.g., slip length) at micro- and nanoscale, there are just few studies that have looked at the combined effect of material properties and the size of the flowing medium on fluid flow. The coupled effect of material properties and the flow at micro- to nanoscale is important in shale rocks and electrokinetic devices where the amount of fluid flux, important for performance, is controlled by the properties of the material in addition to the scale of the flowing medium. Future work on coupled impact of material properties and flowing medium size would be an important contribution for these two applications.



---

*Source: 9565240-2018-12-04.xml*