Multiscale simulation of water flow through laboratory-scale nanotube membranes

Water puri ﬁ cation membranes comprising aligned, dense arrays of carbon nanotubes (CNTs) have been in- vestigated for more than 10 years. Water transport 2 – 5 orders of magnitude greater than Hagen-Poiseuille predictions has been observed in CNTs of diameters 0.8 – 10nm in a small number of experiments. While the measured ﬂ ow rates in di ﬀ erent experiments substantially disagree with each other, there is a clear opportunity for these membranes to impact ﬁ ltration technologies. We propose a multiscale computational ﬂ ow method that combines molecular dynamics (MD) simulations in critical locations of the membrane with a continuum ﬂ ow resistance model. This provides the ﬂ ow resistances in a nanotube membrane con ﬁ guration to enable, for the ﬁ rst time, computationally-e ﬃ cient macroscopic predictions of ﬂ ows through laboratory-scale membranes. Our multiscale simulation results of water ﬂ ow through CNTs are also used to calibrate the Hagen – Poiseuille – Weissberg equation with slip. This study reveals that the slip length, density and viscosity can vary with CNT diameter at sub-2-nm diameters, which would otherwise be challenging to compute using MD alone. Previously published experimental results show either clear agreement or clear disagreement with our multiscale predictions; more work is required to understand this variance for similar ﬂ ow cases.


Introduction
More than a decade ago, membranes of highly-dense arrays of aligned carbon nanotubes (CNTs) were proposed for reverse osmosis water purification, in part due to their observed high permeability and their small pore sizes for ion-rejection [1,2]. CNTs with diameters D < 2 nm are narrow enough to reject the majority of larger salt ions and other contaminants through steric hindrance, while allowing water molecules to flow through at unexpectedly high rates [3]. To our knowledge, only two experiments [4,5] have recorded the water transport through CNTs with these crucial sub-2-nm diameters, although others have investigated larger nanotubes with diameters D = 3-10 nm [6][7][8][9][10][11][12][13]. These experimental reports, however, disagree on the flow enhancement, with a spread of some 2-5 orders of magnitude published for very similar CNT diameters! An effective model for flows within nanotube membranes would be a powerful tool to provide not just useful scientific insight into these experiments, but also a route to exploring a larger parametric space than is currently feasible experimentally. However, computational modelling of these types of fluid dynamics problems is challenging.
Continuum fluid models, such as the no-slip Hagen-Poiseuille flow equation, do not accurately describe the flow within sub-2-nm CNTs due to the dominance of non-continuum flow phenomena, such as high molecular ordering in the radial direction, invalid description of the local viscosity, and large slip at the internal surfaces.
Molecular dynamics (MD) is possibly the most accurate method for simulating these non-continuum flows inside nanotubes [3,[14][15][16][17][18][19][20]. Despite its inherent ability to capture the molecular physics, the major barrier to using MD to design future nanostructured membranes is the immense computational expense when modelling nanostructures, such as nanotubes, with dimensions greater than a few hundred nanometers. Membranes produced in laboratories can have nanotube lengths L from a few micrometres [4,6,10,13] to a few millimetres [5,7]. Full-domain MD simulations are computationally intractable for nanotubes of this length.
The Hagen-Poiseuille (H-P) flow equation with a Navier slip boundary condition could be applied to model flows in nanotubes with diameters down to a few nanometers, as long as MD simulations provide essential corrections to the H-P equation [21]. Recently, Walther et al. [22] (1) where p Δ is the pressure drop across the membrane, Q is the volumetric flow rate through one CNT, D D σ CO = − ⋆ is the hydrodynamic diameter of the CNT, D is the CNT diameter based on the carbon positions, σ CO is the carbon-oxygen characteristic interaction distance, L is the nanotube length (or equivalently the membrane thickness for membranes of aligned nanotubes), L s is the slip length, μ is the bulk viscosity, and C is a prefactor in Weissberg's derivation [23], which has been shown theoretically, experimentally and through simulation (including in this present work) to be close to 3.0 for short pipes or circular orifices [22][23][24]. Eq. (1) is for ideal membranes (i.e. comprising perfectly straight tubes with no defects); the inlet/outlet losses of a membrane are the first term on the right hand side, and Poiseuille slip flow is the second term.
Walther et al. [22] performed very computationally intensive CNTfilling MD simulations of a 2 nm diameter CNT, and developed a regression scheme for Eq. (1) using an additional Young-Laplace pressure term in order to determine the fitting parameters C and L s . Ritos et al. [24] performed a similar regression procedure on Eq. (1), for the same nanotube diameter, using a concurrent multiscale framework instead. Unfortunately, these simulations are still far too computationally expensive to become routine, and in any case have large numerical uncertainties.
In this paper we propose a multiscale flow method that is a more efficient alternative to the methods above, and which is also extensible to the important sub-2-nm diameter carbon nanotubes. We do not consider non-ideal membranes (i.e. defected, highly tortuous tubes), due to the absence of experimental data, but our method could in principle be used for these membranes too.
A multiscale method couples molecular dynamics with macroscopic flow theory in order to balance physical accuracy with computational cost [25]. In this paper we start by drafting a more general form of the H-P-W equation based solely on flow resistances; in this form the multiscale method can be applied to any practical membrane configuration (ideal or non-ideal). We demonstrate this multiscale method on the simple configuration considered in Walther et al. [22] and Ritos et al. [24]: a perfect CNT of D= 2 nm with entrance/exit losses. This enables us to make direct comparisons with previous computational results for validation purposes, and indicates the large computational savings that can be achieved using this new method.
We then investigate CNTs with diameters D = 0.8-4 nm, for which experimental results exist, and use the resolved flow resistances to modify the H-P-W equation with non-continuum terms for density, viscosity and slip length as a function of diameter. This enables us to develop improved models for membrane permeability and flow enhancement for use across the entire range of CNT diameters, lengths, and pressure drops, within the linear flow response regime. Fig. 1(a) illustrates part of an idealised laboratory-scale aligned nanotube (NT) membrane of thickness L, which can typically contain n ≈10 9 NTs. The modelling challenge is to predict the through-membrane mass flow rate ṁ M of water under an applied pressure drop p Δ . In our multiscale method we make the assumptions that the flow of water along one NT is steady-state, and is low Reynolds number. Under these conditions, modelling a single representative NT of the membrane ( Fig. 1(b)) is sufficient; macroscopic flow predictions on the scale of the membrane can be made by summing flow rates for all n NTs.

Methodology
Furthermore, from these assumptions the streamwise flow response within an individual NT of diameter D is likely to be linear (i.e. Stokes flow) and independent of non-continuum flow properties varying in the radial direction. Therefore the flow can be characterised using the following linear flow resistance model: where K is the flow resistance of an arbitrary NT and ṁ is the steadystate mass flow rate through the NT. In order to make accurate predictions of steady flows in general membrane configurations, the main challenge is to determine K for all nanotubes within the membrane, and therefore provide closure to Eq. (2). The flow resistance K may depend on a large number of parameters, including: nanotube length L, diameter D, fluid properties inside the NT (e.g. density, viscosity, temperature), NT material, inlet/outlet configuration (geometry, entrance/exit functionalisation), and NT structure (tapering, tortuosity, defects). A full molecular dynamics (MD) treatment of such a large range of parameters is generally intractable, mainly due to the long NTs and relatively small pressure drops that are typical in filtration.
In this paper we propose treating the flow through the NT using a multiscale method, deploying MD simulations to resolve only small but important elements of the NT, which are then coupled across space in order to solve the macro Eq. (2) and determine K. Our multiscale approach focuses on an individual nanotube ( Fig. 1(b)) and decomposes the computational domain into a number of 'components' based on scale separation in their flow physics. The flow in high-aspect-ratio NTs is highly scale separated because the flow properties (such as pressure and velocity) only vary gradually along the NT. Fabricated membranes may also have non scale separated components (e.g. entrance/exit regions, sharp bends, pinches, blockages, etc.), which could be regions where significant viscous losses occur. So, for practical membranes, the resistance flow model in Eq. (2) can be rewritten more generally to account for any number of resistance-inducing flow components present in each NT: where K i is the i th component's flow resistance, and c is the total number of resistant components in an arbitrary NT in the membrane. The flow resistances K i of all components can be resolved using MD simulations, as long as the fluid properties and the flow through all components are properly coupled together. Furthermore, the MD should not be applied over the full domain if there are to be any computational savings. Flows in components of the NT that are not scale separated should be modelled in their entirety using small MD subdomains, with appropriate boundary conditions. For the long NT components, however, any non-continuum effects mainly occur transverse to the streamwise direction and persist along the full length of the NT. It is therefore reasonable to model the flow through a long NT section using MD for small periodic slices of the NT, connected together through 1D continuity and momentum balances; this is the internal-flow multiscale method (IMM) [26,27].
While our method is for general NTs, for the purpose of proof-ofconcept, and comparison with experiments, this paper focuses on simulating ideal carbon nanotube (CNT) membranes: pristine, aligned single-walled CNTs, with a constant diameter along the CNT length. For this simple case, we decompose the CNT into c 3 = components: a perfect long CNT (treated using IMM), an entrance region and an exit region, as highlighted in Fig. 1(b). The total flow resistance K of one CNT is therefore the sum of the three component flow resistances, as we show schematically in Fig. 1(c). For simplicity (and computational efficiency) we use a single MD "subdomain 1" for the combined entrance/ exit regions. "Subdomain 2" is for the long developed-flow CNT region, which is represented by a smaller and cheaper periodic MD simulation 1 of length L′, as shown in Fig. 1(d). This enables us to simplify the flow resistance model in Eq. (3) to: where K 1 and K 2 are the flow resistances of subdomains 1 and 2, respectively, and k 2 ′ is the flow resistance per unit length measured in subdomain 2.
We propose an optimal way of running the MD simulations in multiscale methods. The concurrent multiscale method used in Ritos et al. [24] iterates ( ∼ 7 iterations) between micro/macro models until convergence; this is very wasteful of computational resource, as data that is generated from very similar MD simulations is only used once and then discarded, despite very similar flow configurations and flow rate measurements occurring in subsequent iterations. Any parameter change (e.g. a change in CNT length) requires the MD simulations to be repeated over again. It is clear that concurrent simulations are too reliant on MD simulations and so limit the parametric space for K that can be explored. Instead, by formulating the multiscale problem using flow resistances, we can run one multiscale simulation of a CNT to determine the K i 's, and these are then sufficient for subsequent macro predictions of flows through similar membranes (e.g. with longer NTs, or different pressure drops) using the macro Eq. (4). This approach is called a sequential multiscale method [28] because the resolved K i 's from the multiscale simulations can be re-used for any other similar membrane with fixed D, without running additional MD simulations. Data/information from a model at one scale feeds into parameters of another model at a larger scale. Therefore in this instance, we can also use the generated data to correct the H-P-W equation so that it incorporates any non-continuum effects occurring radially across the tube.
This multiscale approach therefore has the benefits of being able to resolve longer, complex nanotubes with smaller pressure gradients, thereby enabling quicker parametric studies and more comprehensive comparison with results measured in laboratory membranes.

Results and discussion
We use the multiscale method to simulate water flow through 11 different CNTs with tube diameters in the range D = 0.81-4 nm. All case geometries, intermolecular potentials, input parameters and measurements are given in the Appendices.

Flow through a 2 nm diameter CNT: verification and speed-up tests
To verify our multiscale method, we investigate a (15,15) CNT of diameter D 2.034 = nm. This is the same case published in Ritos et al. [24], who used different simulation approaches, including full MD simulations of CNT membranes with L ranging from 2.5 to 150 nm, and concurrent multiscale simulations with L ranging from 50 nm to 2 μm. The results in this paper were also verified using the data of Walther et al. [22].
The flow resistances obtained from our multiscale simulation are  Figs. 2(a) and (b) we plot results of the mass flow rate and pressure drop, respectivelyin the two separate subdomains of the membrane as they vary with membrane thickness. Note, it is not always possible to measure the two pressure losses (i.e. CNT and entrance/exit) separately in full MD simulations, especially for CNTs, since the watercarbon friction is usually very small, leading to very low pressure drops along the tube. However, obtaining pressure losses in individual elements of the NT membrane system is a natural output of our multiscale method and an important advantage.
There are large computational savings made by the new multiscale method. The computational time taken for this test case using 48 processes on ARCHER 2 is around 2 weeks; the greatest part of this comes from the two MD subdomain calculations. The flow behaviour for any L and any p Δ can, however, subsequently be predicted from Eq. (4) without the need for any further MD simulations. For full MD and concurrent calculations, the computational cost increases with increasing L (although at different rates [24]), and with decreasing p Δ (due to the increase in thermal noise requiring longer sampling time). Although making a fair comparison on methods is hard, if we select a typical experiment membrane of L= 2 μm, the concurrent multiscale simulation of Ritos et al. [24] for this single data point (i.e. fixed L, p Δ ) in Fig. 2 is roughly 8 times more time-consuming than what our multiscale simulation produces for any range of L and p Δ . Our method therefore realises a large saving in computational resource.

Non-continuum flow behaviour inside carbon nanotube membranes
In Fig. 3 we plot the flow resistances K 1 and k 2 ′ predicted by our multiscale simulations (red circles) for all the CNT diameters we have considered (see Appendix B). This ideal nanotube configuration may be equivalently modelled using the H-P-W Eq. (1), substituting mass flow rate for volumetric flow rate, i.e. m ρ Q = ⋆ , where ρ ⋆ is the density inside the tube as defined by the hydrodynamic diameter. This is our preferred expression for the H-P-W equation, since density may vary in sub-2-nm CNTs.
The first term on the right hand side of Eq. (1) is the pressure loss due to entrance/exit effects (which is resolved by subdomain 1 in our simulations). The second term is the loss due to flow through a tube of length L (which is partly resolved by subdomain 2). Comparing Eq. (1) with our multiscale model, Eq. (4), we retrieve continuum-like relations for the flow resistances K 1 and k 2 ′ : kg/m 3 , and viscosity μ 0.855 10 3 = × − Pa s, with C 3.0 = (which is obtained from our MD simulations). The continuum relationship for k 2 ′ is plotted using a slip length L 61 s = nm (which is that of water on graphene, from a separate MD simulation), with the bulk water properties for density and viscosity, as before.
The comparisons in Fig. 3(a) and (b) demonstrate that non-continuum flow behaviour within some CNTs goes beyond just fluid slip. Fig. 3(a) shows a jump in entrance/exit losses (K 1 ) over a small range of CNT diameters ( D 0.95 1.5 < < nm), that is not mirrored in the Weissberg term. Assuming that Weissberg's equation is still valid at this scale, the cause of this discrepancy is a non-continuum effect caused by reordering of the water molecules (i.e. a dependency on ρ ⋆ or/and μ)  (4), and the full MD simulations (empty symbols) and concurrent multiscale simulations (filled symbols) of Ritos et al. [24]. The light grey shaded region indicates 95% confidence intervals of our multiscale results. within this diameter range. To investigate this further we plot the radius-averaged density ρ ⋆ and viscosity μ against CNT diameter in Figs. 4(a) and (b). Density is measured in subdomain 1 (see Appendix A) and is seen to drop with decreasing D, due to molecular reordering in high-confinement (as we show shortly in Fig. 5).
We use two independent approaches to calculate μ; in the first, we rearrange the Weissberg term to determine viscosity as , with all other terms known from our multiscale results. In the second approach, we run independent equilibrium MD simulations of the CNT subdomains (i.e. no flow conditions), and use the Green-Kubo relationship in the streamwise direction to determine the viscosity, as described in Ref. [15]. Despite differences in the methods and geometries (due to the different subdomains in which viscosity is measured) it can be observed in Fig. 4(b) that both approaches show a spike in the fluid viscosity over the same small range of CNT diameters. Here the viscosity increases to almost twice the viscosity in the bulk. The viscosity relatively increases more than the density decreases within D 0.95 1.5 < < nm, as can be seen in the resulting kinematic viscosity (ν μ ρ / = ) in Fig. 4(c). So it is an increase in the viscosity, not a drop in overall density, that is probably the main reason for the increase in flow resistances at the inlets/outlets in this CNT diameter range.
It is still an open question as to why there is a discrepancy between the Green-Kubo measurements of viscosity and the method we propose in this paper. At molecular scales, the definitions of density and viscosity become ambiguous, but we also note that Green-Kubo may be invalid (and noisy) when applied in a region of high molecular ordering, as it is derived for homogeneous systems. In Fig. 4(b), we compare two additional studies measuring viscosity [15,21], alongside our results using the Green-Kubo relationship; they all predict very different behaviour in the viscosity, which could be caused by the different water model, the thermostat, noise, or application of the Green-Kubo equation. From here onwards in this paper we use the viscosity measurements that we inferred from our multiscale simulations. These are indicated by the red circles in Fig. 4(b).
As density and viscosity are dependent properties, the jump in viscosity stems from a local molecular re-ordering of water caused by the high confinement of the CNT. We demonstrate this in the snapshots and radial density measurements in Fig. 5. As also observed by Thomas et al. [14], for smaller CNT diameters, there is a change in the crosssectional arrangement of the constrained water molecules (i.e. from hexagon, to pentagon, diamond, circle, and single file), with more water molecules being forced into the water shell nearest the CNT surface at small diameters. The water density can reach 5 times the bulk water value in this range of small diameters, and this explains the large jump in viscosity and the pressures losses. The entrance of the membrane (rather than the exit) makes the largest contribution to losses as water molecules need to be driven from an isotropic bulk state to a highly ordered structure within a very short development length along the initial part of the nanotube [17]. Fig. 3(b) indicates that the flow resistances (k 2 ′ ) inside the CNT are much lower than H-P predictions with slip for most of the diameters considered (D 2 ≲ nm). In the literature there are some MD results showing enhanced flow at sub-2-nm diameters (e.g. see review [29]), which has been attributed to the curvature-induced smoothing of the carbon-water surface energy landscape as CNT diameters decrease [16]. Our results show very similar trends reflected in the flow resistance. Quantifying slip lengths is challenging for sub-2-nm CNT diameters because the velocity profile for such highly-confined flows is no longer parabolic. However, with the variation of density and viscosity with diameter described above, we can then use the second term on the right hand side of Eq. (5) to infer the slip length L s as it varies with CNT diameter, as all the other terms are provided by our multiscale simulations (i.e. k 2 ′ , ρ ⋆ , μ). The resultant slip lengths are plotted in Fig. 4(c), which shows that smaller diameter tubes have a much higher slip length than CNTs with diameters D < 2 nm, in which the slip length is approximately constant.

'Non-continuum' Poiseuille flow
The H-P-W Eq. (1) is a simple model but in its current form it is only applicable when D≳ 2 nm. It would be useful to have an improved version of Eq. (1) that is east to use by experimentalists or membrane engineers, but incorporates non-continuum fluid effects. We propose here a set of empirical equations for the non-continuum quantities density ρ ⋆ , viscosity μ, and slip length L s in the H-P-W equation by fitting relationships that are functions of CNT diameter to the results from our multiscale simulations. The form of these empirical relationships do not provide physical insight per se, but are proposed here to provide a practical design equation. The fitting functions are:  Figs. 3(a) and (b). Note that these equations are not applicable for other nanotube membrane materials; these would require a fresh round of multiscale simulations and data fitting.

Laboratory-scale membrane flow predictions
In publications of experimental flow measurements through CNT membranes there are currently two parameters that are used for comparison between experiments: a) the membrane permeability κ, which is used to assess membrane flow performance, and b) the nanotube flow enhancement factor E, which describes the degree of departure from classical continuum hydrodynamics within one nanotube. E is equal to the ratio of the observed flow rate in one nanotube to the calculated noslip Hagen-Poiseuille (H-P) flow rate using standard bulk properties, i.e. is the porosity (assuming all n tubes in a membrane are of constant diameter). Note that terms μ L , s and ρ ⋆ in these equations should be calculated using Eqs. (6)- (8).
As the permeability is dependent on three parameters (porosity, diameter, and length), in Table 1 we compare our predictions using Eq. (10) with a selection of experimental results. The general observation is that our predictions agree with the experiments of Holt et al. [4], Kim et al. [10] and Bui et al. [13] to within one order of magnitude, but do not agree very well with the other experiments. The results give an indication, however, of the improvement that may be possible over current commercial reverse osmosis membranes, which report values of κ≈ 1 (Ltr/m 2 -h-bar). Our defect-free CNT membrane predictions indicate there could be ∼ 1-2 orders of magnitude increase in permeability.
Holt et al. [4] also provide the distribution of CNT diameters in their membranes, which enables us to calculate a better estimate of permeability. However, our predictions differ by only 2-3% from permeability calculations that use only the mean diameter in Eq. (10). Care needs to be taken when comparing membranes with small diameters and large standard deviation, as the error in the permeability can be very sensitive. For example, a membrane with 3 nm mean nanotube diameter with a ± 1 nm standard deviation gives an error of ∼ 3% in the permeability, which rises to ∼ 10% if the standard deviation is ± 2 nm.
Inspection of Eq. (11) reveals that the flow enhancement for idealised CNT membranes depends on just two parameters: the CNT diameter D and length L. In Figs. 6(a)-(f) we compare our predictions (solid blue lines) using Eq. (11) with experimental results (symbols) for E varying with L. Every prediction is for a fixed CNT diameter D, as indicated by the arrows and values inset in each figure. For clarity, we distribute the results across six figures, from D= 0.81 nm in Fig. 6(a) to D = 10 nm in Fig. 6(f).
The variation of the flow enhancement E with increasing L is similar for all CNT diameters. When the CNTs in the membrane are short, the losses at the entrance/exit dominate, so there is no flow enhancement. For longer nanotubes, the effect of lower flow resistance inside the CNT means that the observed flow rate does not decrease as rapidly as expected by no-slip H-P theory, and so E increases steadily, and the entrance/exit effects are still important. At large L, the losses in the nanotube itself become dominant, while the inlet and outlet losses become relatively negligible, so E levels off and becomes a constant. The horizontal dashed blue lines in each figure indicate the maximum flow enhancements that can be achieved for the given CNT diameters. This behaviour can also be qualitatively understood through the two competing terms on the right of Eq. (11). The length of nanotube at which inlet/outlet losses can be neglected is given by L πDC L D (1 8 / )/(16 0.01) s ≈ + × , which is modified from [23]. Figs. 6(a) and (b) include the two experimental results for flows in sub-2-nm diameter CNTs that are important for reverse osmosis membranes. Fig. 6(a) shows the results of Qin et al. [5] for water flow through individual CNTs with D ranging from 0.81 to 1.59 nm, all of length L = 1 mm. Our predictions indicate that 1 mm long CNTs are in the constant-E region, but our results do not agree with these experiments by approximately 1-2 orders of magnitude. Fig. 6(b) shows the experimental results of Holt et al. [4] for flows through a membrane of Table 1 Membrane permeabilities κ from a selection of experiments in the literature, with our predictions via Eq. (10); the data is organised in ascending order of CNT diameter.

Experiment
Eq. (10) aligned CNTs with a distribution of diameters, but with an average D = 1.6 nm. We note that there is a difference of 1-2 orders of magnitude in the observed flow rate between Holt's and Qin's results for the same CNT diameter, D= 1.6 nm. Our results agree with Holt et al.'s [4] experiments to within one order of magnitude. In Fig. 6(c) we show the results of two experiments with D = 3.3 nm from Bui et al. [13] and Kim et al. [10]; there is very good agreement with our predictions, and between the two experiments. Finally, Figs. 6(d)-(f) show experimental results [9,6,8,7,30] for CNT diameters D = 4.8 nm, 7 nm and 10 nm, respectively. For each of these laboratory membranes, we predict much smaller flow enhancementsby approximately 3 orders of magnitude.
As most experiments investigate CNTs that are long enough to be operating at their maximum flow enhancements, E is dependent only on D for these cases. Fig. 7 shows our multiscale results for E (solid blue circles), and predictions using Eq. (11) (dashed blue line) at large L, alongside other MD results of flows through periodic CNTs, and experimental results. We divide the figure into three regimes: Regime I (not shown in the figure) where no-slip flow equations can be used (D≳ 1 μm), Regime II where fixed slip-flow can be used, and Regime III (D≲ 2 nm) where diameter-dependent slippage must be accounted for. There is a good agreement between our results and the MD simulations of Thomas et al. [14,15], as well as the experiments of Holt et al. [4], Bui et al. [13] and Kim et al. [10]. We also see good qualitative agreement with the experiments of Qin et al. [5] and Mattia et al. [31]; discrepancies in these cases could be due to the graphitic or imperfect nature of the tubes causing a drop in the slip length, and a subsequent reduction in the flow enhancement.
There is, however, very poor agreement between these results and the experiments of Baek et al.  and Zhang et al. [30], as seen by the collection of data points in the top part of Fig. 7. It is unclear why such differences exist between the various experimental results, and further investigation is needed.

Conclusions
We have proposed a computationally-efficient multiscale method that uses representative MD simulations to provide input to a macroscopic flow resistance model. This has enabled us to make, for the first time, a wide range of predictions of water transport in laboratory-scale membranes comprising carbon nanotubes, which would otherwise be too computationally expensive to perform using full MD simulations. The multiscale data we generated was then used to correct the Hagen-Poiseuille equation with Weissberg entrance/exit losses and slip by calibrating the viscosity, density and slip length. From this analysis, our recommendations for efficient membranes in terms of flow performance are to improve the inlet structure and geometry to reduce inlet/outlet losses.
Our improved H-P-W description was then compared with a range of experimental data without the need for additional computational simulations. In our comparisons of permeability and flow enhancement, experiments fell into two clear types; those that agree reasonably well with our predictions, and those that do not. More investigative work needs to be carried out from both molecular simulation and experimental viewpoints in order to resolve why this is the case.
While the H-P-W equation has its limitations (e.g. it is not applicable to non-ideal configurations), the multiscale method that we proposed in this paper can be used to model complex membrane configurations where the conventional H-P-W approach is no longer sufficient. This can help identify the selectivity or filtration capability, analyse the impact of defects, and help investigate new nanotubes of different materials that are emerging from laboratories. All MD cases of subdomain 1 (see Fig. A.1(a)) are periodic in all three directions, with the following dimensions: x 28.7 = nm, y 10.6 = nm, z 10.3 = nm. The forcing region to drive the flow is X Δ 1 = 2.52 nm, and the length of the carbon nanotube is δ 2 × = 20 nm, where δ is the approximated flow development length calculated from classical fluid dynamics. All cases of subdomain 2 (see Fig. A.1(b)) are periodic in the xdirection, with dimensions L′ indicated in Table B.3.

A.1. Subdomain 1
To enable appropriate conservation of mass between MD subdomains, the reservoir subdomain 1 simulation needs to run first in order to determine the steady-state mass density inside the CNT. The system is initialised by creating the carbon atoms of both the CNT and the graphene sheets with holes that act as the membrane surfaces, while the water molecules are initialised in the reservoirs and inside the CNT. A pressure drop p Δ 1 is applied to the system by imposing a uniform force to all water molecules located in a small forcing region at the periodic boundaries (as shown in Fig. A.1(a)) of magnitude:  where ρ n is the number density in the forcing zone, and X Δ 1 is the x-direction length of the forcing zone. A velocity-unbiased Berendsen thermostat at temperature T 298 = K is applied to the upstream and downstream reservoir regions in 6 bins in the x-direction; no thermostat is applied within the CNT.
Once the MD system reaches a steady-state in the mass flow rate, the simulation is then run again in order to set the absolute pressure in the upstream reservoir. This is achieved using the FADE algorithm [36], which inserts/deletes molecules gradually to reach the target density in the upstream reservoir of 1072 kg/m 3 , corresponding to the target absolute pressure of 200 MPa. Once the system pressure is reached, the density control procedure is turned off and the system is then run one last time in order to measure both the steady-state mass flow rate ṁ 1 at a central cross-section plane through the CNT, and the average mass density ρ 1 in the middle CNT region (see Fig. A.1(a)). The latter is required as an initial condition for subdomain 2 (see below), since water confined inside sub-2-nm CNTs no longer has a bulk-like density [37].
For a CNT diameter D, we evaluate K 1 as the gradient of the variation of pressure drop p Δ 1 with mass flow rate ṁ 1 that is assumed to be linear, i.e.

K d p dm
(Δ ) . The central CNT subdomain 2 simulation is run next. We initialise the system by creating the carbon atoms in the CNT element of length L′ with the same diameter as in subdomain 1. We choose as large L′ as practical in order to improve the simulation statistics. Then we fill the system with water molecules to match the density measured in subdomain 1, i.e. ρ 1 ; this ensures the two subdomains are coupled in terms of water density inside the CNT. A pressure gradient is then imposed on the system by applying a force to all water molecules, of magnitude: that excludes a portion of the gap between the carbon and water molecules [22,24]. In this paper we assume the hydrodynamic diameter to be D D σ OC = − ⋆ , where D is the CNT diameter (carbon-to-carbon distance) and σ OC = 0.319 nm is the oxygen-carbon characteristic length scale in the Lennard-Jones model. Similar to before, we choose an appropriate value of p Δ 2 ′ that is within the linear flow response range, but large enough to produce good statistics. A velocity-unbiased Berendsen thermostat at temperature T 298 = K is also applied to the entire system. Once the system reaches a steady-state, the flow rate ṁ 2 is measured and k 2 ′ calculated from the gradient of the pressure drop vs mass flow rate graph, i.e.