# Molecular Study of Nonequilibrium Transport Mechanism for Proton and Water in Porous Proton Exchange Membranes

Authors: Gang Wang; Zhi Liu; Chang Liu; Wei Chen
Journal: International Journal of Energy Research
 (2023)

---

## Abstract

The transport process of proton and water through a porous proton exchange membrane (PEM) under pressure difference is significant to fuel cell performance. Nonequilibrium molecular dynamics simulations were conducted to investigate the diffusion mechanisms of hydronium ions and water molecules under pressure differences. Here, different driving forces (0.15 kcal/(mol·Å)–0.45 kcal/(mol·Å)) were applied to hydronium ions and water molecules as they transported through Nafion 117 membrane at temperatures between 300 K and 350 K. Results indicated that the transport diffusion coefficient of water molecules was larger than that of hydronium ions under different pressure drops due to the lower molecular weight and diffusion activation energy of water molecules (20.25 kJ/mol). Hydronium ions formed stronger hydrogen bonds with sulfonic acid groups than water molecules, which resulted in a higher diffusion activation energy for hydronium ions (21.15 kJ/mol). The proton conductivity rose with the increase in pressure difference. Moreover, the transport diffusion coefficient of water molecules and hydronium ions were positively correlated with temperature. This is because the kinetic energy of molecules increased and the pore size of the Nafion membrane enlarged at high temperatures. In addition, the dynamics of the hydrogen bond between hydronium ions and sulfonic acid groups/water molecules were accelerated at elevated temperatures, which further promoted proton transfer.


---

## Body

## 1. Introduction

Proton exchange membrane fuel cells (PEMFCs) are extensively used as promising substitutes for fossil energy due to their high power density and low pollution [1–3]. They have great potential in various applications, especially vehicle propulsion systems [4]. As a key component in the PEMFC system, the proton exchange membrane (PEM) plays a critical role in conducting protons and separating electrodes [5]. However, the development of PEMFCs with high durability and reliability in diverse operating conditions is the greatest barrier to commercializing PEMFCs [2]. To achieve a better design of PEMFCs, it is of importance to understand the fundamental physical processes governing water and proton transport through PEMs under different conditions (e.g., pressure gradients and high temperatures).Generally, the Grotthuss mechanism and Vehicle mechanism are the two main mechanisms for proton transport in PEMs (Figure1) [6, 7]. In the Grotthuss mechanism, protons hop from a hydrolyzed ionic site (SO3- and H3O+) to an adjacent one via hydrogen-bonded network chains [8]. In the Vehicle mechanism, protons diffuse from anode to cathode in the form of hydrated hydrogen ions (H+(H2O)x) [9, 10]. The transport of water molecules in PEMs can be induced by the difference in electrochemical (electroosmotic drag), concentration (diffusion), and pressure (hydraulic permeation) between the anode and cathode sides (Figure 1) [11–13]. For ambient pressure PEMFCs, the effect of hydraulic permeation is negligible due to the low permeability of PEMs [14]. Nevertheless, for some PEMFC applications with high back pressure, elevated operating pressure and temperature can significantly affect water balance and proton conductivity in PEMs, and the transport diffusion of proton and water through PEMs driven by pressure gradients cannot be neglected [15, 16]. In these cases, an enhanced pressure difference between anode and cathode is beneficial for the output performance of PEMFCs under the premise of ensuring the safety of PEMs because it can alleviate water flooding in the cathode [17–21]. Moreover, reasonably high temperature can promote the mass transfer of reactants and proton conductivity [2, 22]. Therefore, it is essential to investigate the transport processes of proton and water under different pressure differences and temperatures.Figure 1
Proton and water transport mechanisms in a PEM. GDL: gas diffusion layer; MPL: microporous layer; ACL: anode catalyst layer; CCL: cathode catalyst layer.Molecular dynamics (MD) simulation is a feasible tool for understanding the microscopic mechanism of water and proton transport through PEMs under various working conditions. Various factors (e.g., temperature, water content, and the morphology of PEMs) affect the transport of water and protons through PEMs. By conducting equilibrium MD simulations, Venkatnathan et al. and Chen et al. all reported that the self-diffusion coefficients of hydronium ions and water molecules in Nafion membranes had a positive correlation with temperature because molecular mobility increased and the void of PEMs enlarged at high temperatures [23, 24]. Cui et al. found that the self-diffusion coefficient of hydronium ions reduced with decreasing water content because low water network connectivity led to a reduction in proton transport channels [25]. Moreover, water molecules in the Nafion membrane with a blocky sequence had a higher self-diffusion coefficient than those in the dispersed sequence because of the large hydrophilic clusters and wide water channels in the blocky case [26]. In addition, hydrogen bond interaction has a significant role in proton transport. Simulation results of Zhang et al. indicated that hydrated hydronium complexes were formed due to the hydrogen bond interaction between hydronium ions and water molecules [27]. Paddison et al. found that the nature of hydrogen bonds between water molecules and sulfonic acid groups affected proton transfer [28]. Hydrogen bond arrangement in the carbon nanotube decorated with fluorine atoms and sulfonate groups was analyzed by Song et al., and the result showed that fluorination stabilized the connectivity of hydrogen bonds to maintain continuous proton transport [29].These simulations mentioned above only study the self-diffusion of small molecules in PEMs at equilibrium. However, the transport of water and proton through PEMs is a nonequilibrium process due to the difference in electrochemical, concentration, or pressure between anode and cathode [30]. The transport diffusion of small molecules in porous media can be studied by various simulation techniques of nonequilibrium molecular dynamics (NEMD) such as dual control volume grand canonical MD (DCV-GCMD), boundary-driven NEMD (BD-NEMD), and external field NEMD (EF-NEMD) [31, 32]. Previous studies have shown that these techniques yield the same transport diffusion coefficient of small molecules in nanopores [31, 32]. Yan et al. conducted NEMD simulations in the presence of electric fields to investigate the electroosmotic drag in PEMs, and they found the average transport velocity of water molecules was lower than that of hydronium ions [33]. Li et al. analyzed the interfacial transport characteristics of water molecules at the Nafion-vacuum interface under water concentration difference, and the results of NEMD simulations indicated that interfacial transport coefficients weakly depended on the water contents for high-hydration levels [34]. By imposing external driving forces on water molecules in NEMD simulations, Daly et al. found the mass-transfer resistance at PEM/water interface was negligible [35]. Yet, there are limited studies of water and proton transfer through PEMs under pressure gradients (i.e., driving forces).In this study, EF-NEMD simulations were conducted to elucidate the transport diffusion properties of hydronium ions and water molecules through Nafion membranes at different temperatures and pressure differences (i.e., driving forces). Firstly, the models of the Nafion membrane under varied temperatures (300 K, 315 K, 330 K, and 350 K) were established. Then, the radial distribution function (RDF), the number and strength of hydrogen bonds, the transport diffusion coefficient, and diffusion activation energy of hydronium ions and water molecules under pressure difference were evaluated and compared with relevant experimental and simulation results. Finally, the transport diffusion mechanism of hydronium ions and water molecules in PEMs was revealed and discussed. The simulation results may provide theoretical guidance for improving the performance of PEMFCs during operation.


## 2. Computational Methods

### 2.1. Simulation Details

The transport behaviors of water molecules and hydronium ions through a Nafion membrane under different temperatures and pressure gradients were studied using MD simulations. Nafion membrane is widely used as PEM material due to its reasonable proton conductivity and exceptional mechanical stability [36]. The general basic structure of a monomer of the Nafion polyelectrolyte molecule is shown in Figure 2 [26]. A Nafion 117 membrane with x=10, y=z=1, and n=10 was used in this study. Its polytetrafluoroethylene backbone is hydrophobic, and the sulfonic acid group in the side chain is hydrophilic. This structure is beneficial to phase segregation in the hydrated Nafion system [23]. In previous studies on proton mobility in Nafion membranes, four oligomers of Nafion were generally placed in a simulation box, and reasonable simulation results were obtained [24, 37, 38]. In this work, an initial amorphous structure of ten Nafion 117 was constructed in a simulation box. The steepest descent algorithm was adopted to minimize the energy of the system [39]. The convergence tolerances for energy and force were 1.0×10−4 kcal/mol and 1.0×10−6 kcal/(mol·Å), respectively. Since the structure may include unstable conformations, an annealing procedure was performed to accelerate the attainment of equilibrium [40]. The simulated annealing procedure was described in the paper of Chen et al. [24]. After the annealing procedure was done, an MD simulation using the isothermal-isobaric ensemble was conducted at 300 K for 250 ps. The density of Nafion 117 at 300 K was 1.95 g/cm3, which is in good accord with available experimental data (~2 g/cm3) [41]. In addition, similar procedures were also conducted at different temperatures (315 K, 330 K, and 350 K) to obtain membrane structures at these temperatures. The loss of proton conductivity may occur due to the membrane dehydration at temperatures higher than 373.15 K [42]. Thus, the temperature range studied in this paper was 300 K–350 K.Figure 2
Chemical structure of Nafion [26].After the molecular model of Nafion 117 was obtained, 300 hydronium ions and 3000 water molecules were generated on the left side of the PEM (Figure3). The size of the simulation box was 91.0Å×133.2Å×47.8Å. The detailed composition of the initial structure in the simulation box is listed in Table 1. As the equilibrium and steady-state properties of the system do not depend on the choice of initial conditions, all reasonable initial conditions are in principle acceptable [43]. Hydronium ions and water molecules were initially placed in arrays in the simulation box to ensure that no overlap occurred between the molecules. Two adjacent molecules did not overlap if the distance between them was greater than 3.55 Å (21/6×3.166 Å (the Lennard-Jones size of oxygen)). Moreover, all atoms were initially given the velocities in the Gaussian distribution at the set temperature (300 K, 315 K, 330 K, and 350 K). The water content (λ) was typically expressed as the number of water molecules per sulfonic acid group in membranes [44], which was ten in this study. Chen et al. found a good performance would be obtained for fuel cells when the water content was between 2.5 and 16 [24]. Based on this result, the water content in our study was moderate and reasonable. In addition, the positioned distance between the membrane and water molecules may have a significant influence on the initial stage of simulations, while it may have a relatively small influence on the simulation results in the steady state.Figure 3
Initial configuration of the simulated system. Color code: blue: fluorine atom; yellow: sulfur atom; gray, carbon atom; pink: ether bond oxygen atom; green: sulfonate oxygen atom; red: oxygen atom in water molecule; orange: oxygen atom in hydronium ion; purple: hydrogen atom in water molecule; white: hydrogen atom in hydronium ion.Table 1
Detailed composition of initial structure in the simulation box.
Number of chainsNumber of hydronium ionsNumber of water moleculesSize of the simulation boxDensity of membranes10300300091.0Å×133.2Å×47.8Å1.92 g/cm3–1.95 g/cm3The simple point charge (SPC) model was applied for water molecules [45], and the model of hydronium ions was taken from the paper of Urata et al. [45, 46]. The model of Nafion 117 was based on the Dreiding force field [47]. This force field has been used by other scholars to study the hydrated Nafion system, and reasonable results have been obtained [26, 40, 48]. Detailed parameters were provided in the Supporting Information S1 (Table SI. 1, Table SI. 2, and Table SI. 3). The cutoff distance was set to 15 Å, and the periodic boundary condition was adopted in all directions. The selection of the cutoff distance in this paper refers to the work of Mabuchi and Tokumasu [49]. Choosing 15 Å as the cutoff distance can get reasonable simulation results and save computing resources.In an actual PEMFC system, a pressure difference is generated due to the different stoichiometric ratios of hydrogen and air. As for EF-NEMD simulations, the effect of a transmembrane pressure difference can be mimicked by applying a constant force in thex direction (Figure 4) on each small molecule [31, 50]. Water molecules and hydronium ions transfer through the Nafion membrane under the driving force. The carbon atoms in the Nafion membrane were fixed in their initial positions to prevent the membrane from moving in the simulation box. The simulations were performed under different driving forces (0.15 kcal/(mol·Å), 0.25 kcal/(mol·Å), 0.35 kcal/(mol·Å), and 0.45 kcal/(mol·Å)) and temperatures (300 K, 315 K, 330 K, and 350 K) using the canonical ensemble. There is a positive correlation between driving forces and pressure gradients (i.e., a larger driving force means a higher pressure gradient) [51]. The Nosé−Hoover temperature thermostat was used in EF-NEMD simulations [52, 53]. To avoid the effect of transport velocity, the velocity in the x direction was not coupled for the temperature thermostat [54]. The particle’s motion equation was integrated using the velocity Verlet algorithm [55]. The time step of simulations was set to 1 fs. Each simulation ran 10 ns, and the simulation data in the last 6 ns were used for analysis. The model convergence was assessed by monitoring the average velocity in the x direction of small molecules within the membrane and the temperature of the system (Figure SI. 1 in the Supporting Information S2). All MD simulations were carried out using the Large-Scale Atomic-Molecular Massively Parallel Simulator (LAMMPS) program [56].Figure 4
The snapshot of a EF-NEMD simulation. Color code: blue: fluorine atom; yellow: sulfur atom; gray, carbon atom; pink: ether bond oxygen atom; green: sulfonate oxygen atom; red: oxygen atom in water molecule; orange: oxygen atom in hydronium ion; purple: hydrogen atom in water molecule; white: hydrogen atom in hydronium ion.


### 2.2. Analysis

#### 2.2.1. RDF

RDF was used to analyze the distributions of hydronium ions and water molecules in Nafion membranes. It can explain whether these molecules form clusters or channels in the systems [57]. The pair-wise correlation function describes the probability of finding an atom j at a distance r from a reference atom i, defined as [58]
(1)gr=njVNj4πr2dr,where nj is the number of j particles between r and r+dr from the atom i, V is the volume of the system, and Nj is the total number of j particles.


#### 2.2.2. Transport Diffusion Coefficient

The transport diffusion coefficient is a key parameter to characterize the transport capacity of hydronium ions and water molecules through PEMs under different pressure differences. The theoretical discussion of the transport diffusion coefficient can be found in the Supporting InformationS3 (Figure SI. 2). The transport diffusion coefficient of hydronium ions (DH3O+) and water molecules (DH2O) in EF-NEMD simulations can be calculated as follows [59]:
(2)DH3O+=RTLQn1F,DH2O=RTLQ′n2F,where T is temperature, R is gas constant, L is the thickness of PEMs, F is the external driving force, n1 and n2 are the numbers of hydronium ions and water molecules, respectively, and Q and Q’ are the molar flow rates of hydronium ions and water molecules, respectively. The detailed procedure for calculating the transport diffusion coefficient can be obtained from the previous papers [59, 60]. We performed 6 measurements to evaluate the expanded uncertainty (Δx) of the transport diffusion coefficient.
(3)Δx=kσxn′,where k is the coverage factor, σx is the standard deviation, and n′ is the number of measurements. With a level of confidence of 95%, k=2.57 for 6 measurements.


#### 2.2.3. Diffusion Activation Energy

To analyze the migration barrier of hydronium ions or water molecules through PEMs, the diffusion activation energy was obtained by the following equation [61]:
(4)DH3O+=D0exp−EH3O+RT,DH2O=D0′exp−EH2ORT,where D0 and D0′ denote the preexponential factor, and EH3O+ and EH2O are the diffusion activation energy of hydronium ions and water molecules, respectively.


#### 2.2.4. Proton Conductivity

Proton conductivity is an index to be considered when evaluating the performance of PEMFCs. It can be expressed as the following equation [62]:
(5)σ=DH3O+ne2kBT,where n is the number of hydronium ions per volume, kB is the Boltzmann constant, and e is the elementary charge.


## 2.1. Simulation Details

The transport behaviors of water molecules and hydronium ions through a Nafion membrane under different temperatures and pressure gradients were studied using MD simulations. Nafion membrane is widely used as PEM material due to its reasonable proton conductivity and exceptional mechanical stability [36]. The general basic structure of a monomer of the Nafion polyelectrolyte molecule is shown in Figure 2 [26]. A Nafion 117 membrane with x=10, y=z=1, and n=10 was used in this study. Its polytetrafluoroethylene backbone is hydrophobic, and the sulfonic acid group in the side chain is hydrophilic. This structure is beneficial to phase segregation in the hydrated Nafion system [23]. In previous studies on proton mobility in Nafion membranes, four oligomers of Nafion were generally placed in a simulation box, and reasonable simulation results were obtained [24, 37, 38]. In this work, an initial amorphous structure of ten Nafion 117 was constructed in a simulation box. The steepest descent algorithm was adopted to minimize the energy of the system [39]. The convergence tolerances for energy and force were 1.0×10−4 kcal/mol and 1.0×10−6 kcal/(mol·Å), respectively. Since the structure may include unstable conformations, an annealing procedure was performed to accelerate the attainment of equilibrium [40]. The simulated annealing procedure was described in the paper of Chen et al. [24]. After the annealing procedure was done, an MD simulation using the isothermal-isobaric ensemble was conducted at 300 K for 250 ps. The density of Nafion 117 at 300 K was 1.95 g/cm3, which is in good accord with available experimental data (~2 g/cm3) [41]. In addition, similar procedures were also conducted at different temperatures (315 K, 330 K, and 350 K) to obtain membrane structures at these temperatures. The loss of proton conductivity may occur due to the membrane dehydration at temperatures higher than 373.15 K [42]. Thus, the temperature range studied in this paper was 300 K–350 K.Figure 2
Chemical structure of Nafion [26].After the molecular model of Nafion 117 was obtained, 300 hydronium ions and 3000 water molecules were generated on the left side of the PEM (Figure3). The size of the simulation box was 91.0Å×133.2Å×47.8Å. The detailed composition of the initial structure in the simulation box is listed in Table 1. As the equilibrium and steady-state properties of the system do not depend on the choice of initial conditions, all reasonable initial conditions are in principle acceptable [43]. Hydronium ions and water molecules were initially placed in arrays in the simulation box to ensure that no overlap occurred between the molecules. Two adjacent molecules did not overlap if the distance between them was greater than 3.55 Å (21/6×3.166 Å (the Lennard-Jones size of oxygen)). Moreover, all atoms were initially given the velocities in the Gaussian distribution at the set temperature (300 K, 315 K, 330 K, and 350 K). The water content (λ) was typically expressed as the number of water molecules per sulfonic acid group in membranes [44], which was ten in this study. Chen et al. found a good performance would be obtained for fuel cells when the water content was between 2.5 and 16 [24]. Based on this result, the water content in our study was moderate and reasonable. In addition, the positioned distance between the membrane and water molecules may have a significant influence on the initial stage of simulations, while it may have a relatively small influence on the simulation results in the steady state.Figure 3
Initial configuration of the simulated system. Color code: blue: fluorine atom; yellow: sulfur atom; gray, carbon atom; pink: ether bond oxygen atom; green: sulfonate oxygen atom; red: oxygen atom in water molecule; orange: oxygen atom in hydronium ion; purple: hydrogen atom in water molecule; white: hydrogen atom in hydronium ion.Table 1
Detailed composition of initial structure in the simulation box.
Number of chainsNumber of hydronium ionsNumber of water moleculesSize of the simulation boxDensity of membranes10300300091.0Å×133.2Å×47.8Å1.92 g/cm3–1.95 g/cm3The simple point charge (SPC) model was applied for water molecules [45], and the model of hydronium ions was taken from the paper of Urata et al. [45, 46]. The model of Nafion 117 was based on the Dreiding force field [47]. This force field has been used by other scholars to study the hydrated Nafion system, and reasonable results have been obtained [26, 40, 48]. Detailed parameters were provided in the Supporting Information S1 (Table SI. 1, Table SI. 2, and Table SI. 3). The cutoff distance was set to 15 Å, and the periodic boundary condition was adopted in all directions. The selection of the cutoff distance in this paper refers to the work of Mabuchi and Tokumasu [49]. Choosing 15 Å as the cutoff distance can get reasonable simulation results and save computing resources.In an actual PEMFC system, a pressure difference is generated due to the different stoichiometric ratios of hydrogen and air. As for EF-NEMD simulations, the effect of a transmembrane pressure difference can be mimicked by applying a constant force in thex direction (Figure 4) on each small molecule [31, 50]. Water molecules and hydronium ions transfer through the Nafion membrane under the driving force. The carbon atoms in the Nafion membrane were fixed in their initial positions to prevent the membrane from moving in the simulation box. The simulations were performed under different driving forces (0.15 kcal/(mol·Å), 0.25 kcal/(mol·Å), 0.35 kcal/(mol·Å), and 0.45 kcal/(mol·Å)) and temperatures (300 K, 315 K, 330 K, and 350 K) using the canonical ensemble. There is a positive correlation between driving forces and pressure gradients (i.e., a larger driving force means a higher pressure gradient) [51]. The Nosé−Hoover temperature thermostat was used in EF-NEMD simulations [52, 53]. To avoid the effect of transport velocity, the velocity in the x direction was not coupled for the temperature thermostat [54]. The particle’s motion equation was integrated using the velocity Verlet algorithm [55]. The time step of simulations was set to 1 fs. Each simulation ran 10 ns, and the simulation data in the last 6 ns were used for analysis. The model convergence was assessed by monitoring the average velocity in the x direction of small molecules within the membrane and the temperature of the system (Figure SI. 1 in the Supporting Information S2). All MD simulations were carried out using the Large-Scale Atomic-Molecular Massively Parallel Simulator (LAMMPS) program [56].Figure 4
The snapshot of a EF-NEMD simulation. Color code: blue: fluorine atom; yellow: sulfur atom; gray, carbon atom; pink: ether bond oxygen atom; green: sulfonate oxygen atom; red: oxygen atom in water molecule; orange: oxygen atom in hydronium ion; purple: hydrogen atom in water molecule; white: hydrogen atom in hydronium ion.


## 2.2. Analysis

### 2.2.1. RDF

RDF was used to analyze the distributions of hydronium ions and water molecules in Nafion membranes. It can explain whether these molecules form clusters or channels in the systems [57]. The pair-wise correlation function describes the probability of finding an atom j at a distance r from a reference atom i, defined as [58]
(1)gr=njVNj4πr2dr,where nj is the number of j particles between r and r+dr from the atom i, V is the volume of the system, and Nj is the total number of j particles.


### 2.2.2. Transport Diffusion Coefficient

The transport diffusion coefficient is a key parameter to characterize the transport capacity of hydronium ions and water molecules through PEMs under different pressure differences. The theoretical discussion of the transport diffusion coefficient can be found in the Supporting InformationS3 (Figure SI. 2). The transport diffusion coefficient of hydronium ions (DH3O+) and water molecules (DH2O) in EF-NEMD simulations can be calculated as follows [59]:
(2)DH3O+=RTLQn1F,DH2O=RTLQ′n2F,where T is temperature, R is gas constant, L is the thickness of PEMs, F is the external driving force, n1 and n2 are the numbers of hydronium ions and water molecules, respectively, and Q and Q’ are the molar flow rates of hydronium ions and water molecules, respectively. The detailed procedure for calculating the transport diffusion coefficient can be obtained from the previous papers [59, 60]. We performed 6 measurements to evaluate the expanded uncertainty (Δx) of the transport diffusion coefficient.
(3)Δx=kσxn′,where k is the coverage factor, σx is the standard deviation, and n′ is the number of measurements. With a level of confidence of 95%, k=2.57 for 6 measurements.


### 2.2.3. Diffusion Activation Energy

To analyze the migration barrier of hydronium ions or water molecules through PEMs, the diffusion activation energy was obtained by the following equation [61]:
(4)DH3O+=D0exp−EH3O+RT,DH2O=D0′exp−EH2ORT,where D0 and D0′ denote the preexponential factor, and EH3O+ and EH2O are the diffusion activation energy of hydronium ions and water molecules, respectively.


### 2.2.4. Proton Conductivity

Proton conductivity is an index to be considered when evaluating the performance of PEMFCs. It can be expressed as the following equation [62]:
(5)σ=DH3O+ne2kBT,where n is the number of hydronium ions per volume, kB is the Boltzmann constant, and e is the elementary charge.


## 2.2.1. RDF

RDF was used to analyze the distributions of hydronium ions and water molecules in Nafion membranes. It can explain whether these molecules form clusters or channels in the systems [57]. The pair-wise correlation function describes the probability of finding an atom j at a distance r from a reference atom i, defined as [58]
(1)gr=njVNj4πr2dr,where nj is the number of j particles between r and r+dr from the atom i, V is the volume of the system, and Nj is the total number of j particles.


## 2.2.2. Transport Diffusion Coefficient

The transport diffusion coefficient is a key parameter to characterize the transport capacity of hydronium ions and water molecules through PEMs under different pressure differences. The theoretical discussion of the transport diffusion coefficient can be found in the Supporting InformationS3 (Figure SI. 2). The transport diffusion coefficient of hydronium ions (DH3O+) and water molecules (DH2O) in EF-NEMD simulations can be calculated as follows [59]:
(2)DH3O+=RTLQn1F,DH2O=RTLQ′n2F,where T is temperature, R is gas constant, L is the thickness of PEMs, F is the external driving force, n1 and n2 are the numbers of hydronium ions and water molecules, respectively, and Q and Q’ are the molar flow rates of hydronium ions and water molecules, respectively. The detailed procedure for calculating the transport diffusion coefficient can be obtained from the previous papers [59, 60]. We performed 6 measurements to evaluate the expanded uncertainty (Δx) of the transport diffusion coefficient.
(3)Δx=kσxn′,where k is the coverage factor, σx is the standard deviation, and n′ is the number of measurements. With a level of confidence of 95%, k=2.57 for 6 measurements.


## 2.2.3. Diffusion Activation Energy

To analyze the migration barrier of hydronium ions or water molecules through PEMs, the diffusion activation energy was obtained by the following equation [61]:
(4)DH3O+=D0exp−EH3O+RT,DH2O=D0′exp−EH2ORT,where D0 and D0′ denote the preexponential factor, and EH3O+ and EH2O are the diffusion activation energy of hydronium ions and water molecules, respectively.


## 2.2.4. Proton Conductivity

Proton conductivity is an index to be considered when evaluating the performance of PEMFCs. It can be expressed as the following equation [62]:
(5)σ=DH3O+ne2kBT,where n is the number of hydronium ions per volume, kB is the Boltzmann constant, and e is the elementary charge.


## 3. Results and Discussion

### 3.1. Distribution of Hydronium Ions and Water Molecules in PEMs

#### 3.1.1. The RDF between Hydronium Ion/Water Molecule and Main Chain

Proton transport is closely related to the distribution and transport of water molecules in hydrated PEMs [26]. To better understand the transport diffusion characteristics of hydronium ions and water molecules in Nafion membranes, the RDFs between the carbon atom in the hydrophobic fluorocarbon main chain and the hydrogen atom in the water molecule (gC-Hw(r)) and between the carbon atom in the main chain and the hydrogen atom in hydronium ion (gC-Hh(r)) are calculated and shown in Figure 5. Both gC−Hwr and gC−Hhr increase significantly as the distance between the carbon atom and the hydrogen atom in hydronium ion or water molecule increases from 2 Å to 4.8 Å, then vary in a small range as the distance increases further. It indicates that the maximum probability of a hydrogen atom in a hydronium ion or water molecule being around a carbon atom in the main chain occurs at a distance of 4.8 Å. These broad peaks of gC−Hhr and gC−Hwr indicate that the hydrophobic main chain has a weak effect on the motion of hydronium ions and water molecules [23]. It is worth noting that gC−Hhr and gC−Hwr do not approach 1 at r of 15 Å. Water molecules and hydronium ions tend to adsorb in PEMs due to the van der Waals and Coulomb forces, resulting in a higher local density of small molecules in membranes than their average density in the system [63]. Additionally, the value of gC−Hwr is greater than that of gC−Hhr, which indicates that the water molecule is more likely to appear near a carbon atom in the fluorocarbon backbone than the hydronium ion. It is because the hydronium ions (H3O+) with positive charges tend to gather around sulfonic acid groups (SO3-) with negative charges and stay away from the main chain [64].Figure 5
The RDFs of carbon atom in fluorocarbon backbone and hydrogen atom in water molecule (gC−Hwr) and hydronium ion (gC−Hhr) at 300 K and constant driving force (0.15 kcal/(mol·Å)).


#### 3.1.2. The RDF between Water Molecule/Hydronium Ion and Side Chain

The proton conductivity in PEMs is strongly dependent upon the hydrogen bonding network [65]. The RDF can be used to describe the strength of hydrogen bonds [66]. The RDFs of the oxygen atom in sulfonic acid group (Os) and the hydrogen atom in water molecule (Hw) or hydronium ion (Hh) are given in Figure 6, respectively. It shows that gOs−Hhr exhibits two significant peaks at 2 Å and 3.4 Å, and the third peak is formed near 5.2 Å. The significant peaks of gOs−Hhr correspond to the contact ion pair, and the third peak corresponds to the solvent separated ion pair [23]. Moreover, the dominant peaks of gOs−Hhr and gOs−Hwr appear at roughly the same distance, which indicates a relation between hydronium ions and water molecules. Aggregates or clusters of positively charged hydronium ions and neutral water molecules are formed around negatively charged hydrophilic sulfonic acid groups [57]. These hydrophilic aggregates are embedded in the hydrophobic phase of the fluorocarbon backbone, and they induce the conducting channels in PEMs [67]. In addition, the peak value of gOs−Hhr is much larger than that of gOs−Hwr. Since hydrogen bonds exist between sulfonic acid groups and water molecules/hydronium ions if r is less than 2.7 Å (the position of the first minimum) for gO−Hr [66], hydronium ions form stronger hydrogen bonds with sulfonic acid groups than water molecules. This is because a hydronium ion with a positive charge is more easily attracted by the oxygen atoms with negative charges in sulfonic acid groups than a water molecule [64].Figure 6
The RDFs ofOs and Hw/Hh at 300 K and constant driving force (0.15 kcal/(mol·Å)).


#### 3.1.3. The RDF between Water Molecule and Hydronium Ion

Hydrogen-bonded linkages made of hydronium ions and water molecules could affect the transport of protons [23]. Thus, gOw−Hhr and gOh−Hwr are calculated in this study and shown in Figure 7. The gOw−Hhr is RDF between the oxygen atom in water molecule (Ow) and the hydrogen atom in hydronium ion (Hh), and the gOh−Hwr denotes RDF between the oxygen atom in hydronium ion (Oh) and the hydrogen atom in water molecule (Hw). The values of gOw−Hhr and gOh−Hwr are 0 when r is smaller than 1.6 Å and 2.1 Å, respectively, which means that the least distances are 1.6 Å for gOw−Hhr and 2.1 Å for gOh−Hwr due to strong repulsive force among the atoms at these distances. In addition, the values of these two RDFs approach 1 as r increases to 15 Å, which indicates that the property of the system becomes relatively homogenous when r gets beyond 15 Å. There are three distinct peaks both for gOw−Hhr and gOh−Hwr, and they correspond to different solvation shells [23]. It indicates that water molecules tend to aggregate with hydronium ions to form hydrophilic clusters [5]. The dominant peaks of gOw−Hhr and gOh−Hhr are at about 2.1 Å and 3.7 Å, indicating that hydrogen bonds form more easily between Ow and Hh.Figure 7
The RDFs between water molecule and hydronium ion at 300 K and a constant driving force (0.15 kcal/(mol·Å)). The hydrogen and oxygen in water molecule are denoted by “Hw” and “Ow”, respectively. The hydrogen and oxygen in hydronium ion are denoted by “Hh” and “Oh”, respectively.


### 3.2. Transport Diffusion of Hydronium Ions and Water Molecules

#### 3.2.1. Effect of Pressure Difference on Transport Diffusion Characteristics

To study the effect of pressure difference on the transport diffusion characteristics of hydronium ions and water molecules, different driving forces (0.15 kcal/(mol·Å)–0.45 kcal/(mol·Å)) were applied to hydronium ions and water molecules at 300 K. Water molecules and hydronium ions transported from anode to cathode through the Nafion 117 membrane under these driving forces. The transport diffusion coefficients of hydronium ions and water molecules are calculated based on Equation (2) and are given in Figure 8. As the driving force rises, the transport diffusion coefficient of water molecules increases from 0.88×10−9 m2/s to 3.36×10−9 m2/s, and that of hydronium ions also increases monotonically (0.15×10−9m2/s–0.46×10−9 m2/s). This trend is in line with that in the experiment of Husar et al. [30]. They studied water transfer through PEMs under pressure differences and found that the water mass-transfer flux increased with increasing pressure drop [30]. The higher pressure difference makes it easier for hydronium ions and water molecules to overcome the energy barrier and transport through the Nafion membrane. Additionally, the transport diffusion coefficient of water molecules exceeds that of hydronium ions under driving forces, which is similar to the simulation results displayed by Chen et al. [24]. The lighter water molecules have higher velocity and better transport diffusion capacity than hydronium ions. Moreover, a hydronium ion has an entirely positive charge, while a water molecule is neutral. The movement of hydronium ions is more susceptible to the negatively charged sulfonic acid group. The potential energy between oxygen atom in sulfonic acid group (Os) and hydronium ion (-28.37 kJ/mol) is lower than that between Os and water molecule (-13.17 kJ/mol) based on our statistics, which supports the above point. As a result, the migration capacity of hydronium ions in PEMs is lower than that of water molecules.Figure 8
The transport diffusion coefficients of hydronium ions and water molecules under different driving forces.Proton conductivity is a significant index to evaluate PEM performance. According to Equation (4), the proton conductivity is calculated and given in Figure 9. The proton conductivity increases from 0.005 S/cm to 0.015 S/cm as the driving force rises from 0.15 kcal/(mol·Å) to 0.45 kcal/(mol·Å). The proton conductivity in the simulations is lower than that obtained by the empirical Equation (6) (0.046 S/cm) [68]. This is because the water content of the PEM is lower than that of the system in simulations. The proton conductivity and transport diffusion coefficient of hydronium ions have the same trend with driving force. It indicates that the performance of PEMFCs could be improved by keeping the pressure difference between the anode and cathode in a reasonable range. For an actual PEMFC system, the pressure difference is generally less than 30 kPa [17].
(6)σ=5.139λ−3.261000exp1268303−1268T.Figure 9
The proton conductivity under different driving forces.


#### 3.2.2. Effect of Temperature on Transport Diffusion Characteristics

The hydrogen bond dynamics and proton transfer depend on temperature [24, 69]. The EF-NEMD simulations of water molecule and hydronium ion transport through Nafion membranes were performed under a constant driving force (0.35 kcal/(mol·Å)) and different temperatures (300 K–350 K), and the simulation results are given in Figure 10(a). It shows that the transport diffusion coefficients of hydronium ions (0.31×10−9m2/s–0.95×10−9 m2/s) and water molecules (2.25×10−9m2/s–6.89×10−9 m2/s) both increase as temperature rises. In addition, the transport diffusion coefficient of water molecules is significantly higher than that of hydronium ions.Figure 10
(a) The transport diffusion coefficients of hydronium ions and water molecules under different temperatures. (b) The Arrhenius plots of ln (D) against 1000/T.
(a)(b)At an elevated temperature, the energy of the molecules increases, making it easier for water and protons to cross the energy barrier and transport through PEMs [24]. High temperature also changes the pore structure and thus affects the diffusivity of water molecules and hydronium ions within the proton exchange membrane. The PoreBlazer v4.0 tool was utilized to accurately measure the pore characteristics of Nafion membranes at different temperatures [70, 71]. The pore size distribution, and the pore volume and surface area of Nafion membranes (Figure 11) were calculated by a helium probe. It can be seen that the pore size distribution in Nafion membranes exhibits a bimodal shape (Figure 11(a)). The two peaks shift to the right with increasing temperature, indicating that the pore size increases at high temperatures. As the temperature rises, the volume and surface area of pores also increase (Figure 11(b)). The enlarged pore volume at high temperatures also lowers the energy barrier for water molecules and hydronium ions across PEMs, thereby increasing their transport diffusion coefficients. Moreover, the Arrhenius plots of transport diffusion coefficients for hydronium ions and water molecules are presented in Figure 10(b). According to Equation (3), the diffusion activation energies of hydronium ions and water molecules are 21.15 kJ/mol and 20.25 kJ/mol, respectively. The lower molecular weight and diffusion activation energy of water molecules result in higher transport diffusion ability of water molecules.Figure 11
(a) Pore size distributions of Nafion membranes at different temperatures.Vp is the pore volume, and dp is the pore size. (b) The volume and surface area of pores in Nafion membranes as a function of temperature.
(a)(b)The effect of temperature on proton conductivity is shown in Figure12(a). The proton conductivity is positively correlated with temperature, which is consistent with the correlation between the transport diffusion coefficient of hydronium ions and temperature. To quantify the number of hydrogen bonds in the simulated systems at different temperatures, geometric criteria for determining hydrogen bond formation were used. Two oxygen atoms are thought to form a hydrogen bond via a hydrogen atom when the following two criteria are met simultaneously: (I) the distance between the oxygen atoms is smaller than 3.5 Å and (II) the OHO angle is greater than 120° [67]. The number of hydrogen bonds (hydronium ion acts as donor and water molecule/sulfonic acid group acts as acceptors) at different temperatures is displayed in Figure 12(b). As temperature rises from 300 K to 350 K, the number of hydrogen bonds between hydronium ions and water molecules and that between hydronium ions and sulfonic acid group decreases. Elevated temperature allows the hydronium ions to break existing hydrogen bonds between hydronium ions and water molecules/sulfonic acid groups and readily form new hydrogen bonds with other molecules, which accelerates hydrogen bond dynamics and promotes proton transfer [69, 72]. Moreover, hydronium ions are more active at high temperatures, which increases the probability of them entering the free region (the middle region of the simulation box in Figure 3) through the PEM. Interactions between molecules in the free region are weakened at high temperatures, resulting in a reduced number of hydrogen bonds. High temperature accelerates the passage of protons through PEMs, indicating that it is beneficial to improve the performance of PEMs at a specific moisture content.Figure 12
(a) The proton conductivity and (b) the number of hydrogen bonds (hydronium ion acts as donor, and water molecule/sulfonic acid group acts as acceptors) at different temperatures.
(a)(b)


#### 3.2.3. Comparison with Existing Works

In most molecular dynamics simulations, the self-diffusion coefficients of hydronium ions and water molecules in PEMs were calculated from the Einstein equation based on the mean square displacement [2, 7, 19, 20]. Self-diffusion was used to characterize the diffusion capacity of molecules at equilibrium (i.e., without driving forces such as the gradients in pressure and chemical potential). However, this paper focuses on the transport diffusion of hydronium ions and water molecules in PEMs under pressure differences, and the transport diffusion coefficients of molecules were calculated and discussed. It is found that the transport diffusion coefficients of water (2.252×10−9 m2/s) and hydronium ion (0.309×10−9 m2/s) at 300 K in this work are higher than the self-diffusion coefficients of water (0.591×10−9 m2/s) and hydronium ion (0.158×10−9 m2/s) in the simulation of Jang et al. under the same condition [26]. Moreover, the transport diffusion coefficient of water (2.252×10−9 m2/s) in this study is about an order of magnitude higher than the Fickian diffusion coefficient of water (~1×10−10 m2/s) calculated by Li et al. [34]. In addition, the diffusion activation energy of water (Table 2) in this study (20.25 kJ/mol) is consistent with the experimental results of Yeo et al. (20.2 kJ/mol) and Takamatsu et al. (20.51 kJ/mol) [73, 74]. The diffusion activation energy of hydronium ion (Table 3) in this study (21.15 kJ/mol) is close to the experimental results of Thompson et al. (21.23 kJ/mol) and Cappadonia et al. (21.22 kJ/mol) [75, 76]. The diffusion activation energy of hydronium ion in the experiment of Halim et al. (13.5 kJ/mol) is relatively low [77], which is due to the high water content (λ=23) of the membranes in their experiment. The diffusion activation energy generally decreases with increasing water content [14]. When the water content in PEM is high, the diffusion activation energy of hydronium ion is close to that of hydronium ion in bulk water (8.37–12.56 kJ/mol) [6].Table 2
The diffusion activation energy of water in the Nafion membrane.
EH2O (kJ/Mol)ConditionReference~23298.15K≤T≤373.15KThe experiment of Morris and Sun [78]20.51293.15K≤T≤354.15KThe experiment of Takamatsu et al. [74]20.20K≤T≤382.15KThe experiment of Yeo and Eisenberg [73]20.25300K<T<350KThis workTable 3
The diffusion activation energy of hydronium ion in Nafion membrane.
EH3O+ (kJ/Mol)ConditionReference22.6Dried sample,295.15K≤T≤335.15KThe experiment of Ye et al. [79]13.5λ=23,293.15K≤T≤363.15KThe experiment of Halim et al. [77]21.23λ=3.7,253.15K<T<353.15KThe experiment of Thompson et al. [75]21.224.2 g H2O/100 g Nafion, 225 K < T <300 KThe experiment of Cappadonia et al. [76]21.15λ=10,300K<T<350KThis work


## 3.1. Distribution of Hydronium Ions and Water Molecules in PEMs

### 3.1.1. The RDF between Hydronium Ion/Water Molecule and Main Chain

Proton transport is closely related to the distribution and transport of water molecules in hydrated PEMs [26]. To better understand the transport diffusion characteristics of hydronium ions and water molecules in Nafion membranes, the RDFs between the carbon atom in the hydrophobic fluorocarbon main chain and the hydrogen atom in the water molecule (gC-Hw(r)) and between the carbon atom in the main chain and the hydrogen atom in hydronium ion (gC-Hh(r)) are calculated and shown in Figure 5. Both gC−Hwr and gC−Hhr increase significantly as the distance between the carbon atom and the hydrogen atom in hydronium ion or water molecule increases from 2 Å to 4.8 Å, then vary in a small range as the distance increases further. It indicates that the maximum probability of a hydrogen atom in a hydronium ion or water molecule being around a carbon atom in the main chain occurs at a distance of 4.8 Å. These broad peaks of gC−Hhr and gC−Hwr indicate that the hydrophobic main chain has a weak effect on the motion of hydronium ions and water molecules [23]. It is worth noting that gC−Hhr and gC−Hwr do not approach 1 at r of 15 Å. Water molecules and hydronium ions tend to adsorb in PEMs due to the van der Waals and Coulomb forces, resulting in a higher local density of small molecules in membranes than their average density in the system [63]. Additionally, the value of gC−Hwr is greater than that of gC−Hhr, which indicates that the water molecule is more likely to appear near a carbon atom in the fluorocarbon backbone than the hydronium ion. It is because the hydronium ions (H3O+) with positive charges tend to gather around sulfonic acid groups (SO3-) with negative charges and stay away from the main chain [64].Figure 5
The RDFs of carbon atom in fluorocarbon backbone and hydrogen atom in water molecule (gC−Hwr) and hydronium ion (gC−Hhr) at 300 K and constant driving force (0.15 kcal/(mol·Å)).


### 3.1.2. The RDF between Water Molecule/Hydronium Ion and Side Chain

The proton conductivity in PEMs is strongly dependent upon the hydrogen bonding network [65]. The RDF can be used to describe the strength of hydrogen bonds [66]. The RDFs of the oxygen atom in sulfonic acid group (Os) and the hydrogen atom in water molecule (Hw) or hydronium ion (Hh) are given in Figure 6, respectively. It shows that gOs−Hhr exhibits two significant peaks at 2 Å and 3.4 Å, and the third peak is formed near 5.2 Å. The significant peaks of gOs−Hhr correspond to the contact ion pair, and the third peak corresponds to the solvent separated ion pair [23]. Moreover, the dominant peaks of gOs−Hhr and gOs−Hwr appear at roughly the same distance, which indicates a relation between hydronium ions and water molecules. Aggregates or clusters of positively charged hydronium ions and neutral water molecules are formed around negatively charged hydrophilic sulfonic acid groups [57]. These hydrophilic aggregates are embedded in the hydrophobic phase of the fluorocarbon backbone, and they induce the conducting channels in PEMs [67]. In addition, the peak value of gOs−Hhr is much larger than that of gOs−Hwr. Since hydrogen bonds exist between sulfonic acid groups and water molecules/hydronium ions if r is less than 2.7 Å (the position of the first minimum) for gO−Hr [66], hydronium ions form stronger hydrogen bonds with sulfonic acid groups than water molecules. This is because a hydronium ion with a positive charge is more easily attracted by the oxygen atoms with negative charges in sulfonic acid groups than a water molecule [64].Figure 6
The RDFs ofOs and Hw/Hh at 300 K and constant driving force (0.15 kcal/(mol·Å)).


### 3.1.3. The RDF between Water Molecule and Hydronium Ion

Hydrogen-bonded linkages made of hydronium ions and water molecules could affect the transport of protons [23]. Thus, gOw−Hhr and gOh−Hwr are calculated in this study and shown in Figure 7. The gOw−Hhr is RDF between the oxygen atom in water molecule (Ow) and the hydrogen atom in hydronium ion (Hh), and the gOh−Hwr denotes RDF between the oxygen atom in hydronium ion (Oh) and the hydrogen atom in water molecule (Hw). The values of gOw−Hhr and gOh−Hwr are 0 when r is smaller than 1.6 Å and 2.1 Å, respectively, which means that the least distances are 1.6 Å for gOw−Hhr and 2.1 Å for gOh−Hwr due to strong repulsive force among the atoms at these distances. In addition, the values of these two RDFs approach 1 as r increases to 15 Å, which indicates that the property of the system becomes relatively homogenous when r gets beyond 15 Å. There are three distinct peaks both for gOw−Hhr and gOh−Hwr, and they correspond to different solvation shells [23]. It indicates that water molecules tend to aggregate with hydronium ions to form hydrophilic clusters [5]. The dominant peaks of gOw−Hhr and gOh−Hhr are at about 2.1 Å and 3.7 Å, indicating that hydrogen bonds form more easily between Ow and Hh.Figure 7
The RDFs between water molecule and hydronium ion at 300 K and a constant driving force (0.15 kcal/(mol·Å)). The hydrogen and oxygen in water molecule are denoted by “Hw” and “Ow”, respectively. The hydrogen and oxygen in hydronium ion are denoted by “Hh” and “Oh”, respectively.


## 3.1.1. The RDF between Hydronium Ion/Water Molecule and Main Chain

Proton transport is closely related to the distribution and transport of water molecules in hydrated PEMs [26]. To better understand the transport diffusion characteristics of hydronium ions and water molecules in Nafion membranes, the RDFs between the carbon atom in the hydrophobic fluorocarbon main chain and the hydrogen atom in the water molecule (gC-Hw(r)) and between the carbon atom in the main chain and the hydrogen atom in hydronium ion (gC-Hh(r)) are calculated and shown in Figure 5. Both gC−Hwr and gC−Hhr increase significantly as the distance between the carbon atom and the hydrogen atom in hydronium ion or water molecule increases from 2 Å to 4.8 Å, then vary in a small range as the distance increases further. It indicates that the maximum probability of a hydrogen atom in a hydronium ion or water molecule being around a carbon atom in the main chain occurs at a distance of 4.8 Å. These broad peaks of gC−Hhr and gC−Hwr indicate that the hydrophobic main chain has a weak effect on the motion of hydronium ions and water molecules [23]. It is worth noting that gC−Hhr and gC−Hwr do not approach 1 at r of 15 Å. Water molecules and hydronium ions tend to adsorb in PEMs due to the van der Waals and Coulomb forces, resulting in a higher local density of small molecules in membranes than their average density in the system [63]. Additionally, the value of gC−Hwr is greater than that of gC−Hhr, which indicates that the water molecule is more likely to appear near a carbon atom in the fluorocarbon backbone than the hydronium ion. It is because the hydronium ions (H3O+) with positive charges tend to gather around sulfonic acid groups (SO3-) with negative charges and stay away from the main chain [64].Figure 5
The RDFs of carbon atom in fluorocarbon backbone and hydrogen atom in water molecule (gC−Hwr) and hydronium ion (gC−Hhr) at 300 K and constant driving force (0.15 kcal/(mol·Å)).


## 3.1.2. The RDF between Water Molecule/Hydronium Ion and Side Chain

The proton conductivity in PEMs is strongly dependent upon the hydrogen bonding network [65]. The RDF can be used to describe the strength of hydrogen bonds [66]. The RDFs of the oxygen atom in sulfonic acid group (Os) and the hydrogen atom in water molecule (Hw) or hydronium ion (Hh) are given in Figure 6, respectively. It shows that gOs−Hhr exhibits two significant peaks at 2 Å and 3.4 Å, and the third peak is formed near 5.2 Å. The significant peaks of gOs−Hhr correspond to the contact ion pair, and the third peak corresponds to the solvent separated ion pair [23]. Moreover, the dominant peaks of gOs−Hhr and gOs−Hwr appear at roughly the same distance, which indicates a relation between hydronium ions and water molecules. Aggregates or clusters of positively charged hydronium ions and neutral water molecules are formed around negatively charged hydrophilic sulfonic acid groups [57]. These hydrophilic aggregates are embedded in the hydrophobic phase of the fluorocarbon backbone, and they induce the conducting channels in PEMs [67]. In addition, the peak value of gOs−Hhr is much larger than that of gOs−Hwr. Since hydrogen bonds exist between sulfonic acid groups and water molecules/hydronium ions if r is less than 2.7 Å (the position of the first minimum) for gO−Hr [66], hydronium ions form stronger hydrogen bonds with sulfonic acid groups than water molecules. This is because a hydronium ion with a positive charge is more easily attracted by the oxygen atoms with negative charges in sulfonic acid groups than a water molecule [64].Figure 6
The RDFs ofOs and Hw/Hh at 300 K and constant driving force (0.15 kcal/(mol·Å)).


## 3.1.3. The RDF between Water Molecule and Hydronium Ion

Hydrogen-bonded linkages made of hydronium ions and water molecules could affect the transport of protons [23]. Thus, gOw−Hhr and gOh−Hwr are calculated in this study and shown in Figure 7. The gOw−Hhr is RDF between the oxygen atom in water molecule (Ow) and the hydrogen atom in hydronium ion (Hh), and the gOh−Hwr denotes RDF between the oxygen atom in hydronium ion (Oh) and the hydrogen atom in water molecule (Hw). The values of gOw−Hhr and gOh−Hwr are 0 when r is smaller than 1.6 Å and 2.1 Å, respectively, which means that the least distances are 1.6 Å for gOw−Hhr and 2.1 Å for gOh−Hwr due to strong repulsive force among the atoms at these distances. In addition, the values of these two RDFs approach 1 as r increases to 15 Å, which indicates that the property of the system becomes relatively homogenous when r gets beyond 15 Å. There are three distinct peaks both for gOw−Hhr and gOh−Hwr, and they correspond to different solvation shells [23]. It indicates that water molecules tend to aggregate with hydronium ions to form hydrophilic clusters [5]. The dominant peaks of gOw−Hhr and gOh−Hhr are at about 2.1 Å and 3.7 Å, indicating that hydrogen bonds form more easily between Ow and Hh.Figure 7
The RDFs between water molecule and hydronium ion at 300 K and a constant driving force (0.15 kcal/(mol·Å)). The hydrogen and oxygen in water molecule are denoted by “Hw” and “Ow”, respectively. The hydrogen and oxygen in hydronium ion are denoted by “Hh” and “Oh”, respectively.


## 3.2. Transport Diffusion of Hydronium Ions and Water Molecules

### 3.2.1. Effect of Pressure Difference on Transport Diffusion Characteristics

To study the effect of pressure difference on the transport diffusion characteristics of hydronium ions and water molecules, different driving forces (0.15 kcal/(mol·Å)–0.45 kcal/(mol·Å)) were applied to hydronium ions and water molecules at 300 K. Water molecules and hydronium ions transported from anode to cathode through the Nafion 117 membrane under these driving forces. The transport diffusion coefficients of hydronium ions and water molecules are calculated based on Equation (2) and are given in Figure 8. As the driving force rises, the transport diffusion coefficient of water molecules increases from 0.88×10−9 m2/s to 3.36×10−9 m2/s, and that of hydronium ions also increases monotonically (0.15×10−9m2/s–0.46×10−9 m2/s). This trend is in line with that in the experiment of Husar et al. [30]. They studied water transfer through PEMs under pressure differences and found that the water mass-transfer flux increased with increasing pressure drop [30]. The higher pressure difference makes it easier for hydronium ions and water molecules to overcome the energy barrier and transport through the Nafion membrane. Additionally, the transport diffusion coefficient of water molecules exceeds that of hydronium ions under driving forces, which is similar to the simulation results displayed by Chen et al. [24]. The lighter water molecules have higher velocity and better transport diffusion capacity than hydronium ions. Moreover, a hydronium ion has an entirely positive charge, while a water molecule is neutral. The movement of hydronium ions is more susceptible to the negatively charged sulfonic acid group. The potential energy between oxygen atom in sulfonic acid group (Os) and hydronium ion (-28.37 kJ/mol) is lower than that between Os and water molecule (-13.17 kJ/mol) based on our statistics, which supports the above point. As a result, the migration capacity of hydronium ions in PEMs is lower than that of water molecules.Figure 8
The transport diffusion coefficients of hydronium ions and water molecules under different driving forces.Proton conductivity is a significant index to evaluate PEM performance. According to Equation (4), the proton conductivity is calculated and given in Figure 9. The proton conductivity increases from 0.005 S/cm to 0.015 S/cm as the driving force rises from 0.15 kcal/(mol·Å) to 0.45 kcal/(mol·Å). The proton conductivity in the simulations is lower than that obtained by the empirical Equation (6) (0.046 S/cm) [68]. This is because the water content of the PEM is lower than that of the system in simulations. The proton conductivity and transport diffusion coefficient of hydronium ions have the same trend with driving force. It indicates that the performance of PEMFCs could be improved by keeping the pressure difference between the anode and cathode in a reasonable range. For an actual PEMFC system, the pressure difference is generally less than 30 kPa [17].
(6)σ=5.139λ−3.261000exp1268303−1268T.Figure 9
The proton conductivity under different driving forces.


### 3.2.2. Effect of Temperature on Transport Diffusion Characteristics

The hydrogen bond dynamics and proton transfer depend on temperature [24, 69]. The EF-NEMD simulations of water molecule and hydronium ion transport through Nafion membranes were performed under a constant driving force (0.35 kcal/(mol·Å)) and different temperatures (300 K–350 K), and the simulation results are given in Figure 10(a). It shows that the transport diffusion coefficients of hydronium ions (0.31×10−9m2/s–0.95×10−9 m2/s) and water molecules (2.25×10−9m2/s–6.89×10−9 m2/s) both increase as temperature rises. In addition, the transport diffusion coefficient of water molecules is significantly higher than that of hydronium ions.Figure 10
(a) The transport diffusion coefficients of hydronium ions and water molecules under different temperatures. (b) The Arrhenius plots of ln (D) against 1000/T.
(a)(b)At an elevated temperature, the energy of the molecules increases, making it easier for water and protons to cross the energy barrier and transport through PEMs [24]. High temperature also changes the pore structure and thus affects the diffusivity of water molecules and hydronium ions within the proton exchange membrane. The PoreBlazer v4.0 tool was utilized to accurately measure the pore characteristics of Nafion membranes at different temperatures [70, 71]. The pore size distribution, and the pore volume and surface area of Nafion membranes (Figure 11) were calculated by a helium probe. It can be seen that the pore size distribution in Nafion membranes exhibits a bimodal shape (Figure 11(a)). The two peaks shift to the right with increasing temperature, indicating that the pore size increases at high temperatures. As the temperature rises, the volume and surface area of pores also increase (Figure 11(b)). The enlarged pore volume at high temperatures also lowers the energy barrier for water molecules and hydronium ions across PEMs, thereby increasing their transport diffusion coefficients. Moreover, the Arrhenius plots of transport diffusion coefficients for hydronium ions and water molecules are presented in Figure 10(b). According to Equation (3), the diffusion activation energies of hydronium ions and water molecules are 21.15 kJ/mol and 20.25 kJ/mol, respectively. The lower molecular weight and diffusion activation energy of water molecules result in higher transport diffusion ability of water molecules.Figure 11
(a) Pore size distributions of Nafion membranes at different temperatures.Vp is the pore volume, and dp is the pore size. (b) The volume and surface area of pores in Nafion membranes as a function of temperature.
(a)(b)The effect of temperature on proton conductivity is shown in Figure12(a). The proton conductivity is positively correlated with temperature, which is consistent with the correlation between the transport diffusion coefficient of hydronium ions and temperature. To quantify the number of hydrogen bonds in the simulated systems at different temperatures, geometric criteria for determining hydrogen bond formation were used. Two oxygen atoms are thought to form a hydrogen bond via a hydrogen atom when the following two criteria are met simultaneously: (I) the distance between the oxygen atoms is smaller than 3.5 Å and (II) the OHO angle is greater than 120° [67]. The number of hydrogen bonds (hydronium ion acts as donor and water molecule/sulfonic acid group acts as acceptors) at different temperatures is displayed in Figure 12(b). As temperature rises from 300 K to 350 K, the number of hydrogen bonds between hydronium ions and water molecules and that between hydronium ions and sulfonic acid group decreases. Elevated temperature allows the hydronium ions to break existing hydrogen bonds between hydronium ions and water molecules/sulfonic acid groups and readily form new hydrogen bonds with other molecules, which accelerates hydrogen bond dynamics and promotes proton transfer [69, 72]. Moreover, hydronium ions are more active at high temperatures, which increases the probability of them entering the free region (the middle region of the simulation box in Figure 3) through the PEM. Interactions between molecules in the free region are weakened at high temperatures, resulting in a reduced number of hydrogen bonds. High temperature accelerates the passage of protons through PEMs, indicating that it is beneficial to improve the performance of PEMs at a specific moisture content.Figure 12
(a) The proton conductivity and (b) the number of hydrogen bonds (hydronium ion acts as donor, and water molecule/sulfonic acid group acts as acceptors) at different temperatures.
(a)(b)


### 3.2.3. Comparison with Existing Works

In most molecular dynamics simulations, the self-diffusion coefficients of hydronium ions and water molecules in PEMs were calculated from the Einstein equation based on the mean square displacement [2, 7, 19, 20]. Self-diffusion was used to characterize the diffusion capacity of molecules at equilibrium (i.e., without driving forces such as the gradients in pressure and chemical potential). However, this paper focuses on the transport diffusion of hydronium ions and water molecules in PEMs under pressure differences, and the transport diffusion coefficients of molecules were calculated and discussed. It is found that the transport diffusion coefficients of water (2.252×10−9 m2/s) and hydronium ion (0.309×10−9 m2/s) at 300 K in this work are higher than the self-diffusion coefficients of water (0.591×10−9 m2/s) and hydronium ion (0.158×10−9 m2/s) in the simulation of Jang et al. under the same condition [26]. Moreover, the transport diffusion coefficient of water (2.252×10−9 m2/s) in this study is about an order of magnitude higher than the Fickian diffusion coefficient of water (~1×10−10 m2/s) calculated by Li et al. [34]. In addition, the diffusion activation energy of water (Table 2) in this study (20.25 kJ/mol) is consistent with the experimental results of Yeo et al. (20.2 kJ/mol) and Takamatsu et al. (20.51 kJ/mol) [73, 74]. The diffusion activation energy of hydronium ion (Table 3) in this study (21.15 kJ/mol) is close to the experimental results of Thompson et al. (21.23 kJ/mol) and Cappadonia et al. (21.22 kJ/mol) [75, 76]. The diffusion activation energy of hydronium ion in the experiment of Halim et al. (13.5 kJ/mol) is relatively low [77], which is due to the high water content (λ=23) of the membranes in their experiment. The diffusion activation energy generally decreases with increasing water content [14]. When the water content in PEM is high, the diffusion activation energy of hydronium ion is close to that of hydronium ion in bulk water (8.37–12.56 kJ/mol) [6].Table 2
The diffusion activation energy of water in the Nafion membrane.
EH2O (kJ/Mol)ConditionReference~23298.15K≤T≤373.15KThe experiment of Morris and Sun [78]20.51293.15K≤T≤354.15KThe experiment of Takamatsu et al. [74]20.20K≤T≤382.15KThe experiment of Yeo and Eisenberg [73]20.25300K<T<350KThis workTable 3
The diffusion activation energy of hydronium ion in Nafion membrane.
EH3O+ (kJ/Mol)ConditionReference22.6Dried sample,295.15K≤T≤335.15KThe experiment of Ye et al. [79]13.5λ=23,293.15K≤T≤363.15KThe experiment of Halim et al. [77]21.23λ=3.7,253.15K<T<353.15KThe experiment of Thompson et al. [75]21.224.2 g H2O/100 g Nafion, 225 K < T <300 KThe experiment of Cappadonia et al. [76]21.15λ=10,300K<T<350KThis work


## 3.2.1. Effect of Pressure Difference on Transport Diffusion Characteristics

To study the effect of pressure difference on the transport diffusion characteristics of hydronium ions and water molecules, different driving forces (0.15 kcal/(mol·Å)–0.45 kcal/(mol·Å)) were applied to hydronium ions and water molecules at 300 K. Water molecules and hydronium ions transported from anode to cathode through the Nafion 117 membrane under these driving forces. The transport diffusion coefficients of hydronium ions and water molecules are calculated based on Equation (2) and are given in Figure 8. As the driving force rises, the transport diffusion coefficient of water molecules increases from 0.88×10−9 m2/s to 3.36×10−9 m2/s, and that of hydronium ions also increases monotonically (0.15×10−9m2/s–0.46×10−9 m2/s). This trend is in line with that in the experiment of Husar et al. [30]. They studied water transfer through PEMs under pressure differences and found that the water mass-transfer flux increased with increasing pressure drop [30]. The higher pressure difference makes it easier for hydronium ions and water molecules to overcome the energy barrier and transport through the Nafion membrane. Additionally, the transport diffusion coefficient of water molecules exceeds that of hydronium ions under driving forces, which is similar to the simulation results displayed by Chen et al. [24]. The lighter water molecules have higher velocity and better transport diffusion capacity than hydronium ions. Moreover, a hydronium ion has an entirely positive charge, while a water molecule is neutral. The movement of hydronium ions is more susceptible to the negatively charged sulfonic acid group. The potential energy between oxygen atom in sulfonic acid group (Os) and hydronium ion (-28.37 kJ/mol) is lower than that between Os and water molecule (-13.17 kJ/mol) based on our statistics, which supports the above point. As a result, the migration capacity of hydronium ions in PEMs is lower than that of water molecules.Figure 8
The transport diffusion coefficients of hydronium ions and water molecules under different driving forces.Proton conductivity is a significant index to evaluate PEM performance. According to Equation (4), the proton conductivity is calculated and given in Figure 9. The proton conductivity increases from 0.005 S/cm to 0.015 S/cm as the driving force rises from 0.15 kcal/(mol·Å) to 0.45 kcal/(mol·Å). The proton conductivity in the simulations is lower than that obtained by the empirical Equation (6) (0.046 S/cm) [68]. This is because the water content of the PEM is lower than that of the system in simulations. The proton conductivity and transport diffusion coefficient of hydronium ions have the same trend with driving force. It indicates that the performance of PEMFCs could be improved by keeping the pressure difference between the anode and cathode in a reasonable range. For an actual PEMFC system, the pressure difference is generally less than 30 kPa [17].
(6)σ=5.139λ−3.261000exp1268303−1268T.Figure 9
The proton conductivity under different driving forces.


## 3.2.2. Effect of Temperature on Transport Diffusion Characteristics

The hydrogen bond dynamics and proton transfer depend on temperature [24, 69]. The EF-NEMD simulations of water molecule and hydronium ion transport through Nafion membranes were performed under a constant driving force (0.35 kcal/(mol·Å)) and different temperatures (300 K–350 K), and the simulation results are given in Figure 10(a). It shows that the transport diffusion coefficients of hydronium ions (0.31×10−9m2/s–0.95×10−9 m2/s) and water molecules (2.25×10−9m2/s–6.89×10−9 m2/s) both increase as temperature rises. In addition, the transport diffusion coefficient of water molecules is significantly higher than that of hydronium ions.Figure 10
(a) The transport diffusion coefficients of hydronium ions and water molecules under different temperatures. (b) The Arrhenius plots of ln (D) against 1000/T.
(a)(b)At an elevated temperature, the energy of the molecules increases, making it easier for water and protons to cross the energy barrier and transport through PEMs [24]. High temperature also changes the pore structure and thus affects the diffusivity of water molecules and hydronium ions within the proton exchange membrane. The PoreBlazer v4.0 tool was utilized to accurately measure the pore characteristics of Nafion membranes at different temperatures [70, 71]. The pore size distribution, and the pore volume and surface area of Nafion membranes (Figure 11) were calculated by a helium probe. It can be seen that the pore size distribution in Nafion membranes exhibits a bimodal shape (Figure 11(a)). The two peaks shift to the right with increasing temperature, indicating that the pore size increases at high temperatures. As the temperature rises, the volume and surface area of pores also increase (Figure 11(b)). The enlarged pore volume at high temperatures also lowers the energy barrier for water molecules and hydronium ions across PEMs, thereby increasing their transport diffusion coefficients. Moreover, the Arrhenius plots of transport diffusion coefficients for hydronium ions and water molecules are presented in Figure 10(b). According to Equation (3), the diffusion activation energies of hydronium ions and water molecules are 21.15 kJ/mol and 20.25 kJ/mol, respectively. The lower molecular weight and diffusion activation energy of water molecules result in higher transport diffusion ability of water molecules.Figure 11
(a) Pore size distributions of Nafion membranes at different temperatures.Vp is the pore volume, and dp is the pore size. (b) The volume and surface area of pores in Nafion membranes as a function of temperature.
(a)(b)The effect of temperature on proton conductivity is shown in Figure12(a). The proton conductivity is positively correlated with temperature, which is consistent with the correlation between the transport diffusion coefficient of hydronium ions and temperature. To quantify the number of hydrogen bonds in the simulated systems at different temperatures, geometric criteria for determining hydrogen bond formation were used. Two oxygen atoms are thought to form a hydrogen bond via a hydrogen atom when the following two criteria are met simultaneously: (I) the distance between the oxygen atoms is smaller than 3.5 Å and (II) the OHO angle is greater than 120° [67]. The number of hydrogen bonds (hydronium ion acts as donor and water molecule/sulfonic acid group acts as acceptors) at different temperatures is displayed in Figure 12(b). As temperature rises from 300 K to 350 K, the number of hydrogen bonds between hydronium ions and water molecules and that between hydronium ions and sulfonic acid group decreases. Elevated temperature allows the hydronium ions to break existing hydrogen bonds between hydronium ions and water molecules/sulfonic acid groups and readily form new hydrogen bonds with other molecules, which accelerates hydrogen bond dynamics and promotes proton transfer [69, 72]. Moreover, hydronium ions are more active at high temperatures, which increases the probability of them entering the free region (the middle region of the simulation box in Figure 3) through the PEM. Interactions between molecules in the free region are weakened at high temperatures, resulting in a reduced number of hydrogen bonds. High temperature accelerates the passage of protons through PEMs, indicating that it is beneficial to improve the performance of PEMs at a specific moisture content.Figure 12
(a) The proton conductivity and (b) the number of hydrogen bonds (hydronium ion acts as donor, and water molecule/sulfonic acid group acts as acceptors) at different temperatures.
(a)(b)


## 3.2.3. Comparison with Existing Works

In most molecular dynamics simulations, the self-diffusion coefficients of hydronium ions and water molecules in PEMs were calculated from the Einstein equation based on the mean square displacement [2, 7, 19, 20]. Self-diffusion was used to characterize the diffusion capacity of molecules at equilibrium (i.e., without driving forces such as the gradients in pressure and chemical potential). However, this paper focuses on the transport diffusion of hydronium ions and water molecules in PEMs under pressure differences, and the transport diffusion coefficients of molecules were calculated and discussed. It is found that the transport diffusion coefficients of water (2.252×10−9 m2/s) and hydronium ion (0.309×10−9 m2/s) at 300 K in this work are higher than the self-diffusion coefficients of water (0.591×10−9 m2/s) and hydronium ion (0.158×10−9 m2/s) in the simulation of Jang et al. under the same condition [26]. Moreover, the transport diffusion coefficient of water (2.252×10−9 m2/s) in this study is about an order of magnitude higher than the Fickian diffusion coefficient of water (~1×10−10 m2/s) calculated by Li et al. [34]. In addition, the diffusion activation energy of water (Table 2) in this study (20.25 kJ/mol) is consistent with the experimental results of Yeo et al. (20.2 kJ/mol) and Takamatsu et al. (20.51 kJ/mol) [73, 74]. The diffusion activation energy of hydronium ion (Table 3) in this study (21.15 kJ/mol) is close to the experimental results of Thompson et al. (21.23 kJ/mol) and Cappadonia et al. (21.22 kJ/mol) [75, 76]. The diffusion activation energy of hydronium ion in the experiment of Halim et al. (13.5 kJ/mol) is relatively low [77], which is due to the high water content (λ=23) of the membranes in their experiment. The diffusion activation energy generally decreases with increasing water content [14]. When the water content in PEM is high, the diffusion activation energy of hydronium ion is close to that of hydronium ion in bulk water (8.37–12.56 kJ/mol) [6].Table 2
The diffusion activation energy of water in the Nafion membrane.
EH2O (kJ/Mol)ConditionReference~23298.15K≤T≤373.15KThe experiment of Morris and Sun [78]20.51293.15K≤T≤354.15KThe experiment of Takamatsu et al. [74]20.20K≤T≤382.15KThe experiment of Yeo and Eisenberg [73]20.25300K<T<350KThis workTable 3
The diffusion activation energy of hydronium ion in Nafion membrane.
EH3O+ (kJ/Mol)ConditionReference22.6Dried sample,295.15K≤T≤335.15KThe experiment of Ye et al. [79]13.5λ=23,293.15K≤T≤363.15KThe experiment of Halim et al. [77]21.23λ=3.7,253.15K<T<353.15KThe experiment of Thompson et al. [75]21.224.2 g H2O/100 g Nafion, 225 K < T <300 KThe experiment of Cappadonia et al. [76]21.15λ=10,300K<T<350KThis work


## 4. Conclusions

In this study, the EF-NEMD simulations were performed to investigate the effects of pressure difference (i.e., driving force) and temperature on the transport diffusion of hydronium ions and water molecules through the Nafion 117 membrane. The following conclusions can be drawn:(1)
The influence of sulfonic acid groups on the transport diffusion of hydronium ions is greater than that of water molecules under pressure differences because of the stronger hydrogen bonding between sulfonic acid groups and hydronium ions(2)
The transport diffusion coefficients of hydronium ions and water molecules rise with the increase in pressure difference. Since the diffusion activation energy (EH2O=20.25 kJ/mol and EH3O+=21.15 kJ/mol) and molecular weight of water molecules are lower than that of hydronium ions, the transport diffusion coefficient of water molecules is higher than that of hydronium ions under different pressure differences(3)
The proton conductivity and the transport diffusion coefficient of water molecules are positively correlated with temperature. The pore size of the Nafion membrane enlarges, and the energy of small molecules increases at high temperatures. In addition, the increase in temperature speeds up the dynamics of hydrogen bonds between hydronium ions and water molecules/sulfonic acid groups, which further promotes proton transfer



---

*Source: 1138198-2023-02-07.xml*