Effect of layered water structures on the anomalous transport through nanoscale graphene channels

We analyse the enhanced flow rate of water through nano-fabricated graphene channels that has been recently observed experimentally. Using molecular dynamics simulations in channels of similar lateral dimensions as the experimental ones, our results reveal for the first time a relationship between water structure and the variation of flux in the rectangular graphene channels. The substantial enhancement in the flow rate compared to Poieseuille flow is due to the formation of layered 2D structures in the confined space, which persists up to a channel height of 2.38 nm, corresponding to six graphene layers. The structure of the water shows an intricate crystal of pentagonal and square tiles, which has not been observed before. Beyond six layers we find a sudden drop in flux due to the disordering of the water, which can be understood by classical flow dynamics.


Introduction
Water transport through nanoscale channels has fundamental importance [1] to many biological processes such as molecular binding in an aqueous environment and ion channels in bio membranes, porins and gap junctions [2][3][4][5]. It also has important industrial applications, ranging from drug delivery to gas separation [6]. Unfortunately, the molecular structures and transport properties of biological channels are extraordinarily complex and the dynamical behaviour of water flow through such confined systems is far from understood. As an example, many biological channels are ion specific and selective. They behave like gatekeepers and allow for extremely rapid transit [7] of certain ions in the confined water environment. This behaviour is dramatically different from the conventional fluid-flow theory as described by for example the Hagen-Poiseuille equation [8].
Carbon nanotubes (CNT) [9] have been used as biomimetic models of complex biological channels. These model systems are much simplified and have been analysed in substantial detail by computer simulations [10][11][12][13]. Despite the relative simple structures, understanding their flow properties remains controversial. It is generally observed that small CNT channels can lead to several orders of magnitude enhancement of flux rate, up to 4 or 5 orders of magnitude has been reported for a D=7 nm channel [10,11] with a slip length of 39-68 μm. However, Holt et al [12], who studied smaller nanotubes, found a much smaller 0.14-1.4 μm slip length for diameters of 1.3-2.0 nm. The water flow rates were 2 to 4 orders of magnitude higher than those predicted by classical models. This enhancement is drastically reduced in larger channels. It approaches the classical behaviour of liquid flow in 200-300 nm diameter CNTs [14]. Ultra large CNT experiments have also been reported by several groups and the resulting slip lengths are all less than 100 nm [15][16][17][18][19][20][21][22]. On the other hand, an experiment using ultra-long CNTs with diameters in the range 0.81-1.59 nm [23] reported flow enhancements by a factor of 882 and a slip length of about 53 nm, contradicting the results for the same diameter obtained in the previous studies. Inconsistent results were also obtained from MD simulations with periodic conditions and demonstrate enhancement of flow rates that are 2-3 orders of magnitude lower than the experimental measurements for CNTs with diameters between 0.81-4.99 nm [18,24,25]. All these contradictory results show that the diameter of the CNT channels play a decisive role in water transport in a way that is still incompletely understood. Meanwhile for water itself, an argument of hydrophobic effect towards nanosized liquid-solid contact interfaces has continued for a long time. Both x-ray and simulation studies have shown evidence that water is able to form an ordered skin of up to 12 Å in thickness [26] and thus it is expected that water will show a density depletion near the hydrophobic interfaces. However, results are inconsistent between different groups, e.g. the depletion gap has reportedly ranged from no gap to larger than 10 Å [27]. On the other hand, it has been reported that water in hydrophobic pores only occupies 60% of the pore volume, with the gap filled by a thin layer of vapor film [28]. The density of water in nanoconfinement is controversial and its transport behaviour is not well understood.
Very recently, rectangular graphene channels have been fabricated by Radha et al [29], who were able to control the dimensions of these channels in both heights and widths down to atomic precision. Water transport through these channels was found to depend sensitively on the channel height, h. For large heights, water flow obeys the classical dynamics, which predicts a linear dependence of flux on the channel width, w, as given by the Navier-Stokes equation for the rectangular channel geometry, Here r is the water density, h the viscosity, P the driving pressure, and d the slip length. However, the flux rate deviates from the classical behaviours in the nanometre regime and abnormal flux rate enhancement has been observed for narrow channels with heights between 2 and 6 graphene layers. This enhancement was attributed to the presence of an extremely high channel pressure of ∼1 Kbar (see figure 1) when fitting with equation (1) [29]. However, the microscopic understanding of the underlying mechanism for the abnormal flux enhancement is still missing. It remains unclear how nanometre dimensions affect the water flow, especially in the small height regime. Early simulations given by Liu Ling et al [30] and Liu Bo et al [31] discovered the formation of up to 3 water layers and they also gave enlightening explanations towards the fast follow and layered structures. Another work calculated the potential of mean force of graphene and also showed local energetic minima for water between graphene sheets at distances of 6.75 Å and 9.5 Å which corresponded to the formation of single and double layer of water, respectively [32]. Unfortunately, none of these studies realised the transition between convention Poiseuille flow and the fast flow. Therefore, we perform a series of MD simulations with similar geometries as used in the experiment described by Radha et al. The simulation method is similar to that mentioned in the supporting material is available online at stacks.iop.org/jpco/2/085015/mmedia of Nair et al [33].

Simulations methods
All the simulations were carried out using Gromacs version 5.0.2 [34]. The carbon water interaction was treated as an uncharged Lennard-Jones potential describing a carbon-oxygen interaction with a range of σ co =0.3367 nm, and a potential-well depth ε co =0.4247 kJ mol −1 . All the carbon atoms were held at fixed positions during the simulation. The water-water interactions were modelled using the single point charge/ extended (SPC/E) model [35]. A time step of 1 fs was adopted and the total duration of a simulation is up to 40 ns. An acceleration is added to the water molecules along the channel length direction to mimic the driving pressure, with a value of 0.0023 nm ps −1 [1]. This acceleration driven flow agrees with the pressure driven flow in terms of the velocity profile in the channel as reported before [36,37]. The temperature of the system was maintained at 300 K by using velocity rescaling with a stochastic term [38]. The electrostatic interactions and the van der Waals interaction were both calculated using a cut-off scheme with a cut-off radius of 1 nm. We simulated pressure driven water transport through 2.45 nm, 4.18 nm and infinite width channels with SPC/E model potentials without Particle Mesh Ewald (PME) term. This term is used to include the effects of periodic boundary conditions for a system that is infinite in all three dimensions; the experimental channels we simulated are finite in two directions and practically infinite along the axis of the channel (z-axis). To take into account the length of the channels in the z-direction, we applied periodic boundary conditions in all three dimensions. However, we made sure that the simulation box was larger than the cut-off length of electrostatic interactions such that the repeated channels in the x-and y-directions had no effect on one another. We felt it was inappropriate to apply PME along the z-direction, since it would give rise to physical memory effects in the flow. The results with and without PME were compared and it was shown that using PME or not dramatically effects to the flow of water in the channels as function of height. We also performed simulation with TIP4P model to compare the results with SPC/E model. This is shown in the supporting material (see figure S2). We find a direct correlation between the formation of layered water structures and the transport properties in the twodimensionally confined channels. The layered structures are similar to the ice-like structures observed under nanoscale confinement in earlier studies [39][40][41][42][43][44].

Results and discussion
The mean velocity and flux distributions for the simulated channels with different widths are shown in figure 2. As illustrated in figure 2(a), the velocity has a maximum at a channel height of 2 graphene layers and gradually decreases to up to six layers, after which it suddenly drops in the 2.46 nm and 4.18 nm width channels. However, the mean velocity distribution is much smoother for the infinitely wide channel. In figure 2(b), we see a sudden drop in the flux at the same height for the 2.46 nm and 4.18 nm wide channels, but no decrease of flux for the flow between infinite sheets. In both figures, we see a sudden increase of both velocity and flux from one layer to two-layer high channels. The different layer dependences between the mean velocity ( figure 2(a)) and the flux ( figure 2(b)) indicates that the water density in the 2D confined channels varies sensitively with the dimension of confinement, in particular with the height of the channels. To understand the sharp transition in flux properties shown in figure 2, further analysis of the structures of the water inside of the channel has been performed. Figure 3 shows a snapshot of the water structure (the carbon atoms are not shown) for channel heights ranging from 1 to 8 layers (corresponding to 6.8 Å, 10.2 Å, 13.6 Å, 18 Å, 21.4 Å, 24.8 Å, 28.2 Å and 31.6 Å respectively) for the 4.18 nm wide channel (the results are similar in 2.45 nm wide channels). As we can see, the water molecules in the channel form distinct layered structures in channels up to 6 layers high and these water layers are parallel to the graphene sheets. For a channel height of 7 and 8 layers, the water clearly shows amorphous (or liquid-like) structures. Additionally, the layer structure of water for heights from 2 to 6 layers shows AA stacking with a reasonable stability. For a single layer, water appears to form a crystal with rhombic lattice patterns. Between 2 to 6 layers, the crystal contains pentagonal and square cells, suggesting a possible new water structure in the rectangularly confined channels. Previous studies from Sun et al [45,46] experimentally and numerically predicted that uncoordinated skin water, i.e. water that has less than four neighbours at a contacting surface could form an ice like, low density phase that is hydrophobic and thermally more stable than the bulk. This layer is thought to share characteristics with supersolid, i.e. zero friction at the interface. For 7-and 8-layer high channels water looks much more like a disordered liquid, and no large-scale ordering is visible. The 2D pentagonal and square water structures shown between 2 and 6 layers have not been reported before, but they show some similarity to the structure of Ice III [47], which is a high temperature and high pressure phase of ice, however, our simulation temperature is fixed at 300 K. Additionally, pentagonal shape have also been found in 2D ice simulations by Ji Chen et al [48]. It is possible to further validate the discovery of the water density and structure in our system through experimental techniques, such as x-ray reflectivity [49] and spectroscopic ellipsometry [50]. Figure 4(a) shows the density distribution of oxygen atoms as a function of channel height (i.e. the z-direction). The numbers next to the first peaks are the mean distance of that peak from the graphene wall, i.e. the width of the depletion region. They indicate that water molecules move away from the graphene as the height of the channel increases; the distance between the water layers gradually reduces when the channel width increases. This leads to stronger layer-layer interaction of water molecules when more layers are present. This trend stops at a channel height of six graphene layers, where the depletion region (the region between water and graphene) starts to decrease. The hydrogen bonding between the water layers can no longer hold together additional layers and eventually the ordered water structures collapse, as can be seen for 7 layers where the depletion region suddenly drops to 3.4 Å and disordered liquid-like water is formed. Since the depletion regions are necessary for fast flow nanofluid [51,52], the formation of layered water structures in the channels starting from 2 layer height enhances the hydrogen bond between water and simultaneously weakens the water bonding to the graphene wall. On the other hand, for a single layer, which forms a rhombic ice structure, the water layer is not completely flat. This means water molecules must be interacting with the walls, which significantly suppresses the flow velocity and explains the velocity jump between 1 and 2 layers.
Confined water has been reported to have a liquid-liquid phase transition [53,54] which may also describe the change of phase of water between 6 and 7 layers in our simulations. The reported confinement threshold 2.2 nm is relatively close to our simulation and experiment results 2.7 nm. In bulk ice, the most common layered phase is Ice Ih, a hexagonal crystal that has a interlayer space equal to 3.678 Å [55]. This distance is comparable to the interlayer distance of graphite 3.4 Å. In our simulation, the layered water has an interlayer distance ranging from 3-3.4 Å. For every water layer, there is a small lattice mismatch about 0∼0.4 Å between water layer and graphite layer. When it goes up to 7 water layers, the accumulation of a large gap deforms the layered water structure and leads to the transition to the liquid state.
In an open (infinite width) channel, we do not observe a flow transition in the dependence on channel height of the water flux (blue line in figure 2). We again find a layered structure for one-and two-layer high channels, which explains the flux increase and velocity jump from one to two layers. However, for channels higher than two layers, water molecules no longer form any ordered structure and no flow transition happens. Thus, we conclude that to form multi-layered water structures (more than 3 layers), it is essential to have confinement in 2 directions. The enhanced flow is caused by forming an ordered ice-like layer structure of water inside the graphene channel, which makes the water like a solid with extremely low viscosity compared to a fluid.
In order to demonstrate long range charge-charge effect to the results we obtained, we also compared the results of simulations which included the PME method to those which used a cut-off scheme to handle the longrange effects [56]. The PME describes correctly the effect of periodic boundary conditions along the channel, which is lacking in our calculations above, but also replicates the channels in the other two dimensions. The dipole interactions between water molecules in different channels will be non-negligible and can lead to spurious effects.
Technically, the long-range electrostatic interaction is always of concern in periodic systems containing charged particles. Here we used the PME to calculate the long range electrostatic interaction. As shown in figure 5, the flux is very different from the previous results for the 2.45 nm wide channel. In fact, these results show a behaviour similar to our infinite width channel case ( figure 2(b) blue line). The velocity also increases rapidly from one layer to two and after then increases in a roughly linear behaviour. By considering the structure of water, we find that the layered structures form only for heights 1 and 2 layers which agrees with Yang et al [44]  (the structure for a height of 3 layers is partially disordered). Compared to the same model simulations without PME, even for the two-layer case, the water structure is not as well ordered as before, and although a layer structure can be distinguished, the regular patterns which appeared in the previous simulation are not seen in this simulation. Beyond a height of 4 layers, the layer structure completely disappears, and amorphous water is formed (see supporting material figure S4).

Conclusion
To summarize, we have investigated the phenomenon of water transport through graphene channels with variable widths and heights. We find that water can form ordered layered structures for graphene channels up to six layers in thickness. It is the physical origin for the appearance of a super nano-fluid observed in the recent experiment [29]. We have demonstrated that water confined in nano-scale rectangular channels exhibits different properties compared to bulk water and that its behavior depends strongly on the shape and size of the channel. Most importantly, the water structure is significantly altered due to the confinement at the nanometer scale, which leads to a new ice structure that does not exist in the usual phase diagram. As we show in the supporting information, these conclusions do not change with different potential models. The existence and the dynamical behavior of the new confined water structures have important implications to the understanding of model biometric channels and their transport properties in nature.