1. Introduction
Reverse osmosis (RO) shows its high efficiency (Greenlee et al. Reference Greenlee, Lawler, Freeman, Marrot and Moulin2009; Park et al. Reference Park, Kamcev, Robeson, Elimelech and Freeman2017; Zuo et al. Reference Zuo, Wang, Deshmukh, Jia, Guo, Xin, Elimelech, Ajayan, Lou and Li2020; Di Vincenzo et al. Reference Di Vincenzo, Tiraferri, Musteata, Chisca, Deleanu, Ricceri, Cot, Nunes and Barboiu2021; DuChanois et al. Reference DuChanois, Heiranian, Yang, Porter, Li, Zhang, Verduzco and Elimelech2022; Yao et al. Reference Yao2024) and promising features in desalination (Fornasiero et al. Reference Fornasiero, Park, Holt, Stadermann, Grigoropoulos, Noy and Bakajin2008) and energy applications (Siria et al. Reference Siria, Poncharal, Biance, Fulcrand, Blase, Purcell and Bocquet2013; Janssen et al. Reference Janssen, Härtel and van Roij2014; Macha et al. Reference Macha, Marion, Nandigana and Radenovic2019; Zhang, Wen & Jiang Reference Zhang, Wen and Jiang2021; Baldock & Huang Reference Baldock and Huang2024). Ions are rejected by the steric repulsion, dehydration or dielectric energy barrier at the membrane, enabling water to permeate, thus removing hydrated ions from electrolyte solution via pressure-driven flows across the membrane. The solution-diffusion (SD) model is the most widely used theory to describe the water and ion permeation through the membrane (Paul Reference Paul1973; Wijmans & Baker Reference Wijmans and Baker1995). However, recent studies showed ions are permeable through membranes when the pressure gradient exceeds a threshold (Zhang et al. Reference Zhang, Wei, Xu and Wang2020; Wang et al. Reference Wang, Cao, Dykstra, Porada, Biesheuvel and Elimelech2021, Reference Wang, He, Heiranian, Fan, Song, Li and Elimelech2023; Biesheuvel et al. Reference Biesheuvel, Rutten, Ryzhkov, Porada and Elimelech2023), which is also influenced by salt concentration (Straub et al. Reference Straub, Osuji, Cath and Elimelech2015; Biesheuvel et al. Reference Biesheuvel, Dykstra, Porada and Elimelech2022), thereby compromising desalination selectivity. Besides, water flux has a nonlinear response to the applied pressure (Zhou & Song Reference Zhou and Song2005; Xu et al. Reference Xu, Zhu, Wang, Li, Huang, Shen, Francisco, Zeng and Meng2019; Yong et al. Reference Yong, He and Merlitz2022). These pressure-dependent permeations deviate from the constant that was predicted by the SD model.
The core assumptions of the SD model were challenged with evidence of Ångström-scale pores in membrane (Fujioka et al. Reference Fujioka, Oshima, Suzuki, Price and Nghiem2015; Culp et al. Reference Culp2021; Foglia et al. Reference Foglia, Frick, Nania, Livingston and Cabral2022), which have limitations in describing permeation behaviours (Sahu et al. Reference Sahu, Di Ventra and Zwolak2017; Qasim et al. Reference Qasim, Badrelzaman, Darwish, Darwish and Hilal2019; Wang et al. Reference Wang, Cao, Dykstra, Porada, Biesheuvel and Elimelech2021; Lu et al. Reference Lu2023; Nickerson et al. Reference Nickerson, Antonio, McNally, Toney, Ban and Straub2023; Fan et al. Reference Fan, Heiranian and Elimelech2024). Hence, new theoretical frameworks, including the pore-flow model (Sarbolouki & Miller Reference Sarbolouki and Miller1973; Takeuchi et al. Reference Takeuchi, Tanaka, Nghiem and Fujioka2018; Hegde, Doherty & Squires Reference Hegde, Doherty and Squires2022; Geise Reference Geise2022) based on an extended Nernst–Planck equation and the solution-friction model (Wang & Lin Reference Wang and Lin2021; Wang et al. Reference Wang, Cao, Dykstra, Porada, Biesheuvel and Elimelech2021, Reference Wang, He, Heiranian, Fan, Song, Li and Elimelech2023) based on hydrodynamic friction, were recently developed to explain these deviations. Additionally, coupling between water and ion flux must be considered, as demonstrated by a unified non-equilibrium simulation using a flux matrix under generic forces (Monet, Bocquet & Bocquet Reference Monet, Bocquet and Bocquet2023).
However, recent studies found ‘bizarre’ phenomena when pores/channels come to Ångström-scale regime (Bocquet Reference Bocquet2020; Kavokine, Netz & Bocquet Reference Kavokine, Netz and Bocquet2021; You et al. Reference You, Ismail, Nam, Goutham, Keerthi and Radha2022; Aluru et al. Reference Aluru2023), such as fast transport (Radha et al. Reference Radha2016; Tunuguntla et al. Reference Tunuguntla, Henley, Yao, Pham, Wanunu and Noy2017), ionic Coulomb blockade (Feng et al. Reference Feng, Liu, Graf, Dumcenco, Kis, Di Ventra and Radenovic2016; Kavokine et al. Reference Kavokine, Marbach, Siria and Bocquet2019; Xie et al. Reference Xie, Shi, Wang and Wang2023; Zhang, Fu & Xie Reference Zhang, Fu and Xie2024), proton transport (Comtet et al. Reference Comtet, Rayabharam, Glushkov, Zhang, Avsar, Watanabe, Taniguchi, Aluru and Radenovic2021) and iontronics (Robin et al. Reference Robin, Kavokine and Bocquet2021, Reference Robin2023; Shi et al. Reference Shi, Wang, Liang, Duan, Du and Xie2023), where continuum theories breakdown (Faucher et al. Reference Faucher2019). The molecular-level transport calls for a statistical mechanics framework incorporating molecular thermodynamics for modelling of water and ion permeation (Nickerson et al. Reference Nickerson, Antonio, McNally, Toney, Ban and Straub2023).
In this work, we used non-equilibrium molecular dynamics (NEMD) to investigate ion and water permeation through a single Ångströmchannel under pressure-driven flow, employing a generic Lennard-Jones (LJ) wall nanotube as a toy model. We found that water and ion permeability are pressure-dependent, with ionic fluxes rapidly increasing beyond a threshold pressure, qualitatively consistent with trends in experiments (Zhang et al. Reference Zhang, Wei, Xu and Wang2020; Biesheuvel et al. Reference Biesheuvel, Rutten, Ryzhkov, Porada and Elimelech2023; Wang et al. Reference Wang, Cao, Dykstra, Porada, Biesheuvel and Elimelech2021, Reference Wang, He, Heiranian, Fan, Song, Li and Elimelech2023) but deviating from the constant predicted by the SD model. We identified a negative differential flow resistance (NDFR) phenomenon, where flow velocity decreases with increasing pressure. This NDFR arises from the reduced escape rate of trapped ions at the entrance by joint actions of the energy barrier and hydrodynamic force that provide additional friction in the system. Trapped ions significantly reduce water permeation, emerging as a predominant factor limiting water permeation, even in low-friction systems. Since continuum theories may fail to capture the NDFR behaviour due to neglecting thermodynamic processes, we develop a non-equilibrium statistical theory to describe the ion transport and water flow through an Ångström channel for RO, capturing pressure-dependent permeance and the NDFR phenomenon using phase diagrams as functions of pressure gradient and channel size. Although this NDFR effect may statistically average out in parallel channels with a certain size distribution, the trapped ions still strongly reduce the membrane performance. Our results indicate that even under such ideal conditions, such as neutral and low frictional walls, in Ångström-scale confinement, we observe a strong suppression of water permeability due to the entrance ion trap, suggesting the ion trapping is the intrinsic physical origin for limiting the selectivity-permeability trade-off in RO. Our work offers a molecular-level understanding and statistical framework for RO that identifies ion trapping at membrane interfaces as the fundamental performance-limiting mechanism in RO. This phenomenon can be extended to much broader fields, for instance, electro-osmosis and thermo-osmosis in Ångström-scale systems when ions penetrate into the channel under generic forces.
(a) Simulation system. Pressures on pistons drive water and ion (blue for
${\rm Na}^{+}$
, orange for
${\rm Cl}^{-}$
) transport from the left reservoir through a generic LJ nanotube (cyan) to the right pure water box. Pressures are applied at the pistons, with
$P_1$
at left and
$P_2$
at right side. (b) Averaged PMF for an ion entering the nanotube, with the dashed line marking the nanotube entrance. (c) Water flow velocity over time at 2000 bar nm−1. Low flow velocity (close state) results from a trapped ion at the entrance (inset scheme) and high flow velocity (open state). (d) Flow velocity (orange squares) and ionic flux (green spheres) as functions of pressure gradient. Simulation results agree with non-equilibrium theory (solid lines). The dashed line denotes flow velocity derived from Sampson’s solution.
2. Methods
We performed NEMD simulations using Large-scale Atomic/Molecular Massively Parallel Simulator (Thompson et al. Reference Thompson2022) to study water and ionic transport through a 0.33 nm radius, 10 nm long Ångström tube clamped by two
$3.1 \times 3.1$
nm Face-Centered Cubic (FCC) membranes with 4.48 Å lattice parameter (figure 1
a). The left box contained 0.3 M NaCl, and the right contained pure water, using TIP3P water molecules and an optimised ion model for NaCl (Jorgensen et al. Reference Jorgensen, Chandrasekhar, Madura, Impey and Klein1983; Joung & Cheatham Reference Joung and Cheatham2008). We employed Lorentz–Berthelot mixing rules to calculate the heteronuclear LJ interaction parameters with solid atoms 
. The low frictional wall enables high water permeability in system. Water molecules could fill in the nanotube within 1 ns after thermal equilibrated at 289 K using a Nosé–Hoover thermostat at 1 atm by applying forces at the outermost wall. We computed long-range Coulomb interactions using the particle-particle particle-mesh (PPPM) algorithm with periodic boundary conditions in all directions. In our NEMD simulations, the pressure
$P_1$
and
$P_2$
were applied to the left and right membrane, the pressure drop
$\Delta p = P_1 - P_2$
increased from 4000 to 42 000 bar along the channel’s longitudinal (z) direction by increasing force on the leftmost membrane for hydrodynamic flows. We conducted three independent simulations for each
$\Delta p$
value, each lasting over 5 ns, to obtain water velocity and ionic fluxes. The error bars of velocity and fluxes were the standard deviation calculated from those three independent simulations. To minimise salt concentration changes in the left reservoir, data were collected only when its concentration increased by less than 10 %. In all simulations, the outermost wall was fixed in the radial plane but allowed to move in the
$z$
-direction. The rest of the wall and tube vibrated around the equilibrium position, with a harmonic spring (spring constant
) attached to each solid particle. Here we use umbrella sampling by applying biasing potentials along the axial coordinate of the nanotube, enabling us to characterise the potential of mean force (PMF) for the free energy barrier in our nanotube system. The sampling coordinate of the umbrella sampling calculations ranged between
$z = -6$
Å and
$z = 3$
Å. The
${\rm Na}^+$
and
${\rm Cl}^-$
ions were translocated along the sampling coordination, respectively, to systematically replace individual water molecules, resulting in 12 distinct sampling windows per ion species. In each window, a harmonic bias potential was applied to ensure sufficient overlap with adjacent windows. Molecular dynamics simulations of 2 ns were performed per window, and the final 1.5 ns of each trajectory was used for analysis. The unbiased PMF was reconstructed using the weighted histogram analysis method (WHAM). The theoretical foundations of umbrella sampling are detailed in Kästner (Reference Kästner2011).
3. Results
We first measured the averaged PMF for
${\rm Na}^+$
and
${\rm Cl}^-$
ions along the
$z$
-direction by using umbrella sampling (figure 1
b). The PMF increases significantly at the nanotube vicinity caused by the dielectric and dehydration energy barrier (Esfandiar et al. Reference Esfandiar, Radha, Wang, Yang, Hu, Garaj, Nair, Geim and Gopinadhan2017; Zhou & Xu Reference Zhou and Xu2020; Zhan et al. Reference Zhan, Aydin, Schwegler, Noy and Pham2020), as the basis for ion rejection. For later calculations, we numerically fitted the PMF in the form of a logistic function
$\phi _B= {\phi _0}/({1+ \textrm{e}^{-z/\delta }})$
. The energy barrier
$\phi _0$
can be estimated as
$\phi _0 = \phi _{dh} + \phi _{de}$
, as the summation of
$\phi _{dh}$
arising from ion dehydration (Esfandiar et al. Reference Esfandiar, Radha, Wang, Yang, Hu, Garaj, Nair, Geim and Gopinadhan2017) and dielectric repulsion
$\phi _{de}$
due to the water permittivity contrast from bulk to the confinement, and the
$\delta$
determined profile of the energy barrier at the nanotube entrance that we can fit from PMF in simulation.
Figure 1(c) shows the averaged flow velocity
$v_w(t)$
inside the nanotube as time passed, with a mean velocity every 10 ps.
We observed that the flow velocity alternates between high (light cyan background) values as an ‘open’ state and low values as a ‘close’ state. The ‘close’ state is attributed to a trapped ion near the nanotube entrance (see inset in figure 1
c) that is also noted in previous work (Su et al. Reference Su, Chen, Zhao and Su2019; Li et al. Reference Li, Zhang, Liu and Su2022). The trapped ion, either
${\rm Na}^+$
or
${\rm Cl}^-$
, increases the system friction, thereby reducing the flow velocity. The flow velocity increases again once the trapped ion leaves the entrance. The trapping behaviour of ions at the tube entrance is governed by the joint actions of the energy barrier and the hydrodynamic forces, where the hydrodynamic forces on hydrated ions push the ions toward the nanotube, against the nanotube energy barrier.
The time-averaged ion flux,
$j_{\textit{ion}}(t)= ({1}/{t})\sum _{i=1}^{n}n_{i}(t)$
shown as green dots in figure 1(d), remains negligible at the pressure gradient (
$\boldsymbol{\nabla }{\!P}$
) below
$2000 \,\rm bar\,nm^-{^1}$
, where
$n$
is the number of ions transported through the nanotube within a time interval
$t$
. However, ion flux increases dramatically with
$\boldsymbol{\nabla }{\!P}$
, consistent with the trend of previous experimental observations (Zhang et al. Reference Zhang, Wei, Xu and Wang2020; Biesheuvel et al. Reference Biesheuvel, Rutten, Ryzhkov, Porada and Elimelech2023; Wang et al. Reference Wang, Cao, Dykstra, Porada, Biesheuvel and Elimelech2021, Reference Wang, He, Heiranian, Fan, Song, Li and Elimelech2023) but contradicts the prediction of the SD model. The deviation in critical pressure for ion flux emergence is due to the channel size, friction coefficient and other detailed parameters. We will discuss the size effects in later texts.
The averaged water flow velocity is calculated as
$v_{w}=\langle ({1}/{n})\sum _{i=1}^{n}v_{i}(t)\rangle$
, where
$n$
is the number of water molecules in the nanotube, and
$v_{i}(t)$
is the stream velocity of water molecules measured between 0.2 nm downstream of the tube entrance and 0.2 nm upstream of the exit to minimise the effect of entrance/exit on the velocity. The velocity of water molecules in the first 100 ps was not included to eliminate the transient effects of undeveloped flow. The flow velocity exhibits a more complex behaviour than the ionic flux. For
$\boldsymbol{\nabla }{\!P}\lt 1000 \,\rm bar\,nm^-{^1}$
,
$v_w$
aligns with Sampson’s solution (dashed line in figure 1
d). However, as
$\boldsymbol{\nabla }{\!P}$
increases,
$v_w$
unexpectedly decreases, representing an NDFR effect. When
$\boldsymbol{\nabla }{\!P}$
is over
$2000 \,\rm bar\,nm^-{^1}$
,
$v_w$
rapidly increases with
$\boldsymbol{\nabla }{\!P}$
again meanwhile ionic flux emerges. The pressure gradient (
$\boldsymbol{\nabla \!}P$
) applied in NEMD simulations is orders of magnitude larger than in experiments, as it must drive sufficient water transport within a few nanoseconds despite thermal noise. Nevertheless, ionic blockade at the pore entrance persists even in larger pores under much lower pressures.
4. Theoretical analysis
4.1. Theoretical analysis for ionic flux
We employed the 1D Fokker–Planck equation (FPE) as a non-equilibrium statistic theory framework to capture the ion permeation along the
$z$
-direction:
\begin{equation} \frac {\partial W}{\partial t}=D \left [\frac {\partial }{\partial z}\left (\frac {\partial \phi _z}{k_{\textit{B}} T\partial z} W\right )+\frac {\partial ^2 W}{\partial z^2} \right ] \end{equation}
where
$W(z,t)$
,
$D$
,
$z$
,
$k_{\textit{B}} T$
are probability density of ions, diffusion coefficient of ions, displacement in the
$z$
-direction.
$k_{\textit{B}}$
is Boltzmann constant and
$T$
is temperature, respectively. Here we consider
${\rm Na}^+$
and
${\rm Cl}^-$
as identical ions, although the slightly different hydrated radius may cause little deviation. The potential field
$\phi _z=(\phi _B-\Delta \mu )-f_{h}z$
includes the energy barrier
$\phi _B$
, chemical potential
$\Delta \mu$
, and work done by hydrodynamic friction on a single ion
$f_{h}z$
with Stokes’ drag force
$f_h = 6 {\pi} \eta \sigma _{\textit{ion}} v_{w,o}$
, where
$\eta$
,
$\sigma _{\textit{ion}}$
,
$v_{w,o}$
are the water viscosity, ion’s hydrated radius and flow velocity at nanotube entrance derived from Sampson’s solution. Here we use the bulk water viscosity, since the trapped ion primarily remains in the bulk water, as shown in figure 8 in Appendix C. The Sampson solution in pure water (dashed line in figure 1
d) is validated by the agreement with the measured water velocity in the absence of ion trapping (
$v_{w,o}$
), and the flow velocity at
$\boldsymbol{\nabla \!}P \lt 1000 \,\rm bar\,nm^-{^1}$
where the probability of ion trapping is small. Here we take
${\rm Na}^{+}$
and
${\rm Cl}^{-}$
as identical particles with an averaged hydrated radius of 1.7 Åin calculations. In our system, the chemical potential
$\Delta \mu =\mu _L-\overline {\mu }= R_{a}T \,\textrm{ln} ({C_{L}}/{\bar {C}})$
is much lower than the energy barrier
$\phi _B$
and no ion permeation is observed without applying pressure,where
$R_{a}$
,
$T$
,
$C_{L}$
,
$\bar {C}$
are gas constant, temperature and concentration of the left reservoir and the system. Here
$\mu _L$
,
$\overline {\mu }$
are chemical potentials of the left reservoir and the system, respectively.
We could numerically derive the absorbing boundary (
$d'$
) by setting
$({d \phi _{z}}/{d' z})|_{d'} = 0$
. Thus, the first mean passage time (FMPT) for an ion penetrating the energy barrier and entering the nanotube can be derived as
$\tau _{t,z}$
:
\begin{equation} \tau _{t,z}=\frac {1}{D}\int _{0}^{d'} \textrm{e}^{-\frac {\phi _{z}}{k_{\textit{B}}T}}{\rm d}z\int _{z}^{d'} \textrm{e}^{\frac {\phi _{z'}}{k_{\textit{B}}T}}{\rm d}z' .\end{equation}
Here we consider a unit cell of water containing a pair of
${\rm Na}^+$
,
${\rm Cl}^-$
ions within a volume of
$(N_A C)^{-1}$
, and use it to estimate the probability that there is an ion thermodynamically moving in front of the nanotube, where
$N_A$
and
$C$
are Avogadro’s number and salt concentration. The critical time of flow in the absence of ion trapping is given by
$\tau _o = (\pi\! N_{\!{A}}\! C\! v_{w,o} R^2)^{-1}$
, where the flow velocity
$v_{w,o}$
in the nanotube is derived from Sampson’s solution
$v_{w, o}=\Delta {p}/\lambda _{\textit{nt}}$
with nanotube friction
$\lambda _{\textit{nt}} = (3\pi \eta /R)+(2L\lambda _w/R)$
accounting for the entrance (first term) and nanotube wall (second term), respectively. This
$\lambda _w$
is the friction coefficient of the nanotube wall that we measured from Green–Kubo relations.
$R$
and
$L$
are the radius of tube and length of tube respectively. Thus, the ionic flux correlates with the FMPT of ions trapped at the entrance,
$\tau _{tz}$
, and the critical time of water in the unit cell flowing through the nanotube without ion trapping,
$\tau _{o}$
:
\begin{equation} j_{\textit{ion}}=\frac {2}{\tau _o + \tau _{t,z}}. \end{equation}
Our theoretical results (solid green line in figure 1
d) agree well with the ionic fluxes from the MD simulation (green dots). Here we define a dimensionless factor
$\mathcal{P}_{\mathcal{Z}} \equiv {\tau _{o}}/({\tau _{o}+\tau _{t,z}})$
to characterise the suppression of ionic flux by the energy barrier, hence ionic fluxes can be written as
$j_{\boldsymbol{ion}}=\mathcal{P}_{\mathcal{Z}} j_{\textit{max}}$
. The flux reaches its maximum when
$\mathcal{P}_{\mathcal{Z}}$
approaching 1 (i.e. vanishing barrier), where
$j_{\textit{max}}=2\pi\! N_{\!{A}}\! C v_{w,o} R^2$
is the maximum ionic flux in the absence of an energy barrier. At low
$\boldsymbol{\nabla }{\!P}$
,
$\mathcal{P}_{\mathcal{Z}} \ll 1$
, the hydrodynamic force is insufficient to overcome the energy barrier
$\phi _B$
, resulting in negligible ionic fluxes. As
$\boldsymbol{\nabla }{\!P}$
increases, ion permeation increases significantly due to enhanced hydrodynamic force on the trapped ion.
4.2. Theoretical analysis for water velocity
Now we consider the flow velocity. We observe a high flow velocity
$v_o$
in the absence of ion trapping and a low velocity state
$v_c$
when trapped ions give additional friction in our system in figure 1(c). Simulations show that trapped ions can turn back to the reservoir instead of penetrating through, thus reopening the nanotube for water permeation. Hence, we incorporate the axial symmetric radial escape rate
$\tau _r^{-1}$
of the trapped ion. Trapping occurs due to hydrodynamic friction with radial shear stress, with maximum flow velocity along the axis (figure 2
a) (Heiranian, Taqieddin & Aluru Reference Heiranian, Taqieddin and Aluru2020). According to Langer theory (Langer Reference Langer1969; Berezhkovskii, Pollak & Zitserman Reference Berezhkovskii, Pollak and Zitserman1992; Collins et al. Reference Collins, Kramer, Carpenter, Ezra and Wiggins2014), the escape rate is determined by the saddle-point route, where particles are most likely to escape. In our system, ions experience the smallest hydrodynamic forces at the nanotube’s outer membrane boundary. Thus, we calculate the radial escape rate of a trapped ion shown as the orange dashed line in figure 2(a). Hence, the bidirectional escape rate
$\tau _{t}^{-1}$
for the trapped ion, schematically shown in figure 2(b), is the sum of both radial and axis directions:
\begin{equation} \tau _{t}^{-1} = \tau _{t,z}^{-1} + \tau _{t,r}^{-1}. \end{equation}
(a) The flow velocity distribution at the entrance of the nanotube is derived by Heiranian et al. (Reference Heiranian, Taqieddin and Aluru2020), where colour indicates the velocity from high (yellow) to low (purple). (b) Scheme of trapped ion (blue in figure 2 a) and bidirectional escape from radial and axial directions of nanotube. (c) The FMPT as a function of pressure gradient.
Thus, we need to consider the radial escape rate. By taking the hydrodynamic friction energy of the trapped ion to an ordinary FPE:
\begin{equation} \frac {\partial W}{\partial t}=D \left [\frac {\partial }{\partial r}\left (\frac {\partial \phi _r}{k_{\textit{B}}\! T\!\partial r} W\right )+\frac {\partial ^2 W}{\partial r^2} \right ]\!, \end{equation}
where
$W(r,t)$
and
$r$
are the probability density of the ions in the radial direction and displacement in the
$r$
-direction, respectively. Here
$\phi _r=f_z(r)\times r$
represents the work done by hydrodynamic friction on the trapped ion. We could calculate the FMPT
$\tau _{t,r}$
in the radial direction by defining the friction force on the trapped ion following Stokes’ laws as
$f_z(r)=f_h g(r)$
, where the first term is
$f_h = 6 \pi\! \eta \sigma _{\textit{ion}}$
. Here
$\eta$
is the water viscosity and
$\sigma _{\textit{ion}}$
is the hydrated radius of ions. The function
$g(r)$
represents the radial velocity distribution at the nanotube entrance, calculated from the stream function
$\psi (r,z)$
in cylindrical coordinates. As shown in figure 2(a),
$\psi (r,z)$
was calculated from spheroidal coordinates (
$x$
,
$m$
) according to Heiranian et al. (Reference Heiranian, Taqieddin and Aluru2020), and the details are in Appendix A. Thus the FMPT
$\tau _{t,r}$
for a specific angle
$\theta$
is calculated as
\begin{equation} \tau _{t,r,\theta }=\frac {1}{2\pi D}\int _{0}^{R}\textrm{e}^{-({\phi _{r}}/{k_{\textit{B}}T})}{\rm d}r\int _{r}^{R}\textrm{e}^{({\phi _{r'}}/{k_{\textit{B}}T})}{\rm d}r^{\prime }. \end{equation}
Considering the axisymmetric radial escape from the hydrodynamic trap, we have
\begin{equation} \tau _{t,r}=(2\pi )^{-1}\tau _{t,r,\theta }. \end{equation}
The calculated FMPT of the trapped ion is shown in figure 2(c). The FMPT
$\tau _{t,r}$
increases with
$\boldsymbol{\nabla }{\!P}$
due to enhanced hydrodynamic friction on trapped ions, which slows radial escape. In contrast,
$\tau _{t,z}$
decreases as hydrodynamic forces facilitate ion permeation. At low
$\boldsymbol{\nabla }{\!P}$
,
$\tau _t$
increases, dominated by
$\tau _{t,r}$
. As
$\boldsymbol{\nabla }{\!P}$
rises,
$\tau _t$
peaks when
$\tau _{t,z} \approx \tau _{t,r}$
, corresponding to equal ion escape rates in the radial and axial directions, leading to maximum friction at the entrance and the observed NDFR phenomenon (figure 1
d). Beyond this point,
$\tau _t$
decreases and converges with
$\tau _{t,z}$
. These results suggest that at low
$\boldsymbol{\nabla }{\!P}$
, trapped ions escape radially due to thermal fluctuations, while higher
$\boldsymbol{\nabla }{\!P}$
enhances ion permeation, reducing the ion trapping time
$\tau _t$
.
We derive the flow velocity from the velocity in the open state (
$v_{w,o}$
) and close state (
$v_{w,c}$
), along with their probabilities determined by FMPTs:
\begin{equation} v_{w} = \mathcal{P}_{\mathcal{S}} \boldsymbol{\cdot }v_{w,o}+(1-\mathcal{P}_{\mathcal{S}})\boldsymbol{\cdot }v_{w,c}. \end{equation}
The flow velocity with a trapped ion (close state) is
$v_{w, c} = \Delta p / (\lambda _{\textit{nt}} + \lambda _{\textit{ion}})$
, where the additional friction term
$\lambda _{\textit{ion}} = 6\eta \sigma _{\textit{ion}} R^{-2}$
arises from Stokes’ drag on the trapped ion. The factor
$\mathcal{P}_{\mathcal{S}} \equiv {\tau _o}/({\tau _t + \tau _o})$
represents the probability of the open state, while
$1 - \mathcal{P}_{\mathcal{S}}$
corresponds to the close state probability, where
$\tau _t$
is the bidirectional escape time determined by (4.2), and
$\tau _o$
is the critical time of flow in the absence of ion trapping.
The flow velocity (orange) and ionic fluxes (green) through a 0.33 nm radius nanotube in 0.5 M electrolyte solution. The simulation results (scatters) also demonstrated the NDFR phenomenon in the nanotube, which aligns well with theoretical calculations.
Considering the coupling between chemical potential and pressure gradient, we have the ion flux (
$j_{\textit{ion}}$
) and water flux (
$j_{w}=\rho _z v_w$
) in Onsager matrix (Mason & Lonsdale Reference Mason and Lonsdale1990; Monet et al. Reference Monet, Bocquet and Bocquet2023), where
$\rho _z=N/L$
is the linear density of water molecules along the nanotube axis with total water molecules
$N$
and nanotube length
$L$
:
\begin{equation} \begin{pmatrix} j_{\textit{ion}}\\j_{w} \end{pmatrix} = \begin{pmatrix} m_{11} & m_{12}\\ m_{21} & m_{22} \end{pmatrix} \begin{pmatrix} \Delta \mu \\ \Delta P \end{pmatrix}\!, \end{equation}
where
$m_{11}$
,
$m_{12}$
,
$m_{21}$
and
$m_{22}$
are diffusive diffusion, excess flux under pressure, diffusio-osmotic flow coefficient and water permeability under pressure, respectively. To derive the permeabilities of water and ions, we take the chemical potential
$\Delta \mu =\mu _L-\overline {\mu } =R_{a}T \,\textrm{ln}({C_{L}}/{\bar {C}})$
in the matrix. Thus, we have the ion permeation
$B = ({m_{11}}/{R_a\! T\! A})+m_{12}$
which can be simplified to
$B \simeq m_{12} = ({2 \pi\! N_{\!{A}}\! C\! R^2}/{\lambda _{\textit{nt}}}) \mathcal{P}_{\mathcal{Z}}$
, and water permeation becomes
$A = ({m_{21}}/{R_a T\! A})+m_{22}$
that can be simplified to
$A \simeq m_{22} = ({\rho _z \mathcal{P}_{\mathcal{S}}}/ {\lambda _{\textit{nt}}}) + ({\rho _z (1 - \mathcal{P}_{\mathcal{S}})}/({\lambda _{\textit{nt}} + \lambda _{\textit{ion}}}))$
, since the chemical potential is insufficient to drive transport through the nanotube in our case.
To validate our statistical theories, we approved our theory for a 0.33 nm radius nanotube in 0.5 M NaCl solution – a higher concentration than that in figure 1(d), which is a typical NaCl concentration in seawater. The results shown in figure 3 indicate our theory could also explain the influence of salt concentration under dilute conditions according to (4.3), (4.8) and
$\tau _o$
. The increased salt concentration reduces
$\tau _o$
, thus decreasing the flow velocity and increasing the ion permeability (compared with the solid line in figures 1
d and 3), in agreement with the tendency in experimental observations (Wang et al. Reference Wang, Cao, Dykstra, Porada, Biesheuvel and Elimelech2021). In addition, our theory can also capture the size effect on the flow velocity and ionic flux since the theoretical prediction agrees well with the result from the simulation for a 0.29 nm radius nanotube with 0.3 M shown in figure 9 in Appendix D. However, we neglected the ion–ion interactions in the FPE for ions penetrating the nanotube under dilute conditions, but such interactions may become significant in highly concentrated solutions or ionic liquids where ion pairs are formed in the solution.
As the energy barrier of nanotubes scales inversely with radius,
$\phi _0 \propto R^{-1}$
, according to (Parsegian Reference Parsegian1969; Richards et al. Reference Richards, Schäfer, Richards and Corry2012a
,
Reference Richards, Schäfer, Richards and Corryb
; Fumagalli et al. Reference Fumagalli2018), we can further explore the size effects in RO. By replacing
$\phi _B$
with the corresponding value for nanotubes at different sizes to our previous theories, we derived ion and water permeabilities as functions of
$\boldsymbol{\nabla }{\!P}$
and
$R$
. Our results indicate a monotonic increase in ion permeability (
$B \propto \mathcal{P}_{\mathcal{Z}}$
) with both parameters, ranging from low (yellow) to high (blue) in figure 4(a), as
$\tau _{t,z}$
decreases with increasing
$\boldsymbol{\nabla }{\!P}$
or decreasing
$R$
. When
$ \tau _{o} \ll \tau _{t,z}$
,
$B$
increases exponentially with
$\boldsymbol{\nabla }{\!P}$
, approaching 1 as
$ \tau _{o} \gg \tau _{t,z}$
and the energy barrier becomes negligible. The water permeability
$A$
in figure 4(b) shows a distinct valley (yellow), indicating the NDFR regime. We identify a threshold radius of
$R \lt 0.7$
nm for observable NDFR in our simulation, with the critical
$\boldsymbol{\nabla }{\!P}$
shifting to lower values as
$R$
increases. The NDFR behaviour vanishes when
$\lambda _{\textit{ion}}$
becomes negligible or when ion trapping probability is significantly reduced (
$\mathcal{P}_{\mathcal{Z}} \approx 1$
).
Phase diagram of permeabilities for ions (a) and water (b) as functions of pressure gradient and channel radius. Both phase diagrams share the same
$y$
-axis.
With the water permeability
$A$
and ion permeability
$B$
, we have the selectivity between water and ions as
$\alpha =A/B$
. To demonstrate the impact of NDFR on RO, we calculate the selectivity-water flux (
$\alpha$
vs
$j_w$
) curve for a 0.33 nm radius nanotube (figure 5). Without NDFR effects, selectivity decreases as water flux increases, reflecting the trade-off limit. Incorporating NDFR introduces a concave region in the trade-off curve, indicating unattainable selectivity-water flux combinations. These NDFR effects diminish as the nanotube radius increases, weakening the hydrodynamic trap for ions.
Selectivity-water flux trade-off in a single nanotube with radius of 0.33 nm. The NDFR effect strongly reduces the performance of RO in a single nanotube.
To demonstrate the NDFR is preserved in a more realistic channel, we performed molecular dynamics simulations on a
$m(10,0)$
carbon nanotube (CNT) using the settings detailed in Appendix B, with corresponding results in figure 7. Experimentally observing the NDFR effect in a single Ångström channel remains challenging due to intrinsically low flow rates, even in low-friction systems. However, unconventional approaches may overcome this limitation. For instance, the nanojet method (Secchi et al. Reference Secchi, Marbach, Niguès, Stein, Siria and Bocquet2016) can directly visualise nanoscale flows, and fluorescence correlation spectroscopy (Li et al. Reference Li, Hall, Wang, Li, Byrne, Scammell, Whitney, Allen, Cumings and Noy2024b
) could provide alternative means to characterise such low flow rates.
The NDFR phenomenon vanishes in parallel embedded nanotube membranes due to statistical averaging over a sufficiently broad channel size distribution (inset), which explains its absence in experimental observations. However, ion trapping significantly reduces the selectivity-permeability trade-off, becoming a predominant limiting factor in RO. The inset shows the probability distribution of the parallel channel sizes.
Interestingly, the NDFR phenomenon has not been reported in RO experiments, even if ion trapping could happen in practice. To investigate this discrepancy, we consider a membrane composed of parallel channels with a certain size distribution, as reported in previous experimental systems (Otero et al. Reference Otero, Mazarrasa, Villasante, Silva, Prádanos, Calvo and Hernández2008; Wang et al. Reference Wang2024). The pore size distribution is shown in the inset of figure 6. We calculate the selectivity of a membrane as
$\alpha _{m}=\sum p_{i}\alpha _{i}$
where
$p_{i}$
and
$\alpha _{i}$
are the probability at a specific channel size and the selectivity we derived in the corresponding channel size, according to previous text. The water flux is simply calculated by the summation of all channels
$j_{m,w}=\sum j_{i,w}$
where
$j_{i,w}$
is the water flux from an individual channel.
Our results in figure 6 demonstrate that NDFR vanishes in such a membrane due to statistical averaging in a certain wide channel size distribution. Since the threshold of the NDFR phenomenon is strongly correlated to the channel size, the overall effect becomes suppressed when averaging across many pores with varying diameters. This statistical averaging likely masks NDFR in experimental measurements. Here, we did not account for hydrodynamic interactions between densely packed parallel pores, which may enhance water flow permeation according to previous studies (Jensen, Valente & Stone Reference Jensen, Valente and Stone2014). Our system behaves as a membrane with sparsely distributed pores, enabling us to evaluate its performance from the response of a single channel.
Nonetheless, ion trapping at the channel entrance still leads to a substantial reduction in the selectivity-permeability trade-off, even in the absence of observable NDFR for a membrane. This ion trapping arises from dielectric exclusion as well as ion dehydration and manifests as a size-dependent energy barrier, as previously described by the FMPT framework. Therefore, ion trapping remains an intrinsic and dominant constraint on RO performance, regardless of optimal pore size or low-friction channel conditions.
5. Conclusions
In summary, we investigate the ion and water permeation through a single Ångström channel by pressure-driven flow, using a non-equilibrium statistical framework supported by all atom MD simulations. Our theory explains ion rejection at low pressure gradients due to a dielectric and dehydration energy barrier and ion permeation as pressure increases. Notably, we discovered an NDFR phenomenon, driven by ion trapping at the channel entrance, giving additional friction to the system. Our bidirectional escape model captures the NDFR effect by FMPT and provides quantitative descriptions for ion and water permeation as functions of channel radius and pressure gradient. Furthermore, we reveal a selectivity-water flux trade-off, which deviates from the classical SD model, with a statistical mechanics framework for RO. Although not all the effects like concentrations were explored here, our work provides statistical mechanics insights into RO, explaining previous observations and predicting new phenomena like NDFR.
Acknowledgements
The authors thank the anonymous referees for their helpful comments that improved the quality of the manuscript, M. Sun for the support of figure drawing, and S. Luan for the help of manuscript preparation.
Funding
This work was supported by the National Natural Science Foundation of China (grant nos 12388101 and 12241201).
Declarations of interests
The authors report no conflict of interest.
Appendix A. Calculation of stream function
$\boldsymbol{\psi (r,z)}$
Here, we derived
$\psi (r,z)$
by using the stream function
$\psi$
in spheroidal coordinates (
$x$
,
$m$
) according to Heiranian et al. (Reference Heiranian, Taqieddin and Aluru2020):
\begin{equation} \psi =\begin{cases}\dfrac {1}{2\pi }Qx^3&\textrm{for}\,b\lt x\leqslant 1,\\[6pt] \dfrac {1}{2\pi }Qx^3\left (\dfrac {2x}{b}-\dfrac {x^2}{b^2}\right )&\textrm{for}\,0\lt x\leqslant b,\end{cases} \end{equation}
where
$Q=v_{w,o}\pi R^2$
is the flow rate,
$R$
is the nanotube radius and
$b$
corresponds to
$d$
as in (A2) and (A3) in cylindrical-polar coordinates.
Following Heiranian et al. (Reference Heiranian, Taqieddin and Aluru2020), we take coordinate transformation from the confocal hyperbolas and ellipses oblate spheroidal coordinates, denoted by
$x$
and
$m$
, respectively, to the radial distance from the centre of the orifice (denoted by
$r$
) and distance from the nanotube entrance in cylindrical coordinates by using (A2) and (A3). We presented the stream function
$\psi$
in spheroidal coordinates at the nanotube entrance in both oblate spheroidal coordinates and cylindrical coordinates in figure 2(a):
\begin{align} r^2&=\big[m^2+(a+d)^2\big]\big(1-x^2\big), \end{align}
\begin{align} z^{2}&=m^{2}x^{2}. \end{align}
By substituting (A2) and (A3) into (A1), we obtain the stream function
$\psi (r,z)$
in cylindrical coordinates. The velocity distribution
$g(r)$
is then derived from the partial derivative of the stream function as
$(1/r) ({\partial \psi (r,z)}/{\partial z}) = g(r)$
. We define the absorbing boundary at the nanotube radius
$r = R$
, beyond which ions can no longer be trapped. The hydrodynamic friction energy well is given by
$\phi _r = f_h g(r) r$
.
Appendix B. NDFR effect in CNT
Here we performed additional MD simulations to validate this NDFR phenomena in carbon nanotube. Figure 7 shows velocity (orange squares) and ionic flux (green spheres) as a function of the pressure gradient in an
$m(10,0)$
carbon nanotube (CNT). Following the configuration as in the main text including the ion and water model, we conducted NEMD simulations to examine water and ionic transport through an
$m(10,0)$
carbon nanotube with 10 nm length graphene walls used as membrane to drive the flows. Following settings in Algara-Siller et al. (Reference Algara-Siller, Lehtinen, Wang, Nair, Kaiser, Wu, Geim and Grigorieva2015); Li et al. (Reference Li, Hall, Wang, Li, Byrne, Scammell, Whitney, Allen, Cumings and Noy2024a
), the atoms in graphene walls and carbon nanotubes are modelled by LJ potentials with parameter atoms with
. The water and carbon interactions are also modelled by LJ potentials with parameter 
(Cornellet al. Reference Cornell, Cieplak, Bayly, Gould, Merz, Ferguson, Spellmeyer, Fox, Caldwell and Kollman1995) and other heteronuclear LJ interactions were determined by Lorentz–Berthelot mixing rules. Following thermal equilibration at 289 K under a Nosé–Hoover thermostat and 1 atm pressure, the nanotube was filled with water molecules within 1 ns. Long-range Coulombic interactions were computed using the PPPM algorithm, with periodic boundary conditions applied in all directions. In NEMD simulations, hydrodynamic flow was induced along the longitudinal (
$z$
) axis by imposing a pressure differential
$\Delta p=P_1-P_2$
across the nanotube.
Similar as in the generic LJ systems, we also observed the NDFR phenomena of water flow and pressure-dependent ionic flux in CNT. The statistically counted ion and water fluxes were shown as yellow and green dots in figure 7. The ionic fluxes show nonlinear increases as the pressure gradient
$\Delta p$
, while the water fluxes at
$10^3$
bar nm−1 are apparently higher than the value at
$1.5 \times 10^3$
bar nm−1 representing a NDFR.
The water velocity and ionic flux as a function of pressure gradient in m10 CNT. The dashed lines are used as guidance.
Appendix C. Probability distribution of trapped ion
As shown in figure 8, the trapped ions are typically in front of the nanotube, not inside of the nanotube, that justifies using bulk viscosity and Stokes’ force in the main text.
The distribution probability of trapped ions under three different pressure gradients. The dashed line indicates the entrance of the nanotube.
Appendix D. MD simulation results in 0.29 nm radius nanotube
The simulation data (scattered points) with theoretical predictions (solid lines) across different pressure gradients. The simulation data closely align with the theoretical models. Water velocity (orange squares) matches the ion-free Sampson solution at both low and high pressure gradients, demonstrating an NDFR effect.
Figure 9 shows velocity (orange squares) and ionic flux (green spheres) as a function of the pressure gradient from the molecular dynamics simulation for an LJ tube with 0.29 nm radius. The orange solid line and green solid line calculated from the statistic theory presented in the main text agree with the simulation data (dots). Here we followed all settings as § 2. We conducted NEMD simulations to examine water and ionic transport through a single-Ångström-diameter nanotube with a radius of 0.29 nm and a length of 10 nm. The system was made from an FCC solid with a lattice parameter of 4.1 Å, while the ion and water model is the same as the one inthe main text. The LJ parameters for solid atoms defined as 
were used with Lorentz–Berthelot mixing rules for heteronuclear LJ interaction, hence the water molecules can fill the nanotube within 1 ns after thermal equilibration at 289 K under a Nosé–Hoover thermostat and 1 atm pressure. Long-range Coulomb interactions were calculated using the PPPM algorithm, and periodic boundary conditions were applied in all directions. In the NEMD simulations, the pressure drop
$\Delta p=P_1-P_2$
across the nanotube was varied from 4000 to 42 000 bar along the longitudinal
$z$
-direction by increasing the force on the leftmost membrane to drive hydrodynamic flow. For each
$\Delta p$
value, three independent simulations were performed to determine water velocities and ionic fluxes. Throughout the simulations, the outermost walls were fixed in the radial direction, but enabled to move along the
$z$
-direction. Meanwhile, the remaining walls and the nanotube vibrated around their equilibrium positions, restrained by harmonic springs with a spring constant of 
.
