A Steady-State Monte Carlo Study on the Effect of Structural and Operating Parameters on Liquid Water Distribution within the Microporous Layers and the Catalyst Layers of PEM Fuel Cells

The distribution of liquid water agglomerations within a catalyst layer and a microporous layer of a polymer electrolyte membrane fuel cell employing a Monte Carlo model are simulated. The simulations are based on real material structures and locally resolved operating conditions. The 3D material structures are obtained from focused ion beam tomography capable of nm-scale imaging. The local temperature and the relative humidity proﬁles are provided by a computational ﬂuid dynamics simulation, based on a single channel fuel cell model taking into account all relevant transport and electrochemical processes. The results of this Monte Carlo study conﬁrm the strong dependency of the water saturation in the catalyst layer and microporous layer on the structural parameters and operating conditions, such as: wettability, pore size and local temperature and relative humidity. MC simulation model are combined in order to address the existing challenges in the experiments and the simulations on nm-scale. An MC model, which has been successfully developed to simulate the liquid water distribution in the GDL substrate structure, 2,12 is further developed to be capable of simulations in the CL structures.

Liquid water agglomerations within the gas diffusion layers (GDL) of polymer electrolyte membrane fuel cells (PEMFC) can obstruct the diffusion of the reactant gases from the gas flow channels to the catalyst layer (CL) and therefore, slow down the electrochemical reaction and reduce the overall cell performance. 1 On the other hand, an adequate humidification of the membrane and the ionomer component of the CL is required to ensure the rapid proton transport to the reaction sites at the cathode side. In the recent years, a considerable experimental and simulation progress has been made on understanding the formation of liquid water agglomerations in the GDL fiber substrate, which has helped to decode the mass transport processes of reactants and products. 2,3 However, this knowledge about the microporous layers (MPL) of the GDLs and the CL seems still not to be as expanded and deserves further attention.
In the recent past advanced X-ray and neutron tomography techniques have been successfully employed to elucidate the properties of liquid water in the GDL substrates and fuel cell flow fields. 4,5 However, the nm-scale pore sizes relevant in the MPL and CL, and the interference of the catalyst layer platinum content on the image acquisition constitute fundamental challenges for the aforementioned techniques in these structures. Additionally, the high-resolution tomography techniques capable of reconstructing the nano-structures of these materials are based on electron microscopy and therefore, are not able to image the water distribution in the MPL and CL of an operating cell. Thus, the aim of this work is to employ simulation techniques in order to model the dispersion of liquid water agglomerations in the real MPL and CL structures. Monte Carlo (MC) is a numerical technique, which was used for the first simulation of a liquid by Metropolis et al. [6][7][8] This method, which is based on random sampling of the processes, has been successfully employed to solve different statistical mechanical problems over years. [9][10][11] In this study, focused ion beam -scanning electron microscope (FIB-SEM) nano-tomoghraphy, a computational fluid dynamics (CFD) study and a refined MC simulation model are combined in order to address the existing challenges in the experiments and the simulations on nm-scale. An MC model, which has been successfully developed to simulate the liquid water distribution in the GDL substrate structure, 2,12 is further developed to be capable of simulations in the CL structures.

The MC Simulation Inputs
The inputs of the MC simulations are: the real 3D material structure, the structural and surface characteristics and the boundary conditions. These inputs are provided by three parallel studies, which are explained below.
Tomographic reconstructions of CL and MPL.-For the CL structure a reconstructed dataset based on a FIB-SEM tomography of a Gore PRIMEA A510.1 M710.18 C510.4 membrane electrode assembly (MEA) was taken from Thiele et al. 13 While the CL was taken from previous work, the MPL was reconstructed for this study. The MPL structure of a SGL GDL 25BC was reconstructed via FIB-SEM tomography with pore space filling with ZnO using atomic layer deposition. To enhance the contrast between pores and MPL and facilitate the 3D reconstruction process, the MPL was coated with a ZnO film of 100 nm thickness by allowing diethyl zinc and water to react in a cyclic manner in a vertical-flow, hot-wall reactor (OpAL, manufactured by Oxford Instruments) at 50 • C. 14 The atomic layer deposition process employed allows filling the void space within the porous material without changing the material structure. The FIB-SEM tomography was conducted with a Zeiss Auriga 60 dual beam comprising 245 FIB cuts and SEM images. Each FIB section was conducted at 30 kV accelerating voltage and 20 pA beam current with a cutting distance of 9 nm. SEM images were acquired at 5 kV with a pixel size of 3 nm using an in-lens detector. The total reconstructed MPL volume was 3.7 × 2.2 × 2.2 μm 3 . Images were aligned with MATLAB (Mathworks, Inc.) while image enhancement and segmentation (Otsu threshold and minor manual corrections) were conducted with Fiji, a version of ImageJ. 15 These steps are described in detail in Vierrath et al. 16 The parameters for atomic layer deposition, FIB-SEM and the post-processing procedures were equal to those presented in that work. The total reconstructed MPL volume was 3.7 × 2.2 × 2.2 μm 3 based on 245 FIB cuts and SEM images.
Both datasets for MPL and CL were then down sampled to 45 nm and 40.8 nm voxel size, respectively, in order to be used as the input for the 3D grid based MC simulations employing cubic voxels of the same sizes. The choice of these voxel sizes was a compromise between accuracy of the structure and computational effort. Additionally, in order to employ the smaller voxel sizes provided by tomography the molecular effects would have to be taken into account, which is not the case for the MC model used here. This MC model, however, includes water movement based on surface energy minimization as well as condensation and evaporation based on chemical potential calculations, which is sufficient to cover the relevant effects in the 40 nm voxel size class.
Structure and surface characteristics.-Porosities for the MPL (81%) and the CL (58%) were calculated from the reconstructed structure datasets. The inner contact angle of the GDL material was assessed from an inverse gas chromatography experiment. The details of this method are described by Seidenberger et al. 2 The contact angle to water was found to be 92 • for the complete GDL including the substrate and the MPL. As usually for GDLs featuring a MPL, the MPL inner surface constitutes the GDL inner surface to a large extent, the contact angle of the MPL should not differ significantly from the overall value. Furthermore, no major differences were obtained when measuring fiber substrate-only material (96 • for the SGL GDL 25BA). Consequently, a value of 92 • was used for the MPL structure.
The inner surface of the catalyst layer consists of the ionomer, catalyst support and the catalyst. Therefore, an accurate contact angle could be only estimated by taking into account the contribution of each material on the surface. According to Hiesgen et al., the CL of a pristine Nafion-based MEA has average ionomer surface coverage between 30-60%. 17 Thus, the surface coverage of Nafion for the commercial MEA was considered to be 50% in the simulations. The contact angle of ionomer material, e.g. Nafion, is difficult to measure, since it is dependent on the water content of the material and changes considerably over time. 18 Therefore, the simulations were performed using three different measured contact angles: 25 • , 46 • and 80 • . The first one is the receding contact angle measured by Goswami et al., 19 and the last two are the receding and advancing contact angles, respectively, obtained from a surface contact angle measurement on Nafion 115, which was hydrated for 20 hours before the experiment. The experiments were performed using a Krüss DSA 10 goniometer. Due to the small platinum loading, the overall share of platinum on the inner surface is quite low and furthermore, the Pt particles, which are usually on the scale of a few nanometers, cannot be easily represented in the voxel size of 40-50 nm employed in the simulations. The contact angle of pure catalyst support is measured and used for the other 50% of the surface, which has a value of 82 • . Finally, these contact angle values with the respective surface coverages are randomly distributed on the inner surface of the CL structure.

Operating conditions.-A CFD study with the ANSYS Fluent
Fuel Cell Module 20 based on a single channel fuel cell model was performed for a current density of 1.0 A cm −2 at 70 • C gas inlet temperature, atmospheric pressure and 100% relative humidity at the inlet. The results of this study provided a complete set of local data on relative humidity (RH), temperature (T) and pressure in the relevant regions of the GDL substrate, MPL and the CL. The CFD results are illustrated in Figure 1 also sketching a cross-section of the cathode channel geometry used. More details about the CFD methodology employed can be found in the paper of Enz et al. 21 For the purpose of the MC study, the RH and T values are set to be constant in the y (along channel) direction over the relatively short length considered. However, in the RH profile, moving from the area under rib to under the channel, a remarkable gradient can be observed. This is reflected by defining six different sets of thermodynamic boundary conditions corresponding to the color-coding in Figure 1a for the MPL and CL areas (labeled as #1 to #6). In contrast to that, temperature variation is rather remarkable in the z-direction from CL toward the substrate.
Consequently, six different RH values from the regions marked and the corresponding T values are used to observe the effect of these local thermodynamic conditions on the water agglomerations in the CL/MPL pore structures. These boundary conditions are required to calculate the chemical potentials of water in liquid and vapor phases locally, in order to decide for the condensation and evaporation processes. The employed six boundary condition pairs are shown in Table I.

MC Simulation Model and Results
The voxel-based MC algorithm includes a 3D grid, in which the voxel size can be adjusted appropriately according to the resolution of the tomography dataset. Every voxel can be filled with different types of material, such as graphite, PTFE, ionomer, catalyst support, water and air. If required, for every direction of the 3D simulation volume, periodic boundary conditions can be applied. The liquid water is moved within the open pores of the structure in favor of surface energy minimization. Additionally, the water phase can be changed between liquid and gas, based on chemical potential calculations. Every voxel contributes in surface energy calculations, based on the wetting property of the material assigned to it toward water. However, in these  calculations for each water voxel, only the six direct neighbors are taken into account. In other words, the sum of the surface energies between the moving water voxel and the six direct neighbors at the initial position, and the empty voxel (filled with air) and its own neighbors at the final position is calculated, and then will be subtracted from that of when these two voxels swap their places. If the surface energy decreases by this move, it will take place, and if the surface energy increases, it takes place with a probability inversely proportional to the energy increase. More details about the basic characteristics of the algorithm can be found in Seidenberger et al. 2 In this model, steady-state is reached when: firstly, inner hydration of the structure reaches a stable value, and secondly, the condensation and evaporation events per iteration converge.
The simulation results for the two extreme boundary condition pairs, one and six (cf. Table I) within the MPL and CL structures are presented below. Due to the high computational effort required, it was not possible to simulate a statistically relevant number of structure cutouts. However, to somehow check the reproducibility of the results, the simulations are performed at least in two different MPL structure cut-outs. Since the MC simulation box corresponds to approximately one mesh element of the CFD results, constant temperature and relative humidity values are employed within the whole MPL and CL structures. In Figure 2, representative layers of 3D FIB-SEM image stacks of the MPL and CL, which are used in simulations, are shown. In this figure, which shows the structures in the x-y plane, the material and void are illustrated in white and black, respectively. For the simulations in the MPL, a voxel size of 45 nm is employed, and the simulation volume consists of 35 × 35 × 49 voxels. However, for the  CL the simulation volume consists of 40 × 40 × 41 voxels, with the voxel size of 40.8 nm.
The 3D water distribution results in the MPL for under the rib and channel regions are shown in Figure 3. Considering the fact that the temperatures are almost constant, it is observed that with decreasing relative humidity from 105.8% (p partial = 36.4 kPa and T = 345.435K) to 97.0% (p partial = 33.4 kPa and T = 345.501 K) a considerable change in the water amount in the MPL is obtained. Figure 4 illustrates simulation results on how reducing the RH along the x-axis causes a decrease in the water content of the MPL for the two MPL structure cut-outs. Although the 3D results can perfectly illustrate the qualitative allocation of the liquid water within the pores with respect to the thermodynamic boundary conditions and surface wetting properties, no accurate information on pore sizes and their degree of filling can be extracted directly from these results. Therefore, a quantitative pore analysis has been performed, in which the pore size distribution of the structure and the relation of the pore size and the degree of filling with liquid water are studied. Figure 5 presents the result of this pore analysis for the MPL. In this figure, the same two MPL cut-outs as shown in Figure 4 have been investigated, which are titled as region 1 and 2. It is observed that the pore filling degree increases with the pore size in the MPL. For example, although there are only few pores larger than 0.5 μm in region 1, they are completely filled with water for under the rib boundary conditions. Whereas, the pores with the size of 0.1 μm are almost half filled, although the number of these pores in the structure is more than the larger ones by almost a factor   2). The left axis presents the pore filling degree after performing the MC simulations, and the right axis presents the distribution of the pore sizes of the dry structure. of six. The results under the channel show an upward trend as well, however, the total amount of water in the structure is much less.
The simulation results for the CL are categorized for each contact angle of the ionomer material. These results are illustrated in Figure 6 and Figure 7, where for the two boundary condition extrema the ionomer contact angles of 25 • and 80 • are employed, respectively. By comparing the 3D water distribution results of the CL and the MPL, it is observed that as opposed to the MPL, in the CL water occupies the smaller pores earlier. This could be expected from the different wettability properties and the Young-Laplace equation, according to which the sign of the capillary pressure changes at a contact angle value of 90 • , i.e. from going from hydrophilic to hydrophobic. Thus, for hydrophilic material small pores are filled first, whereas for hydrophobic material, the capillary pressure expelling water from the pores is least for the biggest pores. Furthermore, similar to the results in the MPL, in the CL it is affirmed again that regardless of the exact value of the contact angle, the water amount in the media is strongly affected by the thermodynamic boundary conditions. Reducing the RH by 9% causes a significant drop in the water content from 70% to 40%. However, the water unloading starts from the bigger pores in the CL with an ionomer contact angle of 25 • , while in the material with larger contact angle the larger pores seem to remain filled. In Figure 8, the water saturation of the pores for CL with three different contact angles for the ionomer is presented for all six set of boundary conditions. Although the water content seems to be dramatically dependent on the boundary conditions, it is on a small-scale dependent on the contact angle. By almost doubling the contact angle, the water content decreases only by maximum 2%. However, the distribution of water changes from occupying rather smaller pores to larger ones. As the pore analysis results in Figure 9 show, in the CL for all the ionomer contact angles, which are in the hydrophilic range, the pore filling degree decreases with increasing pore size. However, despite keeping the downward trend, by increasing the ionomer contact angle    the filling degree of the smaller pores (<15 μm) decreases, whereas that of the larger pores increases.

Conclusions
Steady-state water distribution in the MPL and the CL of PEMFCs was studied employing a combination of imaging techniques with nm-scale resolution, CFD simulations and an enhanced MC model. In the MC simulations, real material structures and locally resolved thermodynamic boundary conditions have been used. The model includes liquid water movement based on surface energy minimization, and condensations and evaporations based on chemical potentials. In the MPL, an overall hydrophobic contact angle is considered for the whole material. However, in the CL the contact angles of the ionomer and the catalyst support are randomly distributed on the surface with defined surface coverages. The results show the dependency of the water content and distribution in each structure on parameters such as: wettability, local temperature and relative humidity. Since the CFD results showed a considerable change in the RH values from under the rib to under the channel regions, the MC simulations were performed using six different boundary conditions from selected regions within the cell. The 3D water distribution results for two extreme cases with highest and lowest RH values are presented. By decreasing the RH, the water amount in both MPL and the CL was significantly reduced. As expected, the MPL featuring a larger surface contact angle to water showed a hydrophobic behavior, where water tended to fill the larger pores first. In contrast to that, in the CL preferably the smaller pores were occupied. However, in a sensitivity study by changing the ionomer contact angle in the CL, it was observed that although the wettability had only a slight influence on the water content, it notably changed the distribution of water in the structure. By increasing the contact angle of the ionomer, the liquid water showed a tendency to fill the larger pores in the structure. An additional quantitative investigation of the water distribution with respect to the pore size has shown that in the MPL structure the pore filling degree with liquid water increases with the pore size, whereas in the CL it showed an opposite behavior. Additionally, it is observed that changing the ionomer contact angle affected the water dispersion among the pores in the CL.
Overall, this study shows that the MC technique can be applied to simulate the liquid water content also within the small MPL and CL pores in the range of 10 2 to 10 1 nm size employing real structures obtained by sophisticated imaging techniques. Furthermore, the consequences of variations of the surface chemistry of the material resulting in different hydrophobicity can be investigated. Consequently, detailed properties of the water distribution within such material can be investigated, which to the best of our knowledge, so far is beyond the reach of even advanced imaging techniques, such as synchrotron tomography. Thus, in further studies it would be possible to compare e.g. the implications of material variations in terms of structure and surface chemistry on the water content and distribution, which may assist experimental material development.
The difference in the surface wetting properties between CL and MPL could cause difficulties for water to pass through the interface of these two layers. Since the CL and MPL structures have been modeled individually in this study, the behavior of the interfacial region could not be observed. This topic, which requires some developments in the MC algorithm, will be tackled separately.
Although the thermodynamic boundary conditions applied for the MC simulations here stem from CFD calculations employing a full set of experimental cell operating conditions and thus represent a steady-state water distribution related to these conditions, variations of e.g. the current density cannot be directly represented within the MC technique as applied here due to lack of a time axis. Thus, further development of the simulation code will aim for implementing a kinetic Monte-Carlo (KMC) model that allows simulations of the evolution of a certain system in time. Thus, time-dependent variables as e.g. the cell current and the resulting water production can be directly incorporated and varied. This will open new opportunities to support material development for fuel cells by realistic simulations on the nm scale.