Mechanistic correlation between water infiltration and framework hydrophilicity in MFI zeolites

Hydrophobic zeolites are nanoporous materials that are attracting an increasing interest, especially for catalysis, desalination, energy storage and biomedical applications. Nevertheless, a more profound understanding and control of water infiltration in their nanopores is still desirable to rationally design zeolite-based materials with tailored properties. In this work, both atomistic simulations and previous experimental data are employed to investigate water infiltration in hydrophobic MFI zeolites with different concentration of hydrophilic defects. Results show that limited concentrations of defects (e.g. 1%) induce a change in the shape of infiltration isotherms (from type-V to type-I), which denotes a sharp passage from typical hydrophobic to hydrophilic behavior. A correlation parametrized on both energy and geometric characteristics of the zeolite (infiltration model) is then adopted to interpolate the infiltration isotherms data by means of a limited number of physically-meaningful parameters. Finally, the infiltration model is combined with the water-zeolite interaction energy computed by simulations to correlate the water intrusion mechanism with the atomistic details of the zeolite crystal, such as defects concentration, distribution and hydrophilicity. The suggested methodology may allow a faster (more than one order of magnitude) and more systematic preliminary computational screening of innovative zeolite-based materials for energy storage, desalination and biomedical purposes.


Results
Fitting infiltration isotherms from experiments. It has been experimentally observed that liquid water interacts with the hydrophobic structure of pristine MFI zeolite (also known as silicalite-I) with a three-step process (see black dots in Fig. 1a) 27,50 . First, water cannot intrude the nanoporous framework at pressures lower than the infiltration one, which is typically around 90 MPa. Second, water molecules infiltrate into the silicalite-I structure within a limited range of pressures (between 90 and 110 MPa), giving rise to an endothermic effect. Third, over 110 MPa, water molecules are further compressed in the zeolite pores until the maximum framework capacity of the host structure is eventually achieved, as a consequence of steric hindrance between the molecules 1 .
Humplik and colleagues 27 experimentally characterized MFI zeolites with different concentration of hydrophilic defects, by substituting Si atoms with Al ones and thus creating silanol nests within the pristine MFI structure. On the one side, the internal structure of the pores was unaffected by the introduction of defects, as demonstrated by the similar XRD patterns for all the silicalite-I and defected MFI samples (see the Fig. 1 in ref. 27 ). On the other side, small concentrations of hydrophilic defects led to lower infiltration pressures and dramatically different shapes of infiltration isotherms, due to the alteration of the surface chemistry of the pores. Those previous experimental results are recalled in Fig. 1a (colored symbols): the infiltration isotherms of water in MFI zeolites with up to 0.5% Al/Si substitutions show a type-V shape, namely the typical behavior of hydrophobic frameworks. Under such conditions, water-water interactions are higher than water-zeolite ones and, therefore, water condensation in the MFI pores takes place by the collapse of homogeneously nucleated clusters of infiltrated water molecules 8,34,51 . A progressive transformation from type-V to type-I infiltration isotherms is observed with Al/Si substitutions larger than 0.5%. In these cases, water-zeolite interactions are dominant in the process of condensation, therefore inducing a heterogeneous (and more gradual) nucleation of solvent molecules close to the hydrophilic (i.e., defected) regions of the crystal 8,[52][53][54] . In other words, pore filling in hydrophilic MFI zeolites starts with the water vapor adsorption at pressures eventually lower than p 0 , similarly to what observed in the nanoporous materials for sorption heat storage [55][56][57] . With increasing hydrophobic behavior, instead, the pore filling pressure increases, and it becomes orders of magnitude higher than the saturated vapor pressure 13 ; under these conditions, the intruded water is liquid and pore filling occurs as an infiltration process 58 .
The Dubinin-Astakhov model (D-A) is a correlation that has been demonstrated to underpin a broad variety of adsorption processes [59][60][61] . Here, an empirical correlation similar to the D-A model is introduced for interpolating the infiltration isotherms data with a minimal number of parameters related to the characteristics of the nanoporous material, namely: where ω is the number of intruded water molecules per unit cell of nanoporous material (N/UC), p is the water pressure and T is the system temperature (k B = 1.38 × 10 −23 J K −1 ; N A = 6.022 × 10 23 mol −1 ). Analogously to the D-A model, the parameters E INF and n INF should depend on the sorbate-sorbent (water-zeolite, in this case) interaction energy and the crystal structure of the sorbent, respectively. Furthermore, the infiltration model in Eq. 1 includes also the maximum framework capacities of the adsorption (ω m = ω(p 0 )) and infiltration (ω M = ω(p M )) phase, being p M the water pressure at which ω M is eventually achieved. While ω m , ω M and p M are quantities that can be easily extrapolated from direct measures, E INF and n INF should be obtained by fitting Eq. 1 to ω − p isotherms. The XRD patterns reported in the previous work by Humplik and colleagues 27 confirm that the crystal structure of the zeolite samples in Fig. 1a can be considered as invariant in the considered range of concentrations of defects, at least as a first approximation. Therefore, it is possible to consider n INF as a quantity independent from the concentration of defects, being E INF the sole parameter affected by the increasing framework hydrophilicity. The optimization of E INF and n INF to fit the experimental results with Eq. 1 is then performed, and the best-fitting curves are reported in Fig. 1a as solid lines (R 2 > 0.95). Results show that n INF = 2 is the optimal model exponent for the considered zeolites, while the best-fitted E INF values (black dots in Fig. 1b) clearly highlight their dependence on framework hydrophilicity. In particular, the relation between E INF and the concentration of defects can be accurately fitted by a linear function (R 2 = 0.94, red dashed line in Fig. 1b), namely E INF = a 1 ⋅%Al/Si + a 2 with a 1 = 1632 J mol −1 and a 2 = 857 J mol −1 . Hence, E INF appears as a multiscale parameter that links the fundamental mechanism of water intrusion in the zeolite pores (i.e., water-zeolite nonbonded interactions) with the macroscopic, effective properties of the zeolite sample (i.e., infiltration isotherms).
Zeolite membrane simulations. Molecular dynamics simulations are then carried out to reproduce the infiltration behavior of the defected MFI zeolites observed by experiments and, consequently, to investigate the water intrusion at the atomistic scale. To this purpose, a 4 × 6 × 34 nm 3 computational domain is first simulated, www.nature.com/scientificreports www.nature.com/scientificreports/ where an MFI membrane with 4 × 6 × 4 nm 3 dimensions is placed in the middle of a water box (Fig. 2a). The adopted force field is made of bonded and nonbonded (Lennard-Jones and Coulomb) interactions, and it has been optimized in previous works 38,39 . The pores of the membrane are initially empty; increasing pressures are then applied along z axis to induce water infiltration and thus reproduce the characteristic infiltration isotherm. The hydrophobicity of the pristine MFI framework is progressively decreased by introducing silanol defects. The silanols insertion in the MFI crystal qualitatively mimics the hydrophilicity enhancement obtained in the experiments by substituting silicon atoms with aluminum ones 8,33,38,39 . Partial charges of silanol nests are set to q H = 0.45 e and q O = −0.9 e, where e is the elementary charge 8,38 .
Starting from the preliminary results reported in our previous work 38 , MFI membranes are simulated with random distributions of various concentrations of defects, which are equivalent to 0%, 0.33%, 0.89% and 3.06% substitutions of Si atoms by Al ones (%Al/Si). Coherently with experiments from the literature 35 , simulations (symbols in Fig. 3a) show that more hydrophilic membranes are characterized by lower infiltration pressures and type-I infiltration isotherms.
As previously assumed, the crystal structure of the MFI membranes can be considered as invariant in the simulated configurations. Again, Eq. 1 can be fitted to the numerical infiltration isotherms by considering that n INF is constant for the tested configurations, while E INF depends on the concentration of defects. A genetic algorithm is employed to optimize the model fitting on the results from infiltration simulations, and the best-fitted www.nature.com/scientificreports www.nature.com/scientificreports/ curves (R 2 > 0.90) are reported in Fig. 3a. The fitting procedure finds n INF = 3.14 as the most accurate parameter for the simulated zeolite membranes; the optimized values of E INF , instead, are depicted in Fig. 3b (black dots) for different defects concentration. On the one hand, n INF is larger than the value found by fitting experimental infiltration isotherms. This evidence highlights the non-ideality of the experimental structure, which may lead to discrepancies between numerical and experimental results. In fact, the experimental analysis of zeolite samples presenting pore blockage/narrowing, surface barriers or crystal contaminations may alter the accessible pore volume and thus the n INF value 39 , as also proved by the different maximum framework capacity (ω M ) obtained in the simulations (52 N/UC) with respect to the Humplik's experiments (35 N/UC) 38,50 . On the other hand, E INF is again observed to be proportional to the concentration of defects, because of the enhanced water-zeolite interactions provided by the more hydrophilic surface of nanopores. In Fig. 3b, the E INF values obtained from the MD simulations are accurately fitted (R 2 = 0.94) by a linear function (red dashed line, E INF = a 1 ⋅%Al/Si + a 2 , with a 1 = 478 J mol −1 and a 2 = 2443 J mol −1 ).
Hence, simulation and previous experimental evidence demonstrate that Eq. 1 can accurately fit the infiltration isotherms of water in MFI zeolites at varying hydrophilicity. In particular, while n INF only depends on the geometrical characteristics of the network of nanopores, E INF scales with the magnitude of the interaction potential between water and nanopores, namely zeolite hydrophilicity. Therefore, in principle, the characteristic infiltration isotherms of zeolite membranes should be predictable a priori from the fluid-crystal nonbonded interactions.

Bulk zeolite simulations.
To better investigate the mechanistic relation between water-zeolite interaction potentials and E INF (energy parameter in Eq. 1), the average nonbonded interaction energies per infiltrated water molecule are computed for MFI crystals at different pore hydration (ϑ M = ω/ω M , being ω M = 52 N/UC in the simulated cases) and concentration of defects.
To this purpose, a simulation domain containing a 10.0 × 9.9 × 5.4 nm 3 zeolite crystal is built from the unit cell of silicalite-I, with periodic boundary conditions applied along the three Cartesian axis (see Fig. 2b). Again, hydrophilic zeolites are obtained by introducing silanol nests in the pristine MFI framework, following a random distribution among the possible crystallographic sites (see Fig. 2c-e). The dry crystal is first energy minimized; then, water molecules are introduced into the zeolite pores by means of a Monte Carlo-like algorithm. The considered number of water molecules is chosen to span the whole interval of pressures studied in the infiltration experiments (Fig. 3a), that is ω = 5, 10, 30, 50 water molecules per unit cell. The mean interaction energies arising from both Coulomb and Lennard-Jones potentials are computed at equilibrium conditions.
On the one side, the water-zeolite specific interaction energy can be defined as: wz LJ wz C wz UC where N UC is the number of unit cells in the crystal (100 in the simulated cases); U LJ−wz and U C−wz are the overall water-zeolite interaction energies averaged along the simulated trajectory due to Lennard-Jones and Coulomb potentials, respectively 39 . Note that the interaction energies are computed only for water molecules completely intruded in the zeolite. E wz represents the effective nonbonded potential exerted by the surface of zeolite nanopores on each infiltrated water molecule, on average. E wz shows negative values due to the attractive nature of water-zeolite interactions within the MFI framework; however, for clarity, E wz is reported in absolute terms in the www.nature.com/scientificreports www.nature.com/scientificreports/ following analyses. In Fig. 4a, a linear dependence between E wz and the concentration of defects can be noticed, therefore denoting a clear correlation between E wz and zeolite hydrophilicity. Figure 4a also shows decreasing slopes for the E wz − %Al/Si linear relations with larger pore hydration, because of the dominating effect of water-water interactions at high pore hydration. In fact, the increase in nanopore hydration implies that a smaller fraction of the overall volume of intruded water is in contact with the nanopore surface, therefore progressively lowering the E wz value.
On the other side, the water-water specific interaction energy can be analogously defined as: ww LJ ww C ww UC being U LJ−ww and U C−ww the overall water-water interaction energies averaged along the simulated trajectory due to Lennard-Jones and Coulomb potentials, respectively 39 . Again, E ww shows negative values thus attractive interactions, but it is reported in absolute terms in the followings. Results in Fig. 4b show that the absolute water-water interaction energy tends to increase with pore hydration, mainly because of the higher number of H-bonds between intruded water molecules. In contrast to E wz , E ww appears to be almost insensible to pore hydrophilicity, especially at large hydration regimes (ϑ M → 1).

Mechanistic infiltration isotherms.
The drastically different water infiltration mechanism experimentally and numerically observed in zeolites with different framework hydrophilicity can be ascribed to the water-zeolite interactions and, therefore, E wz appears as the most suitable link between the atomistic details and overall properties of zeolite crystals. Coherently, Fig. 5    www.nature.com/scientificreports www.nature.com/scientificreports/ proportional to the water-zeolite interaction energy dictated by the concentration of defects. Note that the slope of the correlation between E INF and E wz scales with ϑ M . In fact, the E wz /E ww ratio is inversely proportional to pore hydration (see Fig. 4) and, therefore, limited absolute increments of E wz at ϑ M → 1 lead to sharper E INF increases. Hence, an accurate description of the nanoscale properties of zeolite (E wz ) is in principle enough to predict its macroscopic properties (E INF , namely infiltration isotherm of water in zeolite) by means of a multiscale correlation. For example, let us introduce a semi-empirical correlation between E INF , E wz and ϑ M , namely Summarizing, a comprehensive methodology to investigate the effect of the possible degrees of freedom (DOF, e.g. defects concentration and type, pore occlusions, etc.) of zeolite-based membranes on their water infiltration behavior can be then outlined. Here, the more general aim is providing a systematic approach for a fast computational exploration of novel nanoporous materials, with immediate applications in energy or desalination fields. The methodology can be subdivided into two distinct phases, namely (i) tuning the multiscale correlations and (ii) performing the sensitivity analyses. Figure 6a schematically depicts the first phase. In detail, a minimum set of 5 molecular dynamics or Monte Carlo simulations of the considered membrane (see the configuration in Fig. 2a) and 3 molecular dynamics runs of the corresponding bulk nanoporous crystal (see the configuration in Fig. 2b) are needed to tune the correlations allowing a systematic DOF exploration, namely ω = ω(E INF , n INF , ω M , ω m , p) (Eq. 1) and E INF = E INF (E wz , ϑ M ) (Eq. 4). The mechanistic correlation between the atomistic details of the MFI crystal and the corresponding infiltration isotherms can be subsequently determined, namely ω = ω(E wz , n INF , ω M , ω m , p). Second, Fig. 6b shows how sensitivity analyses can be then easily performed by means of a limited amount of molecular dynamics simulations, at least in the limit of small perturbations of the original setup (i.e., n INF , ω M , ω m approximately constant). Note that this hypothesis requires that the geometrical characteristics of pores (that is, zeolite framework) are not significantly altered by DOF variation. Infiltration isotherms can be finally estimated by ω ω ω ω =  = DOF i,2 ). Clearly, the procedure in Fig. 6b  www.nature.com/scientificreports www.nature.com/scientificreports/ reduction in the computational burden otherwise required for the complete procedure, which is instead depicted in Fig. 6a.

Discussion
The methodology outlined in Fig. 6 allows estimating the characteristic infiltration isotherms of water in zeolite crystals by only resorting on a few simulations. This would indeed reduce the computational burden for exploring and optimizing the degrees of freedom of the zeolite, therefore paving the way to a model-driven design of novel materials for RO or energy applications. Let us consider two exemplificative cases to test the prediction capability of the methodology reported in Fig. 6b. In these examples, the distribution and hydrophilicity of the defects in one of the MFI zeolites studied so far (0.89% Al/Si) are modified, to assess their effect on infiltration isotherms.
In the first case, the partial charge of silanol defects are changed from q H = 0.45 e and q O = −0.9 e ("weak" configuration) to q H = 0.65 e and q O = −1.1 e ("strong" configuration). According to Fig. 6b, a sole MD simulation of the zeolite bulk crystal (ϑ M = 0.95, near maximum pore hydration) is needed to compute E wz , which here takes the value of E wz = −20.24 kJ mol −1 . This value can be then used to estimate E INF and thus the complete infiltration isotherm by Eq. 4: E INF = 3.56 kJ mol −1 , which is 17% higher compared to the "weak" case (3.05 kJ mol −1 ). For validating the predicted infiltration isotherm of the MFI crystal with "strong" defects, a complete set of MD infiltration experiments is then performed over the  MPa pressure range. In agreement with the evidence from Cailliez et al. 8 , more hydrophilic zeolites (i.e. presence of stronger dipoles on the pore surface) are characterized by infiltration pressures shifted towards lower values (see Fig. 7a). Noteworthy, the infiltration isotherm directly obtained from MD data is best-fitted by E INF = 3.47 kJ mol −1 , which is only 3% lower than the value predicted by Eq. 4. Considering the workstations used to perform the abovementioned simulations (2x Dual Intel ® Xeon E5-2620v2), the methodology presented in Fig. 6b has the potential to reduce the computational burden needed to compute one infiltration isotherm by more than one order of magnitude, namely from ~5000 to ~300 CPU hours in this case.
Differently from the first case, where the silanol defects were randomly distributed in the MFI crystal (i.e., ρ def (z) = cost, being ρ def (z) the density of defects along the z axis), in the second case we analyze a Stripes Distribution of defects also at 0.89% Al/Si (SD, inset of Fig. 7a). Again, the methodology in Fig. 6b allows estimating E INF by means of a single MD run. In particular, the computed average water-zeolite energy interaction for ϑ M = 0.95 leads to E INF = 3.30 kJ mol −1 through Eq. 4. The complete set of MD infiltration simulations in Fig. 7a confirms that the SD distribution induces a slight enhancement in the zeolite hydrophilicity, with a 13% increase of E INF parameter respect to the random distribution (3.46 vs. 3.05 kJ mol −1 ). The infiltration increase given by SD defect distribution can be due to the localized enhancement of hydrophilicity provided by the defects in the central part of the framework, which promotes the creation of clusters of water molecules easing the water infiltration process (see Fig. 7b). As evident, the prediction capability of Eq. 4 is again demonstrated to be good, being only 5% the discrepancy between the predicted and actual E INF value.

conclusions
In this article, the infiltration of water into MFI zeolites characterized by different concentration of hydrophilic defects is studied by atomistic simulations validated upon previous experimental results. The introduction of defects in an initially hydrophobic MFI crystal (silicalite-I) allows controlling the hydrophilicity of the zeolite and thus the characteristic water infiltration. "weak" partial charges for silanols); black stars and line for the case with more hydrophilic defects (random distribution of defects; "strong" partial charges for silanols); orange rhombus and line for the case with SD defect distribution (stripe distribution of defects; "weak" partial charges of silanols). MD results (symbols) and optimized infiltration models (solid lines, R 2 > 0.85) are both shown. In the inset, the random (green line) and SD (orange line) 1-dimensional distribution of defects are schematically depicted. (b) 2D (y-axis averaged) density distributions of water within MFI crystals (0.89% Al/Si, ω = 10 N/UC) with different defects arrangements, namely random (upper panel) or SD (lower panel). White stars represent the position of defects in the zeolite crystal, whereas the time-averaged water density is colored from blue (lower densities) to red (higher densities). (2019) 9:18429 | https://doi.org/10.1038/s41598-019-54751-5 www.nature.com/scientificreports www.nature.com/scientificreports/ Experimental evidence from the literature showed that the characteristic water infiltration pressure in the MFI zeolites decreases with more hydrophilic pores. The resulting infiltration isotherms can be fitted by a semi-empirical infiltration model similar to the Dubinin-Astakhov's one for adsorption (see Eq. 1), considering E INF and n INF as purely tuning parameters. These previous experiments are here employed also to validate an atomistic model for water infiltration in MFI zeolite. Thanks to this validated simulation setup, this work identifies -for the first time -a mechanistic correlation between the chemical characteristics of the zeolite surface (i.e. defect concentration, distribution and type) and the E INF and n INF parameters, which are no more considered as tuning parameters but take a physical-chemical meaning. In detail, E INF is demonstrated to have a strong dependence on the interaction energy between zeolite surface and infiltrated water molecules; whereas, n INF on the geometrical structure of the zeolite.
This novel mechanistic relationship between the energy of water-zeolite interaction and the parameters of the infiltration model in Eq. 1 is finally employed to explore strategies for regulating the infiltration pressure at a given defects concentration, namely by either introducing more hydrophilic defects or tailoring their local distribution. The suggested methodology is demonstrated to be an accurate tool for reducing more than one order of magnitude the computational time needed to perform extensive sensitivity analyses on geometrical, physical and chemical degrees of freedom of zeolite crystals. The effort, hence, is to provide model-driven guidelines towards the development of advanced materials for zeolite-based devices with the possibility to accumulate, restore and dissipate mechanical energy, as well as for desalination systems based on highly permeable and selective zeolite membranes.

Methods
Molecular dynamics geometries. The framework of MFI zeolite is similar to that of both small-pore LTA (Linde Type A) and large-pore FAU (Faujasite) ones, but it has nanopores with intermediate sizes 62 . MFI has an orthorhombic crystal structure (Pnma space group), with a = 20.022 Å, b = 19.899 Å and c = 13.383 Å lattice constants 63 . Zeolites of MFI type have a 45% porosity arising from a 3-dimensional network of channels, which is given by the superimposition of both zig-zag nanopores parallel to [001] direction and straight nanopores parallel to [010] direction. The average diameter of pores is 5.6 Å, whereas channel intersections present cavities with 6.36 Å diameter.
The infiltration isotherms of water in MFI crystals are computed following the numerical protocol previously described by Fasano et al. 38 . Here, a membrane made of 2 × 3 × 3 crystal cells of MFI zeolite with dimensions 4 × 6 × 4 nm 3 is considered, with periodicity along x and y axis. The pristine MFI crystal without any defects is also known as silicalite-I, and it presents an hydrophobic behavior 50 . Inspired by the "silanol nests model" suggested by Cailliez and colleagues 8,33 , MFI membranes with growing hydrophilicity are here obtained by progressively inserting silanols in the pristine structure, with a random distribution among the possible crystallographic sites. The increased hydrophilicity of the zeolite framework provided by silanols can be related to the concentration of aluminum defects in the MFI structure: the introduction of Al atoms in silicalite-I promotes the presence of dangling oxygens, which in turn form silanol terminals in the structure. Two 4 × 6 × 15 nm 3 boxes of TIP4P water molecules under ambient conditions (T = 300 K, p = 1 bar, ρ = 1 g cm −3 , ≅30000 molecules on average) are then put in contact with the x − y planes of the dry zeolite membrane, thus obtaining the initial computational domain for the infiltration experiments. Note that the most external zeolite surface on the x − y planes is functionalized by silanols to mimic surface oxidation at the membrane-liquid interface.
Concerning the simulations of zeolite bulk crystals, the computational domain is made of 5 × 5 × 4 silicalite-I unit cells (10.0 × 9.9 × 5.4 nm 3 , periodic boundary conditions along the three Cartesian axis) to guarantee good statistics at low pore hydrations. Starting from the pristine hydrophobic framework, silanols are again progressively inserted to increase the pore hydrophilicity. Finally, zeolite pores are hydrated by a Monte Carlo-like algorithm.
Further details on the simulated geometries can be found elsewhere 38,39 .
Molecular dynamics force field. Both bonded and nonbonded interactions are modelled in the considered molecular dynamics force field. Bonded interactions take into account the chemical bonds within the zeolite framework, and are mimicked by stretch and angle harmonic potentials 64 . Nonbonded interactions are instead modelled by Coulomb and 12-6 Lennard-Jones potentials for electrostatic and van der Waals interactions, respectively. Particularly, the partial charges of silanols are tuned to fit the infiltration experiments of water in silicalite-I 33,34 , namely q Si = 1.4 e, q O = −q Si /2 and q H = q Si /4. TIP4P model 65 is used for water molecules, as also reported in previous studies about water infiltration in MFI 8 . A twin-range cut-off with shift function is used for the Lennard-Jones potentials; the Particle-Mesh Ewald algorithm with long range dispersion corrections is instead chosen for the Coulomb interactions 66 . Further details, discussions and the complete list of force field parameters are reported in previous works 8,38,39 .

Molecular dynamics protocol.
Both zeolite membranes and bulk crystals are initially energy minimized (steepest descent algorithm). Velocities of atoms are then initialized according to Maxwell distribution (300 K). The energy minimized structure is subsequently hydrated and equilibrated by means of multiple canonical (300 K; Berendsen thermostat) and isothermal-isobaric (300 K, 1 bar; Berendsen thermostat and barostat) ensembles, with up to 1.5 ns trajectories 67 . Zeolite membrane simulations are finally carried out in the isothermal-isobaric ensemble (300 K, infiltration pressure to be tested; velocity rescaling thermostat with 0.1 ps time constant 68 and isotropic Parrinello-Rahman barostat with 2 ps time constant 69 ), with water molecules progressively infiltrating through the initially empty membrane. Note that only the innermost crystal cells of the membrane are accounted for measuring water uptake, to avoid possible artifacts due to the broken crystallinity at the membrane surface 70 . Simulations are continued up to 10-35 ns, when water uptake converges to a steady state value and thus