OH− and H3O+ Diffusion in Model AEMs and PEMs at Low Hydration: Insights from Ab Initio Molecular Dynamics

Fuel cell-based anion-exchange membranes (AEMs) and proton exchange membranes (PEMs) are considered to have great potential as cost-effective, clean energy conversion devices. However, a fundamental atomistic understanding of the hydroxide and hydronium diffusion mechanisms in the AEM and PEM environment is an ongoing challenge. In this work, we aim to identify the fundamental atomistic steps governing hydroxide and hydronium transport phenomena. The motivation of this work lies in the fact that elucidating the key design differences between the hydroxide and hydronium diffusion mechanisms will play an important role in the discovery and determination of key design principles for the synthesis of new membrane materials with high ion conductivity for use in emerging fuel cell technologies. To this end, ab initio molecular dynamics simulations are presented to explore hydroxide and hydronium ion solvation complexes and diffusion mechanisms in the model AEM and PEM systems at low hydration in confined environments. We find that hydroxide diffusion in AEMs is mostly vehicular, while hydronium diffusion in model PEMs is structural. Furthermore, we find that the region between each pair of cations in AEMs creates a bottleneck for hydroxide diffusion, leading to a suppression of diffusivity, while the anions in PEMs become active participants in the hydronium diffusion, suggesting that the presence of the anions in model PEMs could potentially promote hydronium diffusion.


Description of Systems
A typical AEM or PEM model requires a polymer backbone [1], a cationic/anionic group to ensure charge neutrality [56], and a molecular tether that connects the two via covalent bonds. Hydroxide and hydronium ions serve as charge carriers, respectively in AEMs and PEMs, and these systems are then solvated with water molecules [57]. In previous studies, we explored different GB systems as mimics of different AEM and PEM environments [52][53][54][55] in which the nanoconfined structure in these systems model the layered arrangement recently reported in References [16,57]. Of the various systems studied, we ultimately chose two representative examples that best captured the hydroxide and hydronium ion diffusion mechanisms described above. Each system contains two identical graphane layers aligned in the xy-plane and a number of water molecules (the hydration level) chosen to model the experimentally relevant hydration values from References [16,57]. For the AEM model, the system contains two tetramethylammonium (TMA) cations and two hydroxide ions, while for the PEM model, the system contains two SO − denoted O* 1 and O* 2 throughout the paper. The two cations/anions are attached by (CH 2 ) 2 linkers to fixed points in the GBs but are otherwise free to move in the aqueous solution. The two attachment points define the polymer electrolyte cation/anion spacing in the x and y directions (see Figure 1). As a result, the simulation cell is partitioned into an open region in the center of the cell and constricted regions between the cations/anions, which we refer to as bottleneck regions (BRs) in the discussions of the AEM models. Based on References [16,57], the tunable parameters for the two systems are: (i) the hydration number, λ, chosen to be 4 or 3 for both the AEM and PEM models, respectively, (ii) the polymer electrolyte cation/anion spacing in the x direction, ∆x, as measured between two nitrogen (AEM) or sulfur (PEM) atoms, which is fixed at 10 Å for the two systems, (iii) the polymer electrolyte cation/anion spacing in the y direction, ∆y, as measured between two nitrogen (AEM) or sulfur (PEM) atoms, which is fixed at 8.7 Å and 6.6 Å for the model AEM and PEM systems, respectively, and (iv) the distance between the two carbon sheets, ∆z, fixed at 7.3 Å for the two systems. The distance ∆z is calculated as the distance between the hydrogen atoms on the inner surfaces of the graphane layers and is determined by λ, with a lower bound determined by the height of the cations/anions. These AEM and PEM models were previously referred to as systems a4 and λ3 in References [53,55], respectively. hydration level) chosen to model the experimentally relevant hydration values from References [16,57]. For the AEM model, the system contains two tetramethylammonium (TMA) cations and two hydroxide ions, while for the PEM model, the system contains two SO 3 ─ anions and two hydronium ions. The hydroxide and hydronium ion oxygen cores are denoted O*1 and O*2 throughout the paper. The two cations/anions are attached by (CH2)2 linkers to fixed points in the GBs but are otherwise free to move in the aqueous solution. The two attachment points define the polymer electrolyte cation/anion spacing in the x and y directions (see Figure 1). As a result, the simulation cell is partitioned into an open region in the center of the cell and constricted regions between the cations/anions, which we refer to as bottleneck regions (BRs) in the discussions of the AEM models. Based on References [16,57], the tunable parameters for the two systems are: (i) the hydration number, λ, chosen to be 4 or 3 for both the AEM and PEM models, respectively, (ii) the polymer electrolyte cation/anion spacing in the x direction, Δx, as measured between two nitrogen (AEM) or sulfur (PEM) atoms, which is fixed at 10 Å for the two systems, (iii) the polymer electrolyte cation/anion spacing in the y direction, Δy, as measured between two nitrogen (AEM) or sulfur (PEM) atoms, which is fixed at 8.7 Å and 6.6 Å for the model AEM and PEM systems, respectively, and (iv) the distance between the two carbon sheets, Δz, fixed at 7.3 Å for the two systems. The distance Δz is calculated as the distance between the hydrogen atoms on the inner surfaces of the graphane layers and is determined by , with a lower bound determined by the height of the cations/anions. These AEM and PEM models were previously referred to as systems a4 and 3 in References [53,55], respectively.

Computational Method
AIMD simulations [51] were performed on the aforementioned systems using the CPMD code [58,59]. The two systems were first equilibrated at room temperature using a massive Nosé-Hoover chain thermostat [60], followed by 15-20 ps of canonical (NVT) dynamics, also using a massive Nosé-Hoover chain thermostat, and finally~80 ps of microcanonical (NVE) dynamics. Dispersion forces were accounted for using the Dispersion-Corrected Atomic Core Pseudopotential (DCACP) scheme [61,62] within the Kohn-Sham formulation of Density Functional Theory using the B-LYP exchange-correlation functional [63,64]. The B-LYP+DCACP has been selected for this study, as it has previously been shown to produce satisfactory results for water-acene interactions [65], liquid water [66], and hydronium diffusion in bulk water [38][39][40]67] A detailed description of the construction of the initial structures and of the computational methodology can be found in the Supplementary Information (SI) and in References [52][53][54][55].

Water Distribution
Inspection of the AIMD trajectories reveals that under low hydration conditions, the water distribution is non-uniform in both the AEM and PEM systems. The non-uniformity is pronounced in void areas throughout the system, in which all water molecules in the system are in contact with some part of the "membrane" and inhomogeneous throughout the system. Unlike in bulk solution, in which the water oxygen has an average of a fourfoldtetrahedral coordination pattern [38,39], the non-uniform water distribution results in a first solvation shell that oscillates between zero, one and two for the water oxygens. In order to demonstrate the existence of water-oxygen solvation shells, we present the O w O w radial distribution functions (RDFs) and coordination numbers (CNs) for the two systems in Figure 2. As shown for both cases, the first peak located at~2.8 Å corresponds roughly to that of bulk water. However, as a result of the non-uniform water distribution, the CN values for the first and second solvation shells are lower than those in the bulk [38,39], with values of~1 and~2.7 for the AEM model, and~1.5 and~4 for the PEM model, respectively.

Computational Method
AIMD simulations [51] were performed on the aforementioned systems using the CPMD code [58,59]. The two systems were first equilibrated at room temperature using a massive Nosé-Hoover chain thermostat [60], followed by 15-20 ps of canonical (NVT) dynamics, also using a massive Nosé-Hoover chain thermostat, and finally ~80 ps of microcanonical (NVE) dynamics. Dispersion forces were accounted for using the Dispersion-Corrected Atomic Core Pseudopotential (DCACP) scheme [61,62] within the Kohn-Sham formulation of Density Functional Theory using the B-LYP exchange-correlation functional [63,64]. The B-LYP+DCACP has been selected for this study, as it has previously been shown to produce satisfactory results for water-acene interactions [65], liquid water [66], and hydronium diffusion in bulk water [38][39][40]67] A detailed description of the construction of the initial structures and of the computational methodology can be found in the Supplementary Information (SI) and in References [52][53][54][55].

Water Distribution
Inspection of the AIMD trajectories reveals that under low hydration conditions, the water distribution is non-uniform in both the AEM and PEM systems. The non-uniformity is pronounced in void areas throughout the system, in which all water molecules in the system are in contact with some part of the "membrane" and inhomogeneous throughout the system. Unlike in bulk solution, in which the water oxygen has an average of a fourfold-tetrahedral coordination pattern [38,39], the non-uniform water distribution results in a first solvation shell that oscillates between zero, one and two for the water oxygens In order to demonstrate the existence of water-oxygen solvation shells, we present the OwOw radial distribution functions (RDFs) and coordination numbers (CNs) for the two systems in Figure 2. As shown for both cases, the first peak located at ~2.8 Å corresponds roughly to that of bulk water. However, as a result of the non-uniform water distribution the CN values for the first and second solvation shells are lower than those in the bulk [38,39], with values of ~1 and ~2.7 for the AEM model, and ~1.5 and ~4 for the PEM model respectively.

OH − Solvation Structure and Dynamical Properties
As mentioned in the previous section, the water distribution in these low hydration models is non-uniform. Specifically, for the AEM model, the non-uniformity refers to the formation of separated water clusters in the vicinity of each OH − by roughly 4 Å. As a result, the two hydroxide oxygens in the AEM model may not share the same environment or have the same CN. Therefore, we find it useful to refer to the two hydroxide ions as two different species, O* 1 and O* 2 , and to explore the O*O w RDFs and CNs for the pair (O*), O* 1 , and O* 2 separately (see Figure 3a). The first solvation shell of the hydroxide oxygen is located at 2.7 Å, as was previously shown for bulk solution [38][39][40][41][42][43][44][45][46] and contains three water oxygens (as opposed to four water oxygens as was observed in bulk solution). The CN of the second solvation shell (located at~4 Å in bulk solution [38][39][40][41][42][43][44][45][46]) shows that O* 1 is missing a second solvation shell while O* 2 has a second solvation shell of two water oxygens (see inset of Figure 3a for example). It is well known that the hydroxide ion diffusion mechanism is strongly affected by the solvation structure of OH - [38][39][40][41][42][43][44][45][46]. While in bulk solution, the hydroxide ions typically share a similar solvation patterns, which results in a similar diffusion mechanism in this low-hydration AEM model, the heterogeneity of the environments of each hydroxide suggests that over the time scales of the simulations, each OHmight be governed by different diffusion mechanisms, which we discuss next.
The GB structures have a unique confined geometry that influences the mobility of the waters and hydroxide ions differently along each axis [52]. Therefore, in order to gain a better understanding of the water and hydroxide ion diffusion processes in these structures, we calculate diffusion coefficients along each of the axes separately (see Figure 3b), plot the coordinates of the OHoxygens as a function of time along the xand y-axes (the coordinates in the z-axis are not presented as their contribution to the total diffusion is negligibly), and label proton transfer (PT) events with gray lines (see Figure 3c). As PT events represent a change in the hydroxide oxygen identity, they are used to identify different hydroxide ion diffusion mechanisms. While a sharp jump in the OHoxygens coordinates, caused by a PT event, is associated with structural diffusion [41][42][43][44][45][46] usually known as "Grotthuss diffusion", a continuous change in the OHoxygen coordinates is associate with a vehicular diffusion [52]. As shown in Figure 3b, the hydroxide ion diffuses along the x-axis with a diffusion coefficient of 0.343 Å 2 /ps. However, according to Figure 3c, only O* 2 is diffusive while O* 1 is non-diffusive. The lack of PT events along the trajectory suggests that the diffusion of O* 2 is mainly vehicular.
Based on the results presented in Section 4.2.1, we find that the OHions diffuse only when they have a second solvation shell of at least one water oxygen. The second solvation shell provides the necessary hydration that requires for the OHto shift the competition between its hydrophobic repulsion from the cation and its electrostatic attraction to the cation, which results in vehicular diffusion [38][39][40][41][42][43][44][45][46][47][48]. Hence, we find that the presence of a second solvation shell promotes vehicular diffusion in the case of O* 2, while the absence of a second solvation shell suppresses OHdiffusion as seen for O* 1 .
fusion mechanism is strongly affected by the solvation structure of OH - [38][39][40][41][42][43][44][45][46]. While in bulk solution, the hydroxide ions typically share a similar solvation patterns, which results in a similar diffusion mechanism in this low-hydration AEM model, the heterogeneity of the environments of each hydroxide suggests that over the time scales of the simulations, each OHmight be governed by different diffusion mechanisms, which we discuss next.

OH -Diffusion Mechanism
For the AEM model under low hydration conditions, vehicular diffusion was found to be the dominant diffusion mechanism [53]. The preference for vehicular over structural diffusion relates to the non-uniform water distribution, which results in an insufficient number of water molecules in the vicinity of the hydroxide ions to promote PT events [53]. Additionally, we find that the hydroxide diffusivity is affected by its location in the cell: the hydroxide ion diffuses relatively freely when located at the center of the cell while its diffusion is restricted when it is located in BRs where the diffusion path between a pair of cations is relatively narrow (see Figure 1). In Figure 4, we propose an idealized vehicular diffusion process and corresponding solvation patterns, while emphasizing specific diffusion requirements for each region. For this purpose, we use O* 2 of the model AEM along with five water molecules from the first and second solvation shells. Initially, the hydroxide ion is in a stable, "resting", threefold structure near a cation [52], with two water molecules in the second solvation shell (Figure 4a). Once most of the water molecules are located in the center of the cell, the hydroxide ion drifts towards them via vehicular diffusion, forming a fourfold planar structure in the center of the cell with one water molecule in the second solvation shell (Figure 4b). The hydroxide ion and the five water molecules continue to diffuse via vehicular diffusion towards the nearby cation (Figure 4c), until the hydroxide ion forms, once again, the stable, "resting" threefold structure near the next cation (Figure 4d). To promote hydroxide ion diffusion into the BR, a water molecule must be located near the BR entrance. Once the hydroxide ion is located in the BR, the hydroxide ion forms a threefold tetrahedral structure as part of a well-ordered stable OH − (H 2 O) 4 complex [52], in which three water molecules are part of a threefold tetrahedral structure, where a water molecule is located below the hydroxide ion in the BR, and one water molecules is in the second solvation shell (Figure 4e). Next, the hydroxide ion and the five water molecules diffuse vehicularly through the BR (Figure 4f). We find that the existence of a water molecule above/below the hydroxide ion, as part of the tetrahedral structure, is essential in establishing OHdiffusion through a BR. Additionally, in order to draw the hydroxide ion along the BR, a water molecule is required to be located on the other side of the BR. Once most of the first solvation shell water oxygens diffuse into the next open region, the hydroxide ion is drawn towards the other side of the BR and forms a stable, "resting", threefold structure near the next cation (Figure 4g). The mechanistic cycle is complete once the threefold hydroxide ion structure drifts towards the center of the cell and reforms into a fourfold planar structure (Figure 4h). that the existence of a water molecule above/below the hydroxide ion, as part of the tetrahedral structure, is essential in establishing OHdiffusion through a BR. Additionally, in order to draw the hydroxide ion along the BR, a water molecule is required to be located on the other side of the BR. Once most of the first solvation shell water oxygens diffuse into the next open region, the hydroxide ion is drawn towards the other side of the BR and forms a stable, "resting", threefold structure near the next cation (Figure 4g). The mechanistic cycle is complete once the threefold hydroxide ion structure drifts towards the center of the cell and reforms into a fourfold planar structure ( Figure 4h). As a vehicular diffusion process requires synchronicity between both of the hydroxide ions and water molecules, the similarity in the diffusivities of both species is expected.

H 3 O + Solvation Structure and Dynamical Properties
Next, we turn to explore the hydronium ion solvation structure in the model PEM. We start with the O*O RDF and CNs presented in Figure 5 (O* represents the hydronium oxygens and O represents SO 3 − and water oxygens). As shown, the first solvation shell is located at 2.6 Å, and the CN values are 3.1 and 7.5 for the first and second solvation shells, respectively. Population probabilities for the hydronium ion solvation complexes show that the most common structure is 3A + 0D with 90% (see SI for hydrogen bond (HB) criteria). Additionally, we find that the oxygens taking part in the first solvation shell of As a vehicular diffusion process requires synchronicity between both of the hydroxide ions and water molecules, the similarity in the diffusivities of both species is expected.

H 3 O + Solvation Structure and Dynamical Properties
Next, we turn to explore the hydronium ion solvation structure in the model PEM. We start with the O*O RDF and CNs presented in Figure 5 (O* represents the hydronium oxygens and O represents SO − 3 and water oxygens). As shown, the first solvation shell is located at 2.6Å, and the CN values are 3.1 and 7.5 for the first and second solvation shells, respectively. Population probabilities for the hydronium ion solvation complexes show that the most common structure is 3A + 0D with 90% (see SI for hydrogen bond (HB) criteria). Additionally, we find that the oxygens taking part in the first solvation shell of the hydronium ion are both SO − 3 and water oxygens (see inset of Figure 5). Excluding oxygens from the first solvation shell and the SO − 3 oxygens, the number of water oxygens in the second solvation shell of the hydronium ions is 1.5. This suggests the hydronium ions in the PEM model system are missing a complete second solvation shell as a result of the non-uniform water distribution present in the system. In order to enhance the analysis of the hydronium ion diffusion mechanism in this confined environment, we calculate water and hydronium diffusion coefficients for each of the axes (see Figure 5b). Moreover, in Figure 5c we present the coordinates of the hydronium oxygens along the trajectory. We find that the hydronium ion diffusion occurs along the x-and y-axes, with diffusion coefficients of 0.106 Å 2 /ps and 0.095 Å 2 /ps, respectively. According to Figure 5c, both ions become diffusive (at ~50 ps) after an initial "resting" period. Unlike in the model AEM, the water molecules in the model PEM were found to be non-diffusive, suggesting that water mobility is not necessarily required for In order to gain a better understanding of the conditions that enable the diffusion of the two hydronium ions, we plot in Figure 6 the SO* RDF and CNs. As shown, the first  Figure (a,c) are reprinted with permission from Reference [55]. Copyright 2020 Royal Society of Chemistry.
In order to enhance the analysis of the hydronium ion diffusion mechanism in this confined environment, we calculate water and hydronium diffusion coefficients for each of the axes (see Figure 5b). Moreover, in Figure 5c we present the coordinates of the hydronium oxygens along the trajectory. We find that the hydronium ion diffusion occurs along the xand y-axes, with diffusion coefficients of 0.106 Å 2 /ps and 0.095 Å 2 /ps, respectively. According to Figure 5c, both ions become diffusive (at~50 ps) after an initial "resting" period. Unlike in the model AEM, the water molecules in the model PEM were found to be non-diffusive, suggesting that water mobility is not necessarily required for H 3 O + diffusion.
In order to gain a better understanding of the conditions that enable the diffusion of the two hydronium ions, we plot in Figure 6 the SO* RDF and CNs. As shown, the first peak, located at~1.7 Å, corresponds to SO 3 H, in which an H 3 O + has transferred a proton to an SO − 3 , while the second peak is located at~3.7 Å. The SO* CN values for the first and second solvation shells are 0.7 and 0.8, respectively. This suggests that, as a result of the following reaction: SO − 3 +H 3 O + ↔ SO 3 H + H 2 O , a neutral species exists in the system. To verify this, we calculate the percentage of time that the hydronium ions spent as SO − 3 and SO 3 H, which shows that the hydronium ion appears as SO 3 H for 48.78% of the simulation time (see inset of Figure 6). This leads us to conclude that the protonation state of the sulfonate group is a significant factor in the hydronium ion diffusion mechanism, as the reaction: , is a key component in the diffusion process. In order to dive a bit deeper into the conditions that enable this reaction, we plot in Figure 7a the OnextO RDF and CNs, in which Onext, is the closest water or hydronium oxygen to the SO 3 − oxygens, and O represents all water and hydronium oxygens. As shown, the first peak is located at 2.7 Å , and the CN values for the first and second solvation shells are 1.3 and 5.1, respectively. The CN values found suggest that before a PT occurs between the anion and a nascent water molecule, the water must have a first solvation shell of one water oxygen and an incomplete second solvation shell. This leads us to conclude that the protonation state of the sulfonate group is a significant factor in the hydronium ion diffusion mechanism, as the reaction: , is a key component in the diffusion process. In order to dive a bit deeper into the conditions that enable this reaction, we plot in Figure 7a the O next O RDF and CNs, in which O next , is the closest water or hydronium oxygen to the SO − 3 oxygens, and O represents all water and hydronium oxygens. As shown, the first peak is located at 2.7 Å, and the CN values for the first and second solvation shells are 1.3 and 5.1, respectively. The CN values found suggest that before a PT occurs between the anion and a nascent water molecule, the water must have a first solvation shell of one water oxygen and an incomplete second solvation shell.
and O represents all water and hydronium oxygens. As shown, the first peak is loca 2.7 Å , and the CN values for the first and second solvation shells are 1.3 and 5.1, r tively. The CN values found suggest that before a PT occurs between the anion nascent water molecule, the water must have a first solvation shell of one water o and an incomplete second solvation shell.  In support of this claim, we define a displacement coordinate, δ = |R O a H − R O w H |, where R O a H and R O b H are the distances between a shared proton of SO 3 H and the nearest water oxygen (i.e., O next ). Values of δ > 0.5 are considered to be inactive complexes with respect to PT, while values of δ < 0.1 are considered to be "active" and are associated with PT events [41,42,52,69,70]. In Figure 7b, we present the O next O RDF and CNs for δ < 0.1 and δ > 0.5, where, O next represents the first neighbor oxygen to the SO 3 H oxygen, and O represents water and hydronium oxygens. We find that for δ > 0.5, the peak is located at 2.8 Å with a CN value of 1.54, while for δ < 0.1, the peak is located at 2.7 Å with a CN value of 1.14. This suggests that in order for the reaction to occur, O next is required to have a CN value of~1.
Combining the results presented in Section 4.3.1, we conclude that the factors required to increase the hydronium reactivity in the PEM model relate to the non-uniformity of the water distribution along the membrane, which gives rise to a high probability of obtaining a CN value of~1 for O next , an incomplete second solvation shell for the hydronium ions, and fewer water molecules in the vicinity of the anions oxygens (i.e., SO − 3 ).

H 3 O + Diffusion Mechanism
Based on the results presented in Section 4.3.1 and combined with inspection of configurations from the AIMD trajectory, we propose, in Figure 8, an idealized diffusion mechanism for hydronium ions in the PEM model under idealized hydration conditions (λ = 3) [55]. First, the hydronium ion is located in the center of the cell, solvated by three water molecules (Figure 8a). Next, a PT occurs from the hydronium ion to a neighboring water molecule (Figure 8b). A HB is formed between the nascent hydronium ion oxygen and the anion, SO − 3 , while the hydronium has only one water oxygen in its first solvation shell (Figure 8c). Finally, a PT occurs between the hydronium ion and the anion (i.e, SO − 3 +H 3 O + ), resulting in SO 3 H + H 2 O (Figure 8d). This procedure cycles back to the initial condition and restarts, meaning that the next PT will occur once SO 3 H donates its hydrogen to a neighboring water molecule with a first solvation shell consisting of only one water oxygen (Figure 8f), which will result in SO   (Figure 8d). This procedure cycles back to th tial condition and restarts, meaning that the next PT will occur once SO 3 H dona hydrogen to a neighboring water molecule with a first solvation shell consisting o one water oxygen (Figure 8f), which will result in SO 3 − + H 3 O + (see details in Figu 8h).

Discussions and Conclusions
In the present work, we used AIMD simulations to gain an in-depth atomistic perspective of two idealized AEM and PEM models in confined geometries under low hydration conditions (λ = 3, and 4). We found that for both systems, the water distribution within the simulation cell is not uniform. However, the effect of this unique water distribution on the hydroxide and hydronium diffusion mechanisms is fundamentally different. For the AEM model, as long as the hydroxide ions have both first and second solvation shells, the diffusion mechanism is mostly vehicular [53]. However, for the PEM model, hydronium ion diffusion is structural rather than vehicular, with the participations of the anions according to the reaction: SO − 3 +H 3 O + ↔ SO 3 H + H 2 O [55]. Comparing the water diffusion, we find that for the AEM model, the diffusion coefficients of the hydroxide ions and the water molecules are similar, as vehicular diffusion requires synchronized diffusion of both species. However, for the PEM model, the water molecules were found to diffuse much slower compared to the hydronium ions, suggesting that water mobility is not necessarily required in the hydronium ion structural diffusion process.
Most importantly, we find that the differences between the AEM and PEM models lie in the essence of the membrane materials. The region between each pair of cations in the AEM system was found to create a BR for hydroxide diffusion, as only specific solvation structures are diffusive, leading to a suppression of hydroxide ion mobility [53]. However, we find that this obstacle is not present in the PEM model, as the anions in the membrane play an active role in the hydronium diffusion mechanism via the reaction SO − we have been able to provide atomistic insight, gaining a fundamental understanding of the uniqueness hydroxide and hydronium ion solvation patterns and diffusion mechanisms in this hitherto unstudied regime. We believe that elucidating the key design principles underlying the atomistic differences between the hydroxide and hydronium ion diffusion mechanisms in model AEMs and PEMs can be a first step toward the discovery and characterization of new, stable, highly conductive membrane materials for use in emerging fuel cell technologies.