Enhanced hydronium ion diffusion in proton exchange membranes reinforced with multilayer graphene oxide: new insights into water retention and ion mobility using molecular dynamics simulation

Abstract

Graphene oxide (GO) reinforced perfluorosulfonic acid (PFSA) based proton exchange membranes (PEMs) show enhanced ion diffusion resulting in elevated polymer electrolyte fuel cell (PEFC) performance. However, the mechanisms by which GO influences water dynamics and ion (hydronium) transport are relatively less explored in the literature. In addition, it is expected that the interlayer spacing of multilayer GO plays a crucial role in promoting ion mobility. To this end, this research article explores the possibility of providing new insights into the water/ion dynamics as well as identifying the impact of interlayer spacing of GO on the ion diffusion. Molecular dynamics (MD) simulation is implemented to elucidate the behaviour of multilayer-GO with PFSA structure and to examine the interactions between functional groups (epoxy and hydroxyl) on the GO surface with water molecules and hydronium ions. The retention of water molecules adjacent to the multilayer-GO plays a crucial role in forming transport channels that significantly enhance ion mobility within the membrane structure. The optimal interlayer spacing of 9.5 Å is identified as the critical threshold value where ion diffusion is observed at its peak. In comparison with pristine Nafion®, the ion (hydronium) diffusion coefficient in the multilayer-GO with PFSA polymer shows an improvement of ∼17% and ∼30% at 300 K and ∼9% and ∼12% at 350 K for hydration levels (λ) of 13 and 20, respectively.

---

Introduction

Polymer electrolyte fuel cells, also known as PEFCs, have garnered significant interest among the various fuel cells because of their notable environmental sustainability, lightweight nature, rapid initiation at low temperatures (80 °C), and prompt adaptability to power requirements.^1–4 Consequently, they have become appealing alternatives for portable electronic devices and automotive uses.^5–8 The durability and ionic conductivity of the PEMs, which serve as an ion exchange medium in a PEFC, are crucial factors affecting the fuel cell's efficiency.^9,10 Currently, perfluorosulfonic acid (PFSA) membranes and their commercial form Nafion® are widely used as PEMs in PEFCs. Nafion®, a leading PFSA ion-exchange polymer, was first developed by DuPont in the late 1960s. The PFSA membrane exhibits several features, including chemical stability, thermal stability, hydrophilicity, and strong ionic conductivity for ion transport.^11–13

In PFSA membranes, the ion conductivity is facilitated by the [–SO₃H] group attached to the side chain of the PFSA polymer. Ion mobility within the membrane is responsible for proton conductivity and is strongly influenced by the moisture/water content present in the membrane.^14–16 Numerous research studies have identified several limitations of PFSA membranes, including inconsistent ion conductivity, inferior mechanical properties under high-temperature conditions, degradation of properties in strong oxidative environments and high cost.^17–19 Therefore, PFSA membranes require modification to enhance ion conductivity, mechanical properties, and chemical and thermal stability, making them suitable for diverse fuel cell operating conditions.^20–22 The organic and inorganic reinforcement to these PEMs is a well-known method to enhance the properties including oxidative stability.²³ Different types of fillers, such as metal oxides (zirconia, titania, and silica), heteropoly acids, metal phosphates, clays, zeolites, and graphene oxide (GO), are integrated into PEMs to enhance their characteristics. These fillers enhance mechanical properties through their strong interaction with the polymer matrix, facilitate water absorption due to their hygroscopic characteristics, exhibit high proton conductivity, and provide chemical and thermal stability. These fillers create interlinked ion-conducting channels and maintain significant water content in the membrane, enhancing PEMs' proton conductivity.^24–27 GO, a two-dimensional carbon-based material, has emerged as a promising filler material for PEMs. Numerous studies have explored different forms of GO, including reduced graphene oxide (rGO), oxidized GO, and sulfonated graphene oxide (SGO). Sulfonated graphene oxide (SGO) is formed by introducing sulfonic acid (–SO₃H) groups, enhancing its dispersibility in aqueous solutions and improving its ion exchange capacity, making it an ideal candidate for applications in fuel cells, supercapacitors, and water purification membranes. The oxidized form of GO contains a high concentration of oxygen functionalities, enhancing its hydrophilicity, water retention, and overall mechanical stability, making it a promising option for improving fuel cell performance, particularly as a membrane material in PEFCs.^28–33 It contains oxygen-rich functional groups, including epoxy, hydroxyl, and carboxyl groups. These hydrophilic functional groups enable GO to efficiently absorb water molecules between its layers, with some water molecules remaining trapped in its interlayer spaces even after prolonged drying.^34,35 The robust interaction between the functional groups on the GO layer and water molecules makes it difficult to separate GO from the system.³⁶

On the other hand, the interlayer spacing of GO is a critical parameter in determining the collective functional properties of graphene-based materials, particularly in membrane-based applications, catalysis, and energy storage. Experimental techniques such as X-ray diffraction (XRD) and transmission electron microscopy (TEM) have been widely used to characterize interlayer distances. Studies indicate that oxidation increases the interlayer spacing in GO compared to pristine graphite, primarily due to the introduction of oxygen-containing functional groups. Conversely, reduction processes, such as thermal or chemical treatments, remove oxygen functionalities, leading to a decrease in interlayer spacing in rGO.^37–40 The carbon-to-oxygen (C/O) ratio plays a crucial role in modulating the interlayer spacing of GO-based materials. A lower C/O ratio (higher oxygen content) increases the interlayer spacing. Conversely, as the C/O ratio increases (due to reduction processes), the removal of these functional groups leads to a decrease in interlayer spacing, making rGO more compact compared to GO. These findings emphasize the importance of controlling the oxidation level to tune the interlayer spacing for specific applications. The stacked nanostructure comprises 6 to 7 graphene layers with an average interlayer spacing of 9 Å.²⁸ The interlayer spacing decreases with the reduction of GO, typically to 4 Å depending on the extent of the reduction. Moreover, environmental factors like humidity and temperature also influence the interlayer spacing. High humidity can expand the interlayer spacing to 12–14 Å due to water intercalation, while at a high temperature (1000 K), the spacing can increase by approximately 6.63% to 16.43%, depending on the oxidation.^37,41–46

Experimentally, several notable studies have been carried out to investigate the properties of GO with PFSA membranes. Wang et al.⁴⁷ developed a composite membrane for fuel cell applications by incorporating custom-synthesized GO into Nafion® resin. GO is used as a reinforced filler to improve the mechanical properties of Nafion®. This finding compares the fuel cell performance of a 3 wt% GO/Nafion® composite with pure Nafion®, demonstrating that the composite membrane showed superior mechanical properties. Choi et al.⁴⁸ introduced the synthesis of composite membranes made of Nafion® and GO that can be utilised in direct methanol fuel cells (DMFCs). This study showed GOs’ amphiphilic characteristics resulting from hydrophobic and hydrophilic functional groups. Cha et al.⁴⁹ reduced the cost and improved the performance of PEFCs by impregnating a porous polyethylene substrate with Nafion®, obtaining a reinforced composite membrane with a two-layered asymmetric structure. The reinforced composite membrane shows improved ion conductivity after annealing at 150 °C because of the crystallinity of the Nafion® ionomer. It improves mechanical strength due to the contact between the PE substrate and the Nafion® ionomer. Bayer et al.⁵⁰ examined the electrical and mechanical properties of reinforced membranes with GO and compared them with conventional Nafion® membranes. The findings revealed that the reinforced GO membrane exhibited superior water retention, which is beneficial for improving proton conductivity.

Nevertheless, a limited number of research studies have utilised MD simulation techniques to comprehend the behaviour of reinforced PFSA membranes with GO. MD simulation effectively analyses the interactions between molecules and the system's mechanical, dynamic, and structural properties.⁵¹ Maiti et al.⁵² performed the MD simulation method and showed that adding sulfonic acid (–SO₃H) functionalized graphene oxide (SGO) to PFSA membranes greatly improves their mechanical characteristics, proton conductivity, and glass transition temperature (T_g). The results show that the Nafion/SGO composite membranes perform better in fuel cells. This is mostly because of the strong contacts between SGO and PFSA, which reduce fuel permeability and produce interconnected proton-conducting channels. Kritikos et al.⁵³ observed that hydronium ions are partially adsorbed onto the GO surface, with the adsorption/desorption rate notably increasing with higher hydration levels (hereafter denoted as λ). The translational dynamics of water molecules near the GO surface were significantly slower compared to those at greater distances from the GO surface. Tanaka et al.⁵⁴ explained that the interfacial dynamics between Nafion® and graphene sheets are critical in influencing proton transport processes within Nafion®/carbon nanotube (CNT) composite membranes. The protons near the interface exhibit higher self-diffusion coefficients than those in the bulk Nafion® region across varying water contents. This finding highlights the significance of interfacial water layer development in affecting transport properties. The results suggest that incorporating high-specific surface area CNTs enhances proton transport by optimizing interfacial structures, offering potential improvements in the performance of PEMs. Using MD simulations, Devanathan et al.⁵⁵ examined the interactions of water molecules with GO at different λ values. This work revealed that hydrogen-bonding interactions between the water molecules and the hydroxyl groups of GO cause the diffusion of water molecules in GO to be slowed. While numerous experimental and a few MD simulation studies have focused on Nafion® reinforced with GO composites, the potential of multilayer GO combined with PFSA, featuring varying interlayer spacings, remains an intriguing and largely unexplored research avenue. In this study, MD simulation is conducted to investigate the behaviour of reinforced PFSA (Nafion®) with multilayer GO across varying temperatures and λ values and examine the diffusion of hydronium ions.

2. Computational framework

This study considers a ten-chain Nafion® model with explicit water molecules and hydronium ions to construct the simulation cell. The terminology used to classify simulation cells of Nafion® with multilayer-GO is Nafion/multilayer-GO, where the multilayer-GO loading is 3 wt%. The interlayer spacing (z) of multilayer-GO varies from 7.5 Å to 10.5 Å.^53,55 The λ decides the number of H₂O molecules within each cell, where λ is the ratio of H₂O molecules and (H₃O⁺) hydronium ions to the (SO₃⁻) sulfonic groups.^56,57 The molecular composition of the hydronium ion molecule, water molecule, and Nafion® chain is available in SI Fig. S1. The molecular structural detail of multilayer-GO is shown in Fig. S2. A simulation cell is prepared for the Nafion/multilayer-GO system at varying λ values from 3.5 to 20. Details of H₂O molecules and hydronium ions in the simulation cell at different λ values are presented in Table 1.

Table 1 Details of H₂O molecules and hydronium ions (H₃O⁺) in the simulation cell at varying λ values

λ	Number of hydronium ions	Number of water molecules
3.5	100	250
7	100	600
10	100	900
13	100	1200
16	100	1500
20	100	1900

---

The constituent description of multilayer-GO is shown in Fig. 1. The lateral dimension of the multilayer-GO layer is 14.27 × 16.1 Ų. Multilayer-GO structures exhibit interlayer spacings ranging from 7.5 Å to 10.5 Å, as illustrated in Fig. 1. The ratio of carbon atom to oxygen atom is 5 : 1 and the epoxy group to hydroxyl group is 3 : 2.^53,58,59 The hydroxyl and epoxy groups are randomly distributed. Hydrogen atoms are bonded to the carbon atoms at the edge.⁵³ The multilayer-GO of 3 wt% has 225 carbon atoms, 27 hydroxyl groups, and 18 epoxy groups. The initial configuration of Nafion/multilayer-GO is generated using the Avogadro software.⁶⁰ Next, these molecules are randomly placed in a cubic box using the Packmol software⁶¹.

Fig. 1 Details of the multilayer-GO structure with interlayer spacing (z) ranging from 7.5 to 10.5 Å. The dotted circle highlights the epoxy and hydroxyl functional groups.

The present investigation uses all-atom force fields to conduct classical MD simulations on Nafion®, water, hydronium ions and multilayer-GO. Simulations are performed with the modified Dreiding force field.^62,63 The modified Dreiding force field signifies a notable progression in the realm of ab initio force fields; thus, it provides the precise estimation of various properties in both the gas phase (such as structural, conformational, and vibrational properties) and the condensed phase for a diverse array of molecules. Additionally, this force field represents a pioneering advancement in integrating factors pertaining to organic and inorganic materials.

A modified Dreiding force field's functional form can be categorised into bonded and non-bonded interaction terms. The bonded interaction is composed of a bond-stretching term (E_bond), the angle-bending term (E_angle), and the torsion term (E_torsion) and non-bonded interactions include the van der Waals (E_vdW) term and the electrostatic term (E_{electrostatic}) as explained in eqn (1). Interatomic interactions pertain to the connections between pairs of atoms that are positioned at a distance of two or more atoms in between or are part of separate molecules.

E_total = E_bond + E_angle + E_torsion + E_vdW + E_{electrostatic} (1)

Eqn (1) is expressed in the functional form as

where K_b, K_θ, and V_n are the bond, angle, and torsional force constants, respectively, r₀ and θ₀ are the equilibrium bond lengths and bond angles, respectively, n is the multiplicity, and ϕ is the phase angle for torsional parameters. Regarding van der Waals interaction, ε is the van der Waals well depth (kcal mol⁻¹), σ represents the point at which the potential energy between two particles becomes zero, and r denotes the distance separating these interacting particles, and regarding coulombic interaction, q is the partial charge on atoms as described in eqn (2). The partial charges for the Nafion® chain are adopted from the work of Cha et al.,⁵⁶ while those for the multilayer-GO are taken from Kritikos et al.⁵³ The details of partial charges are presented in Table S1. All interaction parameter-related information leveraged for the Nafion® membrane and multilayer-GO is available in the SI (Tables S2–S5). The F₃C three-site H₂O model is used for water molecules,⁵⁶ and the classical hydronium model⁵⁶ is used for hydronium molecules. For the Lennard-Jones interactions (van der Waals), the Lorentz–Berthelot mixing rule⁶⁴ is adopted to describe how various atom types interact with one another eqn (3). The particle–particle–particle–mesh method (PPPM)⁶⁵ technique is used to handle all coulombic interactions, with a specified cut-off distance of 15 Å.

The simulation cell must be equilibrated to determine the structural and dynamic properties. To attain equilibration, the subsequent steps are as follows: (i) minimisation, (ii) MD simulation for 100 ps at 300 K in an NVT ensemble, (iii) MD simulation for 100 ps at 300 K in an NPT ensemble, (iv) heating of the structure from 300 K to 600 K, for 100 ps in the NVT ensemble, (v) then cooling of the structure from 600 K to 300 K, for 100 ps in the NVT ensemble, and (vi) final MD simulation for 100 ps at 300 K in the NPT ensemble. The iterative process outlined in steps (iv) to (vi) is continued until the target density is achieved.⁶⁶ The resulting density values for λ = 3.5, 7, 10, 13, 16, and 20 are 1.71, 1.68, 1.65, 1.62, 1.58, and 1.56 g cm⁻³, respectively, and are close to the previous studies.^66–68 The simulation box size of the Nafion/multilayer-GO model after equilibration for λ = 20 is 55 × 55 × 55 ų. This accommodates a maximum of three layers of GO. The multilayer-GO is self-tethered throughout the equilibration phase with a spring force of 2 kcal mol⁻¹. This approach maintains realistic interfacial interactions between the multilayer-GO, Nafion®, water molecules, and hydronium ions. The simulation employs three-dimensional periodic boundary conditions with one femtosecond (fs) time step. After achieving equilibration, the system proceeds to the production run phase for 1 ns NPT simulation at 300 K and 350 K, respectively. To ensure statistical convergence and reliability, five independent MD simulations are conducted. Ensemble averaging is performed over these five unique runs to obtain statistically meaningful results. This strategy effectively reduces fluctuations and enhances the accuracy of computed properties such as the radial distribution function (RDF), mean square displacement (MSD), and diffusion coefficient. All MD simulation cases in this work are run in the LAMMPS (large-scale molecular massively parallel simulator) code from Plimpton at Sandia.⁶⁹ The Verlet algorithm is employed to integrate the equations of motion.⁷⁰ Nose–Hoover thermostats are used to maintain the desired temperature, and Nose–Hoover barostats^71,72 are used to maintain the atmospheric pressure during simulation.

Examining configurations derived from MD simulations plays a vital role in understanding the microstructure of hydrated membranes. One method to explore the hydrated membrane's microstructure is RDF. The RDF is a mathematical representation of the likelihood of locating an atom P at a specific distance r from a reference atom A. This probability can be determined using eqn (4).

In eqn (4), V denotes the overall system volume, N_p denotes the total number of P atoms in the simulation box, and n_p is the number of P atoms located at r distance from particle A in a dr-thick shell.

The transportation of water and hydronium ions plays a vital role in the functioning of a PEM. Estimating self-diffusion coefficients in classical MD simulations involves an analysis of MSDs.^56,70,73 The mathematical description of MSD follows eqn (5).

In eqn (5), c_i (t) signifies the positional coordinates of atom i at time t, whereas N denotes the number of atoms undergoing free diffusion. As per Einstein's diffusion law, the self-diffusion coefficient (D) can be determined using eqn (6).⁷⁰

3. Results and discussion

Fig. 2 shows the simulated results of the diffusion coefficient of hydronium ions (D_h) for varying interlayer spacings of GO at λ = 10 and temperatures of 300 K and 350 K. The results show that the highest ion diffusion coefficient occurs at an interlayer spacing of 9.5 Å. This enhancement at 9.5 Å can be attributed to an optimal interlayer transport channel formation within the multilayer-GO, where favourable water structuring facilitates ion mobility. For interlayer spacings below 9.5 Å, the confinement effect restricts the mobility of water and hydronium ions due to reduced free volume and stronger interfacial interactions with GO functional groups, which hinder ion transport. Conversely, when the interlayer spacing exceeds 9.5 Å, the degree of confinement decreases, leading to weaker water–GO interactions and reduced water retention near the GO layers. Notably, this trend is consistently observed at 300 K and 350 K, indicating that the optimal interlayer spacing remains consistent across different temperatures. Therefore, an interlayer spacing of 9.5 Å will be considered for upcoming analyses.

Fig. 2 Diffusion coefficient of hydronium ions (D_h) at λ = 10 for the Nafion/multilayer-GO system at temperatures 300 K and 350 K.

3.1 Radial distribution function (RDF)

To understand the structural organization of water molecules and hydronium ions around the sulfonic acid groups, RDF analysis is conducted between sulphur atoms (from the sulfonic groups) and water oxygen atoms (S–O_w), as well as sulphur and hydronium oxygen atoms (S–O_h). Fig. 3 and 4 show the RDF between S–O_w and S–O_h for Nafion® and Nafion/multilayer-GO systems at 300 K and 350 K, respectively, for varying λ = 3.5 to 20. One may not observe a significant difference in Fig. 3 or 4. Therefore, RDF may not provide exclusive reasons for enhancing the diffusion properties due to reinforcements. However, a trend is observed as a function of temperature: the RDF peak of S–O_w decreases, whereas that of S–O_h increases (for higher λ), for Nafion® and Nafion/multilayer-GO, as the temperature increases from 300 K to 350 K. At 350 K, the RDF peak for S–O_w decreases as shown in Fig. 4(a) and (b), indicating that the hydrogen bonding and structural ordering of water molecules around the sulfonate groups become weaker due to the enhanced motion of water molecules at high temperatures. This reduction in water residence time near the sulfonate sites leads to a less pronounced RDF peak. Conversely, the S–O_h RDF peak increases with temperature as shown in Fig. 4(c) and (d), suggesting that hydronium ions exhibit a stronger localisation around the sulfonate groups when the sulfonate–water interactions are weakened. This phenomenon highlights that water and hydronium ions interact with sulfonate sites, which change significantly at higher temperatures.

Fig. 3 RDF between sulphur atoms (S) and water oxygen atoms (O_w) for (a) Nafion and (b) the Nafion/multilayer-GO system at 300 K. RDF between sulphur atoms (S) and hydronium oxygen atoms (O_h) for (c) Nafion® and (d) the Nafion/multilayer-GO system at 300 K.

Fig. 4 RDF between sulphur atoms (S) and water oxygen atoms (O_w) for (a) Nafion and (b) the Nafion/multilayer-GO system at 350 K. RDF between sulphur atoms (S) and hydronium oxygen atoms (O_h) for (c) Nafion® and (d) the Nafion/multilayer-GO system at 350 K.

The RDF analysis between hydroxyl oxygen atoms (O_Hydroxyl) of multilayer-GO and water oxygen atoms (O_Water) at varying λ values and temperatures provides valuable insights into the local water structuring near functionalized GO surfaces. The RDF is plotted for higher λ = 10, 13, 16, and 20 excluding λ = 3.5 and 7 due to insufficient water molecules. At low hydration, the limited presence of water around hydroxyl groups limits meaningful statistical analysis. As shown in Fig. 5(a), the peak position is around 3.5 Å at λ = 10, while it shifts to approximately 3.1 Å at λ = 20. The peak appears at a shorter distance under higher hydration conditions due to the increased number of water molecules, which leads to a denser and more closely packed local structure near the hydroxyl groups. However, as λ increases, a systematic decrease in peak intensity, along with slight broadening and a shift, is observed, indicating a more localized water structure at high λ. In contrast, the reduced water content at lower λ leads to a more dispersed arrangement, causing the peak to appear at a longer distance with increased intensity. A similar trend is observed at 350 K, as shown in Fig. 5(b). Overall, these findings emphasize that at low hydration, fewer water molecules are available near the hydroxyl group, but at higher λ, the water distribution becomes more homogeneous, which may influence the water retention and ion transport properties.

Fig. 5 RDF between hydroxyl oxygen (O_Hydroxyl) and water oxygen (O_Water) for Nafion/multilayer-GO at (a) 300 K and (b) 350 K.

In addition, the interaction energy (IE) between graphene layers and water molecules, as well as between functional groups (epoxy and hydroxyl) and water molecules, was computed at 300 K and 350 K, including both van der Waals and electrostatic interactions. Negative interaction energy values signify attractive forces between the molecules, with increasingly negative values indicating stronger intermolecular attraction. The total interaction energy is computed using the compute group/group command in LAMMPS. The findings suggest that when the number of water molecules in the simulation cell increases, the interaction energy becomes more negative, indicating strong water retention properties near the functionalized graphene. However, the magnitude of the interaction energy was slightly reduced at 350 K compared to 300 K, due to the increased movement of water molecules at higher temperatures, which weakens the strength of hydrogen bonding and van der Waals forces. This temperature-dependent behaviour suggests that while water retention near the functionalized graphene is robust, thermal agitation at high temperatures slightly reduces the strength of water–graphene interactions. In Fig. 6, the inset bar graphs present the interaction energies for different λ values at 300 K and 350 K, clearly distinguishing graphene–water and functional group–water interactions. Fig. 6 also shows the snapshots of the molecule configuration, including multilayer-GO and water molecules. As the λ increases, water molecules near the GO layer increase.

Fig. 6 Interaction energy (IE) between Nafion/multilayer-GO and water molecules, functional groups (epoxy and hydroxyl) and water molecules at λ = 3.5 to 20 for temperatures 300 K and 350 K.

3.2 Cluster analysis

A prior study indicates that the cluster of H₂O molecules depends on the λ of the Nafion® simulation cell; as the λ increases, the cluster size also increases.^57,74 At low λ values, such as 3.5 and 7, the connectivity between clusters is poor; however, at higher λ values, the connectivity between clusters improves, forming transport channels that facilitate the movement of hydronium ions. Therefore, it is interesting to study specific aspects of this cluster formation behaviour as a function of λ and temperature. For this purpose, ‘agglomerative clustering’ is used to analyze the cluster formation of water molecules near the multilayer-GO. It is the type of hierarchical clustering that works by iteratively merging small clusters into larger ones. This method begins by treating each data point as an individual cluster; the two closest clusters are merged at each step. The process continues until all points are part of a single cluster or until a specified stopping criterion, such as a distance threshold, is met. Some parameters are defined before implementing this algorithm, such as metric, distance threshold, linkage, etc. In this study, the ‘metric = Euclidean’ metric is used to compute the pairwise distance between data points. The ‘linkage = single’ (distance between two clusters) is defined as the shortest distance between any pair of points in the two clusters. For a distance threshold = 4 Å, the maximum distance between clusters will be merged, beyond which the algorithm will stop merging clusters.

The ‘agglomerative clustering’ is applied to the simulation data obtained from the production run stage at a particular time step, at varying λ = 10, 13, 16 and 20 for 300 K and 350 K. Fig. 7 compares the cluster of H₂O molecules in a Nafion/multilayer-GO simulation cell as well as in a Nafion® simulation cell at 300 K. The size of the water clusters in the Nafion/multilayer-GO is larger than that of the water clusters in Nafion®. This difference is consistent across various λ values. In the Nafion® case, the small clusters are more in number, whereas in the Nafion/multilayer-GO, larger clusters form instead of smaller ones. At 350 K, a similar phenomenon of water cluster formation is observed, as illustrated in Fig. 7(c) and (d) for the Nafion/multilayer-GO system and Nafion®, respectively. The main observation is that the size of the cluster somewhat decreases at 350 K compared to the 300 K case because of increased water molecule mobility at higher temperatures.

Fig. 7 Water cluster formation at varying λ = 10, 13, 16 and 20 for (a) Nafion/multilayer-GO and (b) Nafion® at 300 K and for (c) Nafion/multilayer-GO and (d) Nafion® at 350 K.

Furthermore, the critical observation is that the mobility of water molecules near the GO layer is much less than that of the water molecules away from the GO layer. The functional groups of GO help bind the water molecules together due to its hydrophilic nature.⁷⁵Fig. 8(a) shows two regions in the simulation cell. Region 1 covers the volume near the multilayer-GO and region 2 covers the volume far from the multilayer-GO. Fig. 8(b) shows the total path traversed by water molecules for a particular time, 0.25 ns. For this study, out of a total of 80 water molecules, 40 are considered in region 1 and 40 in region 2. The total path traversed by all 40 water molecules in region 1 covered much less distance in 0.25 ns than the other 40 water molecules in region 2.

Fig. 8 (a) Schematic representing region 1 (near the multilayer-GO) and region 2 (far from the multilayer-GO) inside the simulation box. (b) Representation of the total path traversed by water molecules in 0.25 ns in region 1 and region 2.

3.3 Diffusion coefficient

The robustness of the Nafion® model is tested for the diffusion coefficient of hydronium ions (D_h) at various hydration levels and validated against previously reported studies.^66,76Fig. 9 shows the comparative values of D_h for the Nafion® model. For this work, the D_h values obtained (in the units of ×10⁻⁵ cm² s⁻¹) of Nafion® at 300 K for λ = 10, 13, 16 and 20 are 0.10066, 0.24602, 0.3067 and 0.40232, respectively. As observed from the figure, these values are in good agreement with the simulation results for the Nafion® two hundred-chain model⁷⁶ and four-chain model.⁶⁶

Fig. 9 Diffusion coefficient of hydronium ions (D_h) for varying λ values as obtained by the present work (triangle), Tsung et al.⁷⁶ (square), and Devanathan et al.⁶⁶ (circle), respectively.

Fig. 10(a) shows the D_h at 300 K for Nafion® and Nafion/multilayer-GO, with error bars indicating the variation from the mean values. The D_h values (in the units of ×10⁻⁵ cm² s⁻¹) of Nafion/multilayer-GO at 300 K for λ = 10, 13, 16 and 20 are 0.23417, 0.29018, 0.36915 and 0.52272, respectively. Similarly, Fig. 10(b) shows the D_h at 350 K. The D_h values(in the units of ×10⁻⁵ cm² s⁻¹) of Nafion® at 350 K for λ = 10, 13, 16 and 20 are 0.21115, 0.4004, 0.61525 and 0.7945 respectively. The D_h values (in units of ×10⁻⁵ cm² s⁻¹) of Nafion/multilayer-GO at 350 K for λ = 10, 13, 16 and 20 are 0.36604, 0.43507, 0.68235 and 0.89183, respectively. The D_h at 350 K is significantly higher than that at 300 K, primarily due to the enhanced mobility of hydronium ions at higher temperature. The D_h values at 300 K and 350 K for Nafion® and Nafion/multilayer-GO at different λ values are detailed in Tables S6 and S7.

Fig. 10 Diffusion coefficient of hydronium ions (D_h) of Nafion® and Nafion/multilayer-GO at (a) 300 K and (b) 350 K. The error bars represent the variation from the mean value.

For low λ = 3.5 and 7 at 300 K, the D_h in the Nafion/multilayer-GO system is observed to be lower than that of Nafion®. This reduction in hydronium ion mobility can be attributed to the lack of well-connected transportation channels at low water content. A similar trend is observed at 350 K, where the D_h remains lower for the Nafion/multilayer-GO system at low λ, confirming that insufficient water connectivity effects persist even at higher temperatures. However, at high λ = 10, 13, 16, and 20, the D_h of the Nafion/multilayer-GO is higher than that of Nafion® at 300 K and 350 K. The Nafion/multilayer-GO system exhibits a significantly enhanced diffusion coefficient compared to Nafion® at temperatures of 300 K and 350 K. This improvement is attributed to the multilayer-GO, which promotes water retention and facilitates the formation of an efficient transport channel, thereby enhancing the mobility of hydronium ions within the membrane structure.

4. Conclusion

In this study, MD simulations are conducted on a reinforced PFSA membrane, specifically Nafion® with multilayer-GO, to investigate the effect of hydration on hydronium ion diffusion. The optimal interlayer spacing for enhancing ion diffusion is identified as 9.5 Å. The presence of functional groups on the GO surface, due to their hydrophilic nature, promotes water retention by holding water molecules close to the multilayer-GO. The same phenomenon is observed in the RDF between O_Hydroxyl and O_Water. This phenomenon led to the formation of continuous water channels, which are vital for efficient hydronium ion transport. Moreover, the agglomerative clustering algorithm shows that the Nafion/multilayer-GO system forms larger water clusters than Nafion®. The clustering effect is more pronounced at higher λ, where the interaction between the water molecules and the hydrophilic functional groups on the multilayer-GO is most intense. These larger water clusters near the GO layers contributed to water retention and created a stable environment for hydronium ion transport. Interestingly, the water molecules near the GO layers exhibited minimal movement compared to those away from the multilayer-GO. Furthermore, the D_h at 300 K and 350 K is higher in the Nafion/multilayer-GO system compared to that in Nafion®. Specifically, the D_h in the Nafion/multilayer-GO system is ∼132% and ∼73% higher than that of Nafion® at λ = 10 and 300 K and ∼20% and ∼11% higher at λ = 16 and 350 K. This demonstrates the significant enhancement in ion mobility due to the multilayer-GO reinforcement, particularly at high λ. These findings offer valuable insights into designing high-performance PEMs for PEFCs and could be extended to other ion transport applications.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information: Dreiding force field parameter information and the diffusion coefficient of hydronium ions. See DOI: https://doi.org/10.1039/d5ma00766f.

Acknowledgements

This work is supported by the DST-SERB-Core Research Grant (File No.: CRG/2022/008594), the DST-INSPIRE Faculty Award (File No.: DST/INSPIRE/04/2016/000735), the DST-SERB Early Career Research Award (File No.: ECR/2018/002635) Schemes by the Department of Science and Technology, Government of India and the CPRI-RSoP scheme of Ministry of Power (MoP) (File No.: RSoP/21-26/GD/17), Government of India. Furthermore, the authors acknowledge the use of the PARAM UTKARSH HPC facility of the National Supercomputing Mission project and thank the Centre for Cyber Physical Systems (CCPS), National Institute of Technology Karnataka, Surathkal for the financial support for the same. In addition, PKK acknowledges the financial support of Petronet LNG Limited under its CSR sanction head for upgrading the facilities in the Fuel Cell Research Center at NITK Surathkal.

References

1. T. Shimoaka, C. Wakai, T. Sakabe, S. Yamazaki and T. Hasegawa, Phys. Chem. Chem. Phys., 2015, 17, 8843–8849.
2. K. Zhang, W. Yu, X. Ge, L. Wu and T. Xu, J. Membr. Sci., 2022, 661, 120922.
3. Y. Wang, D. Y. C. Leung, J. Xuan and H. Wang, Renewable Sustainable Energy Rev., 2016, 65, 961–977.
4. V. Kumar, P. K. Koorata, U. Shinde, P. Padavu and S. C. George, Polym. Degrad. Stab., 2022, 205, 110151.
5. Y. Wang, K. S. Chen, J. Mishler, S. C. Cho and X. C. Adroher, Appl. Energy, 2011, 88, 981–1007.
6. X. Ren, E. Gobrogge and F. L. Beyer, J. Membr. Sci., 2021, 637, 119645.
7. P. K. Koorata, J. Electrochem. Soc., 2022, 169, 104505.
8. U. Shinde and P. K. Koorata, J. Power Sources, 2022, 527, 231179.
9. D. E. Curtin, R. D. Lousenberg, T. J. Henry, P. C. Tangeman and M. E. Tisack, J. Power Sources, 2004, 131, 41–48.
10. P. K. Koorata and S. D. Bhat, Int. J. Hydrogen Energy, 2021, 46, 5570–5579.
11. T. Taner, J. Therm. Eng., 2019, 5, 456–468.
12. P. Chaturvedi, N. K. Moehring, T. Knight, R. Shah, I. Vlassiouk and P. R. Kidambi, Mater. Adv., 2023, 4, 3473–3481.
13. K. K. Poornesh and C. Cho, Fuel Cells, 2015, 15, 196–203.
14. P. W. Majsztrik, A. B. Bocarsly and J. B. Benziger, Macromolecules, 2008, 41, 9849–9862.
15. R. Ding, S. Zhang, Y. Chen, Z. Rui, K. Hua, Y. Wu, X. Li, X. Duan, X. Wang, J. Li and J. Liu, Energy AI, 2022, 9, 100170.
16. Y. Wang, B. Seo, B. Wang, N. Zamel, K. Jiao and X. C. Adroher, Energy AI, 2020, 1, 100014.
17. K. I. Ozoemena, RSC Adv., 2016, 6, 89523–89550.
18. J. Peron, A. Mani, X. Zhao, D. Edwards, M. Adachi, T. Soboleva, Z. Shi, Z. Xie, T. Navessin and S. Holdcroft, J. Membr. Sci., 2010, 356, 44–51.
19. A. Parekh, Front. Energy Res., 2022, 10, 1–13.
20. Z. Yang, H. Peng, W. Wang and T. Liu, J. Appl. Polym. Sci., 2010, 116, 2658–2667.
21. P. Padavu, P. K. Koorata, S. Kattimani and D. N. Gaonkar, Int. J. Hydrogen Energy, 2024, 84, 435–446.
22. J. Savage, Y.-L. S. Tse and G. A. Voth, J. Phys. Chem. C, 2014, 118, 17436–17445.
23. A. Ibrahim, O. Hossain, J. Chaggar, R. Steinberger-Wilckens and A. El-Kharouf, Int. J. Hydrogen Energy, 2020, 45, 5526–5534.
24. R. Yadav, A. Subhash, N. Chemmenchery and B. Kandasubramanian, Ind. Eng. Chem. Res., 2018, 57, 9333–9350.
25. P. P. Singh and R. Ranganathan, ACS Omega, 2024, 9, 9063–9075.
26. Y. Hu, J. Li and S. Wang, Int. J. Hydrogen Energy, 2023, 48, 31734–31744.
27. A. Kumar, K. Sharma and A. R. Dixit, Mol. Simul., 2020, 46, 136–154.
28. L. Stobinski, B. Lesiak, A. Malolepszy, M. Mazurkiewicz, B. Mierzwa, J. Zemek, P. Jiricek and I. Bieloshapka, J. Electron Spectrosc. Relat. Phenom., 2014, 195, 145–154.
29. N. M. S. Hidayah, W. W. Liu, C. W. Lai, N. Z. Noriman, C. S. Khe, U. Hashim and H. C. Lee, in AIP Conference Proceedings, American Institute of Physics Inc., 2017, vol. 1892, p. 150002.
30. W. Li, J. Liu and C. Yan, Carbon, 2013, 55, 313–320.
31. R. K. Joshi, S. Alwarappan, M. Yoshimura, V. Sahajwalla and Y. Nishina, Appl. Mater. Today, 2015, 1, 1–12.
32. S. Miao, H. Zhang, X. Li and Y. Wu, Int. J. Hydrogen Energy, 2016, 41, 331–341.
33. J. Song, X. Wang and C. T. Chang, J. Nanomater., 2014, 2014, 276143.
34. A. Lerf, A. Buchsteiner, J. Pieper, S. Schöttl, I. Dekany, T. Szabo and H. P. Boehm, J. Phys. Chem. Solids, 2006, 67, 1106–1110.
35. N. V. Medhekar, A. Ramasubramaniam, R. S. Ruoff and V. B. Shenoy, ACS Nano, 2010, 4, 2300–2306.
36. T. Cheng, M. Feng, Y. Huang and X. Liu, Ionics, 2017, 23, 2143–2152.
37. Z. Yang, Y. Sun and F. Ma, Appl. Surf. Sci., 2020, 529, 147075.
38. A. B. Peressut, M. Di Virgilio, A. Bombino, S. Latorrata, E. Muurinen, R. L. Keiski and G. Dotelli, Molecules, 2022, 27, 1507.
39. R. S. Rajaura, S. Srivastava, V. Sharma, P. K. Sharma, C. Lal, M. Singh, H. S. Palsania and Y. K. Vijay, Int. J. Hydrogen Energy, 2016, 41, 9454–9461.
40. M. R. Karim, M. S. Islam, K. Hatakeyama, M. Nakamura, R. Ohtani, M. Koinuma and S. Hayami, J. Phys. Chem. C, 2016, 120, 21976–21982.
41. C. Vallés, F. Beckert, L. Burk, R. Mülhaupt, R. J. Young and I. A. Kinloch, J. Polym. Sci., Part B:Polym. Phys., 2016, 54, 281–291.
42. A. M. Dimiev and J. M. Tour, ACS Nano, 2014, 8, 3060–3068.
43. J. Chen, X. Zhang, X. Zheng, C. Liu, X. Cui and W. Zheng, Mater. Chem. Phys., 2013, 139, 8–11.
44. H. Hou, X. Hu, X. Liu, W. Hu, R. Meng and L. Li, Ionics, 2015, 21, 1919–1923.
45. S. Zhou, S. Kim, E. Di Gennaro, Y. Hu, C. Gong, X. Lu, C. Berger, W. De Heer, E. Riedo, Y. J. Chabal, C. Aruta and A. Bongiorno, Adv. Mater. Interfaces, 2014, 1, 1300106.
46. J. A. L. Willcox and H. J. Kim, J. Phys. Chem. C, 2017, 121, 23659–23668.
47. L. Wang, J. Kang, J. Do Nam, J. Suhr, A. K. Prasad and S. G. Advani, ECS Electrochem. Lett., 2015, 4, F1–F4.
48. B. G. Choi, Y. S. Huh, Y. C. Park, D. H. Jung, W. H. Hong and H. Park, Carbon, 2012, 50, 5395–5402.
49. J. E. Cha, S. Jang, D. J. Seo, J. Hwang, M. H. Seo, Y. W. Choi and W. B. Kim, Chem. Eng. J., 2023, 454, 140091.
50. T. Bayer, S. R. Bishop, M. Nishihara, K. Sasaki and S. M. Lyth, J. Power Sources, 2014, 272, 239–247.
51. Q. Zheng, Q. Xue, K. Yan, L. Hao, Q. Li and X. Gao, J. Phys. Chem. C, 2007, 111, 4628–4635.
52. T. K. Maiti, J. Singh, S. K. Maiti, J. Majhi, A. Ahuja, M. Singh, A. Bandyopadhyay, G. Manik and S. Chattopadhyay, Eur. Polym. J., 2022, 174, 111345.
53. G. Kritikos, R. Pant, S. Sengupta, K. Karatasos, A. Venkatnathan and A. V. Lyulin, J. Phys. Chem. C, 2018, 122, 22864–22875.
54. R. Tanaka, T. Mabuchi, Y. Zang, B. Hinds and T. Tokumasu, ECS Trans., 2021, 104, 309–316.
55. R. Devanathan, D. Chase-Woods, Y. Shin and D. W. Gotthold, Sci. Rep., 2016, 6, 29484.
56. J. Cha, Sci. Rep., 2020, 10, 22014.
57. A. Venkatnathan, R. Devanathan and M. Dupuis, J. Phys. Chem. B, 2007, 111, 7234–7244.
58. K. Karatasos and G. Kritikos, RSC Adv., 2016, 6, 109267–109277.
59. G. Kritikos and K. Karatasos, Mater. Today Commun., 2017, 13, 359–366.
60. M. D. Hanwell, D. E. Curtis, D. C. Lonie, T. Vandermeersch, E. Zurek and G. R. Hutchison, J. Cheminf., 2012, 4, 17.
61. L. Martinez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
62. I. H. Hristov, S. J. Paddison and R. Paul, J. Phys. Chem. B, 2008, 112, 2937–2949.
63. T. Mabuchi and T. Tokumasu, J. Chem. Phys., 2014, 141, 104904.
64. M. K. Petersen, F. Wang, N. P. Blake, H. Metiu and G. A. Voth, J. Phys. Chem. B, 2005, 109, 3727–3730.
65. D. E. Breen, D. H. House and M. J. Wozny, Text. Res. J., 1994, 64, 663–685.
66. R. Devanathan, A. Venkatnathan and M. Dupuis, J. Phys. Chem. B, 2007, 111, 8069–8079.
67. Y.-L. S. Tse, A. M. Herring, K. Kim and G. A. Voth, J. Phys. Chem. C, 2013, 117, 8079–8091.
68. S. Ban, C. Huang, X.-Z. Yuan and H. Wang, J. Phys. Chem. B, 2011, 115, 11352–11358.
69. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
70. S. S. Awulachew and K. N. Nigussa, Chem. Phys. Lett., 2023, 816, 140380.
71. W. G. Hoover, Phys. Rev. A, 1985, 31, 1695–1697.
72. W. G. Hoover, Phys. Rev. A, 1986, 34, 2499–2500.
73. R. Devanathan, A. Venkatnathan and M. Dupuis, J. Phys. Chem. B, 2007, 111, 13006–13013.
74. G. Zhang, G. Yang, S. Li, Q. Shen, H. Wang, Z. Li, Y. Zhou and W. Ye, Membranes, 2021, 11, 695.
75. R. Kumar, C. Xu and K. Scott, RSC Adv., 2012, 2, 8777–8782.
76. A.-Tsung Kuo, W. Shinoda and S. Okazaki, J. Phys. Chem. C, 2016, 120, 25832–25842.

---