Introduction

The shortage of clean water remains one of the main issues in the global water scarcity regions. Seawater desalination has attracted extensive attention for large-scale clean water supply. Compared with conventional multi-stage flash distillation, the membrane separation technology has been proven as one of the most effective methods to produce clean water, owing to its advantages in energy conservation and high efficiency1. Common membrane technologies for desalination include reverse osmosis (RO), membrane distillation, and pervaporation (PV). RO is a mature desalination technology in industry but still suffers from several drawbacks, such as membrane fouling2, and the limited feed concentration of brine solution (3–4 wt%)3. When low-cost heat energy is available, vapor-pressure-driven membrane distillation and PV may be competitive with pressure-driven RO. PV desalination integrates excellent salt rejection (>99.5%) and low energy consumption (~2 kWh/m3)4, with separation efficiencies that are not affected by the membrane wetting issues or the salinity5.

The PV desalination membrane can preferentially permeate water molecules by simply applying a dense and hydrophilic polymer membrane, by which salt ions are rejected. The active layer of the PV desalination membrane is usually a thin (∼1 μm) polymer matrix formed on a porous support substrate by solvent evaporation method. The materials for PV desalination should exhibit hydrophilic properties to enable the preferential permeation of water, thereby facilitating the concentration of salt solutions6. Polyvinyl alcohol (PVA), known for its exceptional hydrophilicity and mechanical strength, is one of the most widely used polymer materials for fabricating PV active layers7. By modulating the PVA cross-linked network using suitable cross-linking agents, an optimal balance between membrane permeation flux and stability can be achieved8. Additionally, other polymers such as polyetherimide, cellulose derivatives, and block copolymers have also been explored in PV desalination processes5,9.

The fundamental theories of molecular transport within the active layer of PV membranes includes the pore-flow (P-F) mechanism and the widely accepted solution-diffusion (S-D) model10. In the S-D model, some assumptions have been proposed to comprehend the molecules transport process10,11,12,13: (a) a stable concentration gradient of water exists within the membrane10, (b) the polymer membrane is a nonporous dense phase11, (c) molecules are uniformly dispersed in the membrane11,12, (d) molecules motility is associated with polymer dynamics13, and (e) molecules diffuse within the polymer membrane by intermittently jumping from one free volume cavity to another13. Given the complex interaction between penetrants and membranes during the actual separation process, the S-D model has also been modified by introducing some empirical parameters. For example, Brun et al. proposed a “six coefficient model” to predict the diffusion coefficient of penetrants based on its concentration14.

Currently, PV technology is mainly applied for solvent dehydration, recovery of volatile organic compounds (VOCs) from aqueous solution, and desalination of saline solution. Molecular transport processes in membranes for solvent dehydration and VOCs recovery are consistent with the assumptions of the S-D model. The positronium annihilation lifetime spectroscopy research results demonstrated sub-nanometer pores existence in dense polymer membranes15,16, which are scarcely possible to observe by the existing electron microscopes. The swelling experiment indicated that the membrane is slightly swollen by solvent17,18,19. Moreover, our previous molecular dynamic (MD) simulation demonstrated that the size and number of these sub-nanometer pores are less affected by the dissolved molecules20. Furthermore, the S-D model-based equations also predict well in slightly swollen membrane. For instance, Mutto et al. modified the S-D models to simulate the relationship between temperature and concentration of penetrants21, which deduces a ratio of ternary feed with minimum coupling.

However, within a highly-swollen PV desalination membrane, the water molecules can easily couple with themselves to form large-sized clusters. This dispersion form fundamentally affects the separation performance of the PV membrane. The obvious fact is that the membrane fluxes for PV desalination (e.g., 10.0–100.0 kg/(m2·h) at 313 K) are great higher than that of other reports about dehydration and VOCs recovery (e.g., 0.5–2.0 kg/(m2·h) at 313 K)4,15,16,17,18. Moreover, the S-D model poses challenges in evaluating the separation performance of PV desalination, as the relationship between diffusion coefficient and the confined water concentration is experimentally difficult to measure. Interestingly, recent experiments and simulation studies have revealed that water transport in highly swollen RO membranes is not governed by the S-D model12,22, reminding us that the typical S-D model might not be entirely appropriate for highly swollen PV desalination membranes. Furthermore, the insufficient understanding of molecular transport mechanisms has become a bottleneck in the development of transport models for PV desalination.

In this study, we proposed a modeling method to build the cross-linked PVA atomistic models using four different cross-linkers. Based on those models, the structural properties of dry and swollen PVA membranes were systematically studied by the equilibrium MD simulations. The non-equilibrium MD (NEMD) simulation was also employed to study the transport mechanism of water molecules and ions through the membrane. Combined with the calculation of solvation free energy, the compensation mechanism of ion transport through the membrane was clarified. The mathematical evidence, which proves the permeation flux is still related to feed temperature by the Arrhenius equation, was illuminated. More importantly, a “near-equilibrium discrete” method was proposed to eliminate the chasm in time and space between simulations (∼1 ns, ∼1 nm) and experiments (∼1 s, ∼1 μm). Using this method, an analytical expression describing the relationship between transport diffusion coefficient and solubility was regressed from simulation results. A modified S-D model was thus established using this analytical expression, which provided the calculated membrane permeability close to the experimental value.

Results

Cross-linking produces additional pores in dry PVA membranes

The PVA chains need to be cross-linked to restrain swelling in aqueous solutions. We put forward a procedure suitable for most kinds of cross-linkers to build the cross-linked PVA atomistic model. The details of this procedure are shown in Supplementary Methods 1. Four common cross-linkers (i.e., maleic acid (MA), fumaric acid (FA), glutaraldehyde (GA), and diethoxydiphenylsilane (DEDPS)) were selected to build the atomistic models, and the molar ratio of PVA monomers to cross-linker was adjusted to control the degree of cross-linking. The molecular structures of the cross-linkers and the cross-linking reaction equations are shown in Fig. 1a and Supplementary Fig. 1, respectively.

Fig. 1: Structural properties of dry cross-linked PVA matrix.
figure 1

a Cross-linkers type. b FFV. The red dot line was the FFV of uncross-linked PVA. The error bars are present stand deviation (SD) of the FFV values for each cross-linked PVA atomistic model, and SD values obtained by three independent MD simulations. c PSD. The scatter points of simulated PSD are fitted by the Gaussian function. d Visualization of membrane pores. The colors denotes the size of pore radius. Light yellow represents the minimum value (0.3 nm) and dark blue represents the maximum value (0.6 nm). e Effect of cross-linking on HB number. The box plot was created from 1,000 data points and displays the first quartile, third quartile, whiskers (1.0×interquartile range), and median. f HB distance distribution. The lines in this figure serve as visual guides. g HB autocorrelation function and HB lifetime (inset figure). h Schematic diagram of the effect of cross-linking on the HB network of PVA. In c, d, f, g, the ratio of PVA monomer to cross-linker is 960 to 90.

Full size image

The structural properties of dry cross-linked PVA membranes were studied for comparison with the swollen membranes. Various cross-linkers cause changes in membrane density (Supplementary Fig. 2). The molecular weight distribution of cross-linked PVA membrane is shown in Supplementary Fig. 3, indicating that the chain length becomes uniform with increasing in cross-linking degree. The fractional free volume (FFV) and pores radius size distribution (PSD) are characteristics of the number of pores and the size distribution of pores, respectively. The mean and variance of PSD, and FFV increase with the increase of the degree of cross-linking (Fig. 1bc and Supplementary Fig. 4). Visualization of membrane pores enables the intuitive observation of the cross-linking effect on the PSD (Fig. 1d and Supplementary Fig. 5).

The backbone chain motion rate of PVA has a very low root mean square fluctuation (RMSF) (Supplementary Fig. 6), due to the hydrogen bond (HB) formed by hydroxyl groups on the side chain restricts the motion of the backbone chain. The sub-nanometer pores in dense polymers are voids created by the stacking of polymer chains. The solubility parameter was computed to measure the mutual interactions between PVA chains. The solubility parameter value decreases by 39.36% after removing the HB interaction (Supplementary Fig. 7), indicating that the HB interaction affects the stacking form of the PVA chains and thus affects the membrane pores.

As shown in Fig. 1e, the cross-linking would cause a significant decrease in HB number. The strength of HB between PVA chains is highly dependent on the bond length of HB and the lifetime of HB. Increasing the temperature from 303 K to 363 K can reduce the lifetime of HB from 68.74 ps to 48.32 ps, and significantly reduce the variance of the bond length distribution of HB (Supplementary Fig. 8). However, the bond length distribution (Supplementary Fig. 9 and Fig. 1f) and the lifetime of HB (Supplementary Fig. 10 and Fig. 1g) change indistinctively with the increase of the degree of cross-linking.

The equilibrium snapshots (Supplementary Fig. 11) show that the GA and DEDPS with a relatively non-polar functional group can completely block the HB between PVA chains, but the MA and FA containing carbonyl groups can form HB with PVA chains. Therefore, more additional pores in the PVA membrane can be produced using GA and DEDPS as the cross-linkers, resulting in a significant increase in FFV and a larger mean and variance of PSD (Fig. 1b–d). Collectively, the cross-linker disrupts the HB network of the PVA matrix at the cross-linking site to produce the additional pores, and the schematic diagram is shown in Fig. 1h.

Swelling resizes the membrane pores into mesoporous

Swelling in water solution is the nature of the PVA matrix. At the molecular level, the dual control volumes grand canonical molecular dynamics method23 was used to study the characteristics of membrane pores during the swelling process. The simulation details, the method to build an ultra-thin 2-dimension PVA atomistic model, and the parameters to perform swelling simulation, are shown in Supplementary Methods 2.

The snapshots of the swollen cross-linked PVA membrane (Fig. 2a) show a clear phase interface between the membrane and water solution compared to the uncross-linked PVA membrane. The corresponding density profile of the snapshots (Fig. 2b) indicated that the solubility of water molecules within the PVA membrane is uniform. In Fig. 2c, the swelling dynamic within 1000 ns shows that the swelling rate is close to zero after 200 ns. The time scale of simulation (μs) is significantly shorter than that of the experiment (days), making it challenging to obtain the degree of swelling (DS) by the simulation method in agreement with the experimental value. Therefore, a widely accepted proposal (Supplementary Methods 3) has been adopted to accelerate the swelling process24,25,26, and the DS of the membrane increased dramatically after using this proposal (Fig. 2c).

Fig. 2: MD simulation of PVA swelling in water solution.
figure 2

a Snapshot of swollen PVA. The water molecules are visualized as a cyan transparent surface. The cross-section size of those models is 4 nm × 4 nm, and the length of those models ranges from 24.0 nm to 24.4 nm. b Density profile of PVA and water. c Swelling dynamic. The left part is the result within 1000 ns, and the right part is the result adopt the proposal for accelerating the swelling process. The error bars are present SD of the DS values for the last five annealing cycle. d The HB number between PVA chains. e The HB number between PVA chains and water molecules. f PSD of swollen PVA. The scatter points of simulated PSD are fitted by the Gaussian function. g Visualization of pores with various degrees of cross-linking. The colors denotes the size of pore radius. Light yellow represents the minimum value (0.3 nm) and dark blue represents the maximum value (2.0 nm). h PSD with various degrees of cross-linking. The scatter points of simulated PSD are fitted by the Gaussian function. i DS with various degrees of cross-linking. The error bars are present SD of the DS values for the last five annealing cycle. The temperature is 300 K for all simulations. In Fig. 2cde, the lines in these figures serve as visual guide.

Full size image

As DS of MA-cross-linked PVA increases from 13.7 wt% to 43.5 wt% (Fig. 2c), the HB number between PVA chains decreases apparently from 0.13 NA/mol to 0.04 NA/mol (Fig. 2d and Supplementary Fig. 12, NA is the Avogadro constant). Before the membrane reaches a swelling equilibrium state, the water molecules can continuously form HB with PVA chains (Fig. 2e). A series of snapshots for the swelling process of the PVA membrane in water solution were captured in Supplementary Fig. 13, and the visualization of membrane pores at different times showed that the water molecules accumulated within the membrane pores can restructure the size and shape of these pores. When the swelling reaches an equilibrium state, it is found that the pore radius could be beyond 1 nm (PSD results in Fig. 2f and the visualization of pores in Supplementary Fig. 14). As a result, the size of some membrane pores transforms from microporous to mesoporous.

Although the impact of the degree of cross-linking on DS was studied using the experimental method (Supplementary Table 1), it is still quite difficult to figure out the shape and size of the swollen membrane pores. In Fig. 2g–i, a lower degree of cross-linking leads to a significant increase in the number of mesopores and DS, proving the elasticity of the cross-linked network plays an important role in restricting the growth of water clusters at the molecular level.

Confined water transitions from cluster to single molecules

A NEMD simulation was performed to understand the molecular-level mechanisms of water and ion molecule transport in the PVA membrane. As shown in Fig. 3a, the simulation box consists of a cross-linked PVA membrane (8 nm), a brine (15 nm, 3.5 wt% NaCl), one opening graphene sheet (four holes with a diameter of 1 nm), and two graphene sheets without holes. The position of those graphene sheets was fixed, and the opening graphene sheet could avoid the drift of the membrane. The temperature was set as 300 K, and the permeate side was kept vacuum to drive those molecules in the feed side through the membrane. More details of the NEMD simulation are shown in Supplementary Methods 4.

Fig. 3: NEMD simulation of water and ion transport through a cross-linked PVA membrane.
figure 3

a Setup of NEMD. Color code: membrane (cyan), opening graphene (gray), water molecules (cyan transparent surface), sodium ions (purple balls), chloride ions (green balls). b Volume fraction profile of water molecules. The grey areas represent the membrane. The inset figure shows the volume fraction profile of water molecules within membrane. The lines in this figure serve as visual guides. c Visualization of pores at different positions within the membrane. The colors denotes the size of pore radius. Light yellow represents the minimum value (0.3 nm) and dark blue represents the maximum value (1.8 nm). d PSD at different positions within the membrane. e Density profile of ions with different positions. The grey areas represent the membrane. The lines in this figure serve as visual guides. f Water numbers pass through the membrane at different times. g Schematic diagram of the difference between simulation and experiment in concentration gradient. The inset figure briefly show the concentration gradient from MD simulation.

Full size image

In Fig. 3b, the volume fraction of water in the simulation box declines quickly around the membrane surface. As the simulation time increases, the water molecules gradually show a gradient of volume fraction within the membrane. The curve of volume fraction gradient at 800 ns is similar to that at 1000 ns for cross-linked PVA with various cross-linkers (Fig. 3b and Supplementary Fig. 15). This curve remains almost unchanged until 2000 ns (Supplementary Fig. 16), confirming the assumption of the S-D model (i.e., a stable concentration gradient of water exists in the membrane).

As the volume fraction of water decreases, the dispersion form of water molecules within the membrane transforms from nano-sized clusters to single molecules (Fig. 3c). In Fig. 3d, the PSD of the membrane surface part is close to that of the equilibrium swollen membrane, and the PSD of the membrane bottom part is almost the same as the dry membrane. In Fig. 3e, the density profile of ions declines to zero around the membrane surface, indicating the saline ions are hard to penetrate into the membrane. But water molecules can leave from the bottom part of the membrane (Fig. 3f).

Our simulation revealed that the concentration gradient exists within an ultra-thin (8 nm) membrane. However, this thickness is 80∼100 folds lower than those membranes used in the experiments. Such a huge difference would result in the concentration gradient decreasing more slowly in the experiments than in the simulation (Fig. 3g). It can be noted that both the pore characteristic and DS change continuously with the concentration gradient. Therefore, we proposed a “near-equilibrium discrete” method that spatially discretizes the membrane along the concentration gradient to eliminate this chasm between simulation and experiment.

Compensatory mechanisms for ion transport via membrane pores

As shown in Fig. 4a, b, the Maxwell-Stefan (MS) diffusion coefficient of PVA backbone chains is much smaller than that of the confined water, but close to that of ions. These MS coefficients were calculated from their mean squared displacement curves (Supplementary Fig. 17), and the calculation methodology was shown in Supplementary Notes 1. The great differences in solubility (Fig. 3e) and diffusion resistance between water molecules and ions explained well the excellent desalination rate (>99.9 wt%) of the PV membrane.

Fig. 4: MD simulation of ions interaction with swollen PVA membrane.
figure 4

a MS diffusion coefficient of water and backbone chain of PVA. b MS diffusion coefficient of sodium and chloride ions. In Fig. 4ab, the error bars are present SD of MS diffusion coefficient of water and ions, and 3∼10 times independent MSD analysis were used to compute the value MS diffusion coefficient, more details were shown in Supplementary Fig. 17. c SFE of ions within swollen PVA membrane. The error bars denotes SD of ions, calculated from five independent MD simulations. d Coordination number of ions with water or PVA chains. e The HB characteristic of confined water. The value of free energy is based on the value of ions in bulk water solution, and the value of the latter is −381.32 kJ/mol and −363.54 kJ/mol for sodium and chloride ions, respectively. f RDF and coordination number of ions with hydrogen atoms on methylene of PVA chains. The DS of the PVA membrane is 2 wt%. g Electro static potential (ESP) of the compound. Those compounds consist of a molecule (water or ethanol) or an ion (sodium or chloride). The colors denotes the value of ESP. Light yellow represents the minimum value (−40.0 kJ/mol) and dark blue represents the maximum value (40.0 kJ/mol). h Schematic diagram of the compensation mechanism of ion transport in the membrane. In Fig. 4abc, the lines in these figures serve as visual guide.

Full size image

Dehydration-based energy barrier will manifest at the sub-nanometer pore inside dense polymer membranes27. The wall of membrane pores interacts directly with ions to compensate for the loss of water molecules in the hydration shell, which is well pronounced to reduce the energy barrier for ions transport through the membrane28,29. The PSD and its visualization results of PVA membrane at various DS were shown in Supplementary Fig. 18, the size and shape of membrane pores exhibited significant variations at different DS. The solvation free energy (SFE) of ions was used to measure the extent of compensation in a swollen PVA membrane (Fig. 4c). Interestingly, the SFE of chloride ions declines monotonically with DS, but the SFE curve of sodium ions has an inflection point.

For the membranes at various DS (Supplementary Fig. 19), the radial distribution functions (RDF) of both the hydroxyl groups on PVA chains around ions and the confined water molecules around ions were computed. Based on the RDF analysis, the coordination number (CN) of ions was further calculated (Supplementary Fig. 20). As shown in Fig. 4d, the total coordination number of sodium ions reduces from 5.59 to 3.75 as DS declines from 40 wt% to 2 wt%, but the total CN of chloride ions remains almost unchanged. This difference in CN results suggests that the wall of membrane pores can well compensate for chloride ions, but not for sodium ions.

When the DS declines from 40 wt% to 24 wt%, the CN of ions is almost unchanged but the SFE of ions reduces. Meanwhile, the water molecules are mainly coordinated with ions, because the membrane pores are too large for the dehydration of solvated ions. In Fig. 4e, the HB lifetime of confined water increases with the decrease of DS, which suggests that the ions and water molecules can bind more tightly by strengthening the HB between water molecules. This mechanism of strengthening the HB by reducing the size of the confined space should be the reason for the decrease in SFE.

When the DS is low (DS = 2–8 wt%), the hydration shell of ions needs to lose a part of the water molecules to adapt to the size of membrane pores, and the exposed ions will contact the wall of the membrane pores. Compared with sodium ions, the negatively charged chloride ions can coordinate more with the methylene hydrogen of the PVA chains (Fig. 4f). The weak electrostatic interactions also help the PVA chain close to chloride ions, but repulsion for sodium ions. Considering that the hydroxyl groups on ethanol and PVA have similar chemical environments, the ethanol molecule was used to study the interaction between hydroxyl groups on PVA and ions to overcome the trade-off between accuracy and computational resources. The difference in binding energy between water and ethanol to ions is only around 10 kJ/mol (Fig. 4g), indicating the hydroxyl group on the PVA chain can be well compensated for the loss of water in the hydration shell of ions. The schematic diagram of the transport mechanism of ions through the membrane is shown in Fig. 4h.

Modified S-D model for PV desalination

When the concentration gradient reaches a steady state, the equation of classical Fick’s first law was used to calculate the permeation flux:

$$J(T)=-{D}_{{\mbox{T}}}(S,T)\frac{{\mbox{d}}S}{{\mbox{d}}L}$$
(1)

where J is the permeation flux of water; DT is the transport diffusion coefficient of water molecules; S is the solubility of water molecules, L is the position along the membrane thickness; and T is the temperature. Subsequently, integrating the differential of Eq. 1 to yield Eq. 2:

$$J(T)=-\frac{1}{{L}_{{\mbox{m}}}}{\int }_{{S}_{{\mbox{f}}}}^{{S}_{{\mbox{p}}}}{D}_{{\mbox{T}}}(S,T){\mbox{d}}S$$
(2)

where Lm is the membrane thickness; Sf and Sp are the solubility of water molecules at the surface part and bottom part of the membrane, respectively. Recently, the concept of permeability originated from gas separation and has been used for evaluating the separation performance of PV membranes30,31,32,33. In PV desalination, the relationship between permeation flux and permeability is expressed as Eq. 3:

$$P(T)=\frac{1}{x\exp (A-\frac{B}{T+C})}J(T){L}_{{\mbox{m}}}=\frac{1}{x\exp (A-\frac{B}{T+C})}{\int }_{0}^{{{\mbox{S}}}_{{\mbox{f}}}}{D}_{{\mbox{T}}}(S,T){\mbox{d}}S$$
(3)

where x is the molar fraction of water in saline solution; both A, B, and C are constants of the Antoine equation for pure water solution; and P is the permeability of the water component. A detailed derivation of Eq. 3 is shown in Supplementary Notes 2.

Equation 3 can be employed to calculate permeability upon solving the definite integral. However, certain molecular-level information (e.g., diffusion modes, dispersion form) for describing molecule transport in the PVA membrane was simply encompassed within the symbol (i.e., DT) of the classical first Fick’s law, thereby making it difficult to obtain the DT of water molecules in the membrane with specific DS through experimental methods. The study of combining neutron spectroscopy and MD simulation revealed that the water molecules have various diffusion modes in highly swollen polyimide membranes34, which include a long-range diffusion of water molecules between two membrane pores and a local diffusion of water molecules inside membrane pores. We track the diffusion trajectory of water molecules within the PVA membrane (Fig. 5a, DS = 2 wt%), which shows that the diffusion of water molecules conforms to these diffusion modes.

Fig. 5: Validation of the characteristic of water transport through the membrane with reported experimental data.
figure 5

a Motion trajectory of water molecules within PVA membrane. The various colors represent the difference between RMSF. b MS diffusion coefficient of water within PVA membrane. The snapshot was captured from the equilibrium swollen PVA atomistic model (DS = 40 wt%). c The probability of long-range diffusion mode of water molecules within PVA membrane. To compare the results clarity, the logarithm of the probability was used. d Schematic diagram of the diffusion characteristics of water molecules at various DS. The grey color represents the membrane matrix. And blue arrow represents the direction of water transport. e Schematic diagram of cross-linked membrane pore types. f Schematic diagram of solving the definite integral in Eq. 3. g MS diffusion coefficient vs solubility. The scatter points are fitted by Eq. 4. h Comparison of simulated results (right part of Eq. 7) with reported experimental results (left part of Eq. 7). In Fig. 5bg, the error bars are present SD of MS diffusion coefficient of water and ions, and 3–10 times independent MSD analysis were used to compute the value MS diffusion coefficient, more details were shown in Supplementary Fig. 32. In b, h, the lines in these figures serve as visual guide.

Full size image

Considering PV desalination is a deswelling process along the transfer direction of molecule permeation, and the MS diffusion coefficient is highly affected by DS (Fig. 5b). In the swollen membranes with different DS, it is crucial to know the diffusion modes of confined water to obtain the analytical expression between the MS diffusion coefficient and solubility. The RMSF value of displacement marked every water molecule was further employed to distinguish the diffusion modes of water molecules. Unlike long-range diffusion, the local diffusion of water molecules is typically less than 1 nm. Furthermore, we programmed an in-house code (Flow-chart in Supplementary Fig. 21) to investigate the probability (PLR) of water undergoing long-range diffusion within 10 ps. The trajectory analysis of water molecules shows that 10 ps is essential to complete a long-range diffusion. Compared with the uncross-linked membrane, the cross-linked one has little impact on the value of PLR as a membrane with various DS (Supplementary Fig. 22). In Fig. 5c, PLR for water molecules in bulk, membrane (DS = 40 wt%), and membrane (DS = 2 wt%) are 0.18, 0.04, and 6.325×10-5, respectively. Such a magnitude difference in PLR reveals the reason why the MS diffusion coefficients of confined water molecules vary with DS.

Moreover, the HB autocorrelation function of confined water molecules will be gradually similar to that of bulk water as the DS increases (Supplementary Fig. 23). Moreover, the snapshot of water molecules inside the membrane was also checked (Supplementary Fig. 24). These snapshots indicate that the diffusion characteristics of water molecules inside the highly swollen membranes (DS = 24–40 wt%) are like bulk water, i.e., water molecules travel freely in water clusters. However, the diffusion characteristics of water molecules inside the lower swollen membrane (DS = 2–8 wt%) are like gas diffusion in membranes, i.e., individual water molecules jump from one pore to another pore (Fig. 5d).

Except for the effect of DS on the MS diffusion coefficient, the cross-linking also needs to be considered. As shown in Fig. 5e, the cross-linked membranes have three types of pores. The wall of pore types 1, 2, and 3 are formed by cross-linker, polymer chain, and both cross-linker and polymer chain, respectively. According to the RDF results from Supplementary Fig. 25, the water molecules keep away from the cross-linkers. Meanwhile, the cross-linking decreases the interaction energy between water molecules and PVA chains (Supplementary Fig. 26). For these pores, it is rational to deduce that the smoothness order is type 1 > type 3 > type 2 and the capability to catch water order is type 2 > type 3 > type 1. When the PVA membrane is slightly swollen, the MS diffusion coefficient of water molecules in cross-linked PVA membrane is close to that in uncross-linked membrane (Supplementary Fig. 27). However, within a highly swollen PVA membrane, cross-linking can reduce the size of water clusters by restricting the movement of the PVA chain, which significantly decreases the PLR and MS diffusion coefficient of confined water molecules.

Before an analytic expression between DMS and solubility is derived, DS needs to be transformed into solubility. The density of confined water plays a vital role in determining the relationship between DS and solubility. However, the pore structure within the PVA membrane is irregular, making it challenging to measure the occupied volume of water molecules. In Supplementary Methods 5, we propose a method for calculating the occupied volume of water molecules and derived the analytic expression between the DS and the density of confined water (Supplementary Fig. 28). Furthermore, we also derived the analytical expression between DS and solubility in Supplementary Notes 3, and this analytical expression is in good agreement (R2 = 0.99) with our simulation results (Supplementary Fig. 29).

Employing the “near-equilibrium discrete” method, the MS diffusion coefficients of the confined water within a membrane with specific solubility were computed in detail. The free volume theory revealed that the analytic expression between FFV and the diffusion coefficient of molecules within a slightly swollen polymer membrane is an exponential function35. However, this theory is insufficient to account for the diffusion of components through a highly swollen polymer membrane, as the free volume of a highly swollen polymer is experimentally challenging to measure36. Fortunately, the FFV has a linear relationship with solubility (Supplementary Fig. 30), indicating that the FFV in the free volume theory can be substituted with the solubility. Therefore, by regressing the scatter of data points, these computed MS diffusion coefficients and solubility can be related by Eq. 4:

$${D}_{{\mbox{MS}}}(S,T)={D}_{0}(T)\exp [K(T)S]$$
(4)

where D0 and K are constants related to temperature, and the numerical values of these constants at different temperatures are shown in Supplementary Table 2.

The transport diffusion coefficient is corrected to MS diffusion coefficient by Eq. 5:

$${D}_{{\mbox{T}}}(S,T)={\varGamma \left(S,T\right)D}_{{\mbox{MS}}}(S,T)$$
(5)

where the Γ is the thermodynamic factor, and this factor can be solved by using Flory-Huggins theory. The analytical expression between thermodynamic factors and solubility is shown in Eq. 6, and the proof of this equation is described in Supplementary Notes 1.

$$\varGamma (S,T)=\frac{1}{(\rho {(S,T)})^{2}}[{\rho }^{2}(S,T)-(2\chi (T)+\beta (T))\rho (S,T)S+2\chi (T){S}^{2}]$$
(6)

where χ is the Flory-Huggins parameters; and ρ is water density; β is a coefficient related to the molar volume ratio of water molecules to PVA. Combined Eqs. 3–6, the Eq. 3 can be rewritten as Eq. 7

$$J(T){L}_{{\mbox{m}}}={\int }_{0}^{{{\mbox{S}}}_{{\mbox{f}}}}{\varGamma (S,T)D}_{0}(T)\exp [K(T)S]{\mbox{d}}S$$
(7)

Herein, Eq. 7 is a bridge to connect the experiment results (i.e., \(J(T){L}_{{\mbox{m}}}\)) and simulation results (i.e., \({\int }_{0}^{{S}_{{\mbox{f}}}}{\varGamma \left(S,T\right)D}_{0}(T)\exp [K(T)S]{\mbox{d}}S\)) to avoid solving the permeability directly. Figure 5f is the schematic diagram of solving the definite integral of Eq. 7. Figure 5g shows the impacts of temperature on the MS diffusion coefficient. Kudo et al. used the swollen experiment proven the temperature has little effect on DS as the alcoholysis degree of the PVA matrix is larger than 98%37, and the 100% alcoholysis degree in this work lead to a small effect of temperature on DS (Supplementary Fig. 31). All the MS diffusion coefficients of confined water molecules were computed by mean squared displacement curves (Supplementary Fig. 32).

For membranes with different DS, the MS diffusion coefficients of confined water molecules are all consistent with the Arrhenius equation with temperature (Supplementary Fig. 33). The activation energy of confined water molecules within a highly swollen membrane is independent of solubility (Supplementary Fig. 34), and the slightly swollen region has a small contribution to the definite integral in Eq. 7 (Fig. 5f and Supplementary Fig. 35). Therefore, the permeation flux is still related to temperature with the Arrhenius equation. The mathematical proof of this conclusion is presented in detail in Supplementary Notes 4. The numerical calculation method for solving the integration in Eq. 7 is presented in detail in Supplementary Notes 5.

Finally, a fundamental S-D model, which was generated by substituting the symbol of DT in the classical Fick’s first law with an analytic expression relating the diffusion coefficient and solubility, was employed to compare the simulation results with the experimental results. As shown in Fig. 5h, the calculated permeability values are also linear with temperature and are close to the experimental values. The error mainly lies in the presence of salt crystallization on the membrane surface4,38, concentration polarization39,40, and transport barrier of support substrates41, which have adverse impacts on the permeation flux during the experiment. The numerical values and references of these experimental data are listed in Supplementary Table 3.

Discussion

This study unraveled the mechanism of water transport through the PV desalination membrane by implementing the MD simulation method and computed the water permeability via a modified S-D model. A more general method for building cross-linked PVA atomistic models was developed. Based on those models, one could find that the cross-linkers produce additional membrane pores by disrupting the HB network in the membrane. Swelling simulations revealed that water molecules could reconstruct the shape and size of the membrane pores, yielding the pore size from microporous to mesoporous. The NEMD simulation demonstrated that a concentration gradient exists within the membrane, and the dispersion form of water molecules transfers from nano-sized clusters to single molecules. According to the SFE analyses, the PVA membrane compensated the hydration layer of chloride ions more than that of sodium ions. A “near-equilibrium discrete” method was proposed to solve the analytical expression between the solubility and transport diffusion coefficient of confined water molecules. Meanwhile, this analytical expression was applied to build a modified S-D model for calculating the water permeability. This S-D model was validated with the experimental results. However, the differences in separation mechanism and mass transfer model between PV desalination and other PV membranes need further discussion, and the “near-equilibrium discrete” method would be applied in other potential fields.

It is contradicted to our simulation findings (Fig. 1bc) that some PV dehydration reports explained that the cross-linking induces decreases in both FFV and membrane pore size, and their conclusion deduced from the experimental results, in which the membrane flux decreases with the degree of cross-linking increases42,43. In these reports of PV dehydration, the membrane was slightly swollen by feed solution, and water molecules jumped from one pore to another pore in the dispersion form of single molecules or small clusters. The cross-linking resulted in the decline of the transport diffusion coefficient, thus reducing the permeation flux. However, the permeation flux of the slightly swollen PV membrane was controlled by solution rather than diffusion, and the single water molecules across the membrane were assisted by the movement of the polymer chain44. Therefore, it is a reasonable interpretation that cross-linking reduces the permeation flux by decreasing the solubility (Fig. 2i) and restricting the movement of polymer chains (Supplementary Fig. 6), rather than decreasing FFV and membrane pore size.

The solubility at a certain position within the membrane can be obtained by solving Eq. 7. The mathematical proof of this solution process is shown in Supplementary Notes 6. Recent work has mathematically proven both the S-D model and the P-F mechanism are all valid to describe the mass transfer of swollen polymer membranes45. Okada et al.46 proposed a P-F model for the prediction of PV flux, whose key presumption is a vapor-liquid interface existence within the membrane. This presumption was validated by our simulation (Fig. 3), inspiring us to derive a P-F model suitable for PV desalination in subsequent research.

MD simulation is quite difficult to directly observe the molecule transport through a realistic membrane, due to the differences between simulation and experiment in membrane thickness and operation time. In general, the membrane thickness should be reduced to about 10 nm, and some acceleration methods (e.g., higher hydraulic pressure on the membrane surface) are used to understand the transport behavior of molecules47,48,49. However, the physical parameters obtained by these methods are usually not conducive to solving the mass transfer equation (e.g., Maxwell-Stefan equation), due to the setting parameters of those simulation methods beyond the real conditions. The “near-equilibrium discrete” method aims to solve the equation by fitting a series of discrete points produced from MD simulation. The membrane thickness can be discretized spatially to solve the flux by combining the differential equation, and this method is versatile for molecular separation membrane technologies (e.g., RO, PV, gas separation). With the emergence of a more accurate force-field and improvement of MD algorithms, the “near-equilibrium discrete” method will shine in the field of molecular separation.

Methods

Simulation descriptions

All MD simulations were performed by GROMACS (version 2021) software50. The non-bond interactions were described by the Lennard-Jones (LJ) and coulomb methods, whereas the LJ interactions were truncated at 1.4 nm in conjunction with a switching function. The equations of motion were integrated by the leap-frog algorithm, and the time step was 1 fs. The temperature was maintained by the Nosé–Hoover thermostat method. The initial velocity distribution was generated by the Gaussian distribution. The Berndsen method was used to control the pressure to achieve equilibrium, and then the Parrinello-Rahman Barostat method was used for sampling.

The PVA chains were described by the OPLS-AA force-field51. The atomic charges were obtained by the RESP2 method through Multiwfn (version 3.7) software calculation52,53, and the water molecules were described by the TIP3P model54. The force-filed parameters of ions were taken from our previous work55. The fitting results of DEDPS dihedral angle parameters and the simulation density of cross-linkers are shown in Supplementary Fig. 36 and Supplementary Table 4, respectively. As shown in Supplementary Fig. 37, the density, X-ray diffractometer, and glass transition temperature of the PVA atomistic models were examined with different molecular weights. These simulation results are in agreement with the experimental results, guaranteeing the reasonableness of the force-field parameters. The force-field parameters are shown in Supplementary Tables 5-10.

Analysis methods

Pores size distribution (PSD) of all molecular models was calculated by ZEO + + (version 0.3) software56. The sample steps were 100000 and the radius of the probe was 0.01 nm. These parameters passed the convergence testing in our previous work and ensured the accuracy of the simulation results20,26. The discrete points of simulated PSD are fitted by the Gaussian function. An in-house code was converted as the output file of ZEO + + to the input file of VMD (version 1.9.3) software57, which could visualize the PSD within a certain range.

Fractional of free volume (FFV) is defined by:

$${FFV}=\left( 1 {-}1.3\frac{{V}_{{\mbox{vdw}}}}{{V}_{{\mbox{sp}}}}\right)\times 100\%$$
(8)

where Vsp is the volume of a polymer membrane and Vvdw is the van der Waals volume of a polymer matrix. The most widely used formula was reported by Lee et al.58, and the Vvdw of all molecular models in this work was calculated by Material Studio software.

HB number formed by per hydroxyl group of PVA is defined by:

$${{\mbox{HB}}}{{\mbox{number}}}=\frac{N}{{N}_{{\mbox{OH}}}-2{N}_{{\mbox{c}}}}{N}_{{\mbox{A}}}$$
(9)

where N is the total number of HB inside the PVA molecular model formed by the hydroxyl group of the PVA chain. NOH is the number of hydroxyl groups in the PVA molecular model. Nc is the number of cross-linkers. The definition of HB in this work was that the distance between donor and acceptor is smaller than 0.30 nm and the angle formed by donor, acceptor, and donor is less than 30 degrees. The characteristics of HB were recorded every 1 ps, and the HB lifetime was the integration of the HB autocorrelation.

Root mean square fluctuation (RMSF) is defined by:

$${{\mbox{RMSF}}}=\sqrt{\frac{1}{t}{\sum} _{t}\frac{1}{N}\left({\sum }_{i=1}^{i=N}{\left[{{{\boldsymbol{x}}}}_{i}\left(t\right)-{\overline{{\boldsymbol{x}}}}_{i}\right]}^{2}\right)}$$
(10)

where t and N are the simulation time and atom numbers, respectively. xi(t) is position of atom i at time t. \({\bar{x}}_{{\mbox{i}}}\) is the average position of atom i.

The Density profile was calculated by counting the number of molecules in each bin (0.2 nm). The density in each bin was attained using the mass of the molecule included in this bin to divide the total volume of this bin. The volume fraction profile and mass fraction profile also could be acquired by using the density profile.

Radial distribution functions (RDF) were used to identify the interaction sites between ions, water, and PVA, and further calculate the CN based on RDF. The CN can be obtained by:

$${N}_{{{\rm{AB}}}}(r)=4\pi \frac{{N}_{{{\rm{B}}}}}{{V}_{{{\rm{S}}}}}{\int }_{0}^{r}g(r){r}^{2}{{\rm{d}}}r$$
(11)

where r is the average cartesian distance between atom type A to atom type B, NAB(r) is the CN of atom type B to atom type A, g(r) is the RDF, NB is the total number of atom type B in simulation box, and Vs is the volume of simulation box.

Electro static potential (ESP) of ions, water, and ethanol molecules was calculated by Gaussian (version 09) software at M062X/6-311 + G(d,p) level. Analysis and visualization of ESP were performed by Multiwfn software and VMD (version 1.9.3) software, respectively. The binding energy was calculated under basis set superposition error correction.

The solvation free energy (SFE) of a single ion was calculated via the free energy perturbation (FEP) method and subsequently evaluated with the Bennett acceptance ratio (BAR) method59. The simulations were conducted by randomly positioning one ion in swollen PVA atomic models with diverse DS, and this atomic model was derived from previous diffusion simulations. Prior to the FEP calculation, a 10 ns NPT MD simulation was executed at 300 K and 1 bar to relax the simulation box, and then the CN between ions and the confined environment was calculated based on the results of RDF, where the confined environment was composed of water molecules and the functional groups of the PVA chain. Considering that the membrane pore size governs the confined environment of ions, for each DS value, at least five independent NPT MD simulations for ions at various positions were conducted to generate the initial atomistic model. In these atomistic models, the CN between ions and the confined environment must all be in close to each other, ensuring that the pore where the ions are located represents the typical confined environment of the entire simulation box. Subsequently, five independent FEP calculations were separately performed using these atomistic models to guarantee adequate sampling for each SFE value. The FEP calculation was performed in the production run for 6 ns, with the first 1 ns equilibration data being discarded. The process of simulating the solvation free energy was divided into 3-steps: (1) the neutral LJ interaction was created; (2) the electrostatic potential parameter (λcoul) was scaled from 0 to 1 with 4 replicas; (3) the van der Waals potentials parameter (λvdW) was scaled from 0 to 0.5 with 5 replicas and from 0.5 to 1.0 with 10 replicas.

The “near-equilibrium discrete” method originated from sampling and calculus concepts. When the differential equation is under the Dirichlet boundary condition, the system will attain a steady state after relaxation for a sufficiently long duration. Subsequently, a stable concentration gradient exists within the system, and Fick’s first law can be used to calculate the stable permeate flux. For each position within the system, their chemical potentials are close to those of the adjacent positions. Therefore, when we discretize the space of the macroscopic system to the nanometer level, these sub-systems should be in a near-equilibrium state. More importantly, by calculating the physical properties of these sub-systems, it is feasible to estimate the profile of the entire system. For an open system, the diffusion coefficient of the molecules is related to their concentration or position, and concentration and position can be interconverted by solving the first Fick’s law.

Employing MD simulation, the sub-system can be regarded as an atomistic model at the nanometer scale. In this study, we built a series of atomistic models consisting of various ratios of water molecules and PVA monomers, causing the DS of each atomistic model to vary. The atomistic models were built through the following steps: (1) A 10 ns NPT MD simulation at 800 K and 1 bar was conducted to relax an equilibrated dry dense PVA atomistic model to a porous state; (2) A certain amount of water molecules was inserted into the porous PVA atomistic model, and the DS of it was calculated; (3) A 10 ns MD simulation was performed at 100 bar, with the temperature decreasing from 800 K to 300 K, and the 100 bar was utilized to compress the system to prevent void generated; (4) A 10 ns NPT MD simulation was carried out at 1 bar and 300 K; (5) A 10 ns MD simulation was conducted at 100 bar. The temperature rose from 300 K to 500 K in the first 5 ns and dropped from 500 K to 300 K in the last 5 ns; (6) A 10 ns NPT MD simulation was implemented at 300 K and 1 bar. (7) Steps 5 to 6 were repeated until the density difference between the adjacent atomistic models was less than 1 kg/m3. The values of DS were selected as 2 wt%, 8 wt%, 16 wt%, 24 wt%, 32 wt%, and 40 wt%. The number of water molecules and PVA monomers are shown in Supplementary Table 11.