Virtual Liquid Water Intrusion in Fuel Cell Gas Diffusion Media

A cluster based full morphology (CFM) model is developed to predict liquid water intrusion in GDLs. The CFM model is used to simulate water intrusion into a dry GDL microstructure obtained from μ-CT. Numerical water saturation distributions are compared to experimental μ-CT reconstructions of the same GDL sample at varying saturation levels. A quantitative validation of the CFM model results is then provided by studying the number of voxels in the image that contain water in both the CFM model results and the μ-CT reconstructions. Results reveal that CFM simulations showed 56–95.7% agreement in the liquid water voxels when compared to the μ-CT simulations for saturations in the range of 29–92.3%. Gas transport simulations are then performed on μ-CT and CFM partially saturated GDLs in order to study the validity of the CFM model images to estimate transport properties of partially saturated GDLs. A maximum error of 20% was observed between the predicted effective diffusivities obtained from CFM simulations and those obtained directly from simulations on the μ-CT data for saturations below 40%. Effective diffusivity predictions from the CFM simulations agree well with in-plane effective diffusivities in literature while the through-plane effective diffusivities were underpredicted by a factor of 2. © The Author(s) 2018. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0921807jes]

Gas diffusion layers (GDLs) are porous media made of either carbon paper or carbon cloth that are often treated with PTFE to introduce hydrophobicity. 1 They are an integral part of proton exchange membrane fuel cells (PEMFCs) as they are responsible for the distribution of reactants from the channels to the catalyst layer (CL), providing an electrical connection between the bipolar plates and CL, removing water and heat produced in the CL and providing mechanical stability to the membrane electrode assembly (MEA). 2 At high current density as well as during cold and wet operation, GDLs are required to provide efficient removal of liquid water in order to sustain reactant transport to the CL and prevent cell shutdown. A better understanding of the effects of the GDL morphology on the liquid water movement, local saturations and gas transport could be crucial in improving the water management in PEMFCs at high current densities.
For GDLs, a typical value for the capillary number (Ca), the ratio of the viscous forces to the surface tension, is of the order of 10 −8 . 3,4 Further, for air-water system in a GDL the value of the viscosity ratio (M), defined as the ratio of dynamic viscosity of the non-wetting phase (water) to the wetting phase (air), is of the order of 10 1 . Based on the results in References 5,6, the values of Ca and M indicate that the liquid water pore filling in the GDL is capillary driven and leads to a regime called capillary fingering. Experimental investigations of the liquid water transport in GDLs have primarily focused on measuring the saturation as a function of capillary pressure for various commercially available GDLs. [7][8][9] Although useful for characterizing GDLs, these curves do not provide any information about the liquid water distribution. Numerical models extrapolating the average saturation from such curves without accounting for the local water distributions to obtain saturation dependent properties might not necessarily be meaningful as explained in the recent work by Garcia-Salaberri et al. 10 Direct measurement of the effective diffusivity for various commercial GDLs under different compressions and PTFE loadings has also been performed in previous studies. 8,[11][12][13][14][15][16][17] Most of these studies however have focused on measuring the effective diffusion coefficient on dry GDL samples. Correlations between effective GDL diffusivities and local saturations in the GDL were not determined. These correlations are however essential to aid GDL selection for fuel cell operation at high current density under wet conditions. Very few experimental studies [18][19][20] have investigated the effective diffusivity of the GDL as a function of the local saturation. This is likely due to the difficulty associated with controlling the liquid water saturation through the thin GDL for the duration of the experiment.
Microstructure based modeling provides a feasible alternative to study the liquid water and gas transport in the GDL. Advancements in microscopy techniques, such as micro-computed tomography (μ-CT), [21][22][23][24][25][26][27] have enabled the visualization and 3D reconstruction of the GDL morphology. Numerical models can be used to simulate the liquid water distributions on the 3D reconstructions and determine saturation dependent effective transport properties which are more challenging to obtain from experiments. However, experimental data is required to validate the simulated liquid water distributions and effective diffusivities. μ-CT has recently been used to obtain the 3D reconstructions of partially saturated GDLs 10,28-32 thereby enabling the study of the changes in the liquid water distributions with compression, 28 capillary pressure 10,29,30 and temperature gradients. 31,32 These reconstructions provide an ideal way to validate and assess the applicability of the numerical models used to study liquid water intrusion and gas transport in GDLs.
Among these, PNM is the most popular approach to study pore filling and liquid water distribution in the GDLs and to evaluate the corresponding effect on the transport properties. In the PNM approach, the physical domain is idealized as a computational domain comprised of a network of pores and throats with certain radii. The computational domain is generated either from microstructure images 34,44,48,49 or random networks 3,36-39 calibrated against experimental data such as porosity, mercury intrusion porosimetry or saturation-pressure profiles. Due to the low computation cost of PNMs, it has recently been integrated into a full MEA model to predict the local saturation and effective properties of the GDL. 50,51 Despite their usefulness, PNMs do not faithfully represent the pore space, so their predictions of water configurations are only rough estimations. This further reduces the reliability of the gas transport predictions as a function of the local saturation in the GDL.
A more accurate picture of the water distribution can be obtained by simulating water configurations directly on tomography images. The most commonly used method to estimate water configurations is an image analysis technique known as full morphology or morphological image opening. Full morphology (FM) models use the Young-Laplace

F554
Journal of The Electrochemical Society, 165 (7) F553-F563 (2018) equation to describe pore filling in the porous media 33,46,47,52,53 similar to PNMs. However, FM models can be applied to real microstructures of the porous media. FM models use image erosion with spherical balls of a given size followed by image dilation to estimate the liquid water distribution in the microstructure. The radius of the spherical balls is calculated using the Young-Laplace equations based on the capillary pressure.
Image analysis techniques suffer from a few key limitations, largely because a detailed momentum conservation equation is not applied to the gas and liquid domains and a detailed treatment of surface tension forces is not included. The Lattice-Boltzmann method is the most commonly used multi-physics based technique to study water transport in fuel cell porous media. 4,33,[41][42][43][44][45] Eulerian-Lagrangian formulations [54][55][56] and volume of fluid (VOF) methodologies 57-59 could however also be applied. LBM is able to capture the intricate dynamics of liquid water and water droplets with the porous media as opposed to the other methods which assume a quasi-static water front.
Vogel et al. 60 compared PNM, FM and LBM to estimate the capillary pressure-saturation relationship for sintered glass. Their results showed that the PNM approach is sufficiently accurate to predict the capillary pressure-saturation relationship. Further, Vogel et al. 60 also reported that the computational time for PNM, FM and LBM was of the order of 1, 1.5 and 10 4 -10 5 respectively. However, Vogel et al. 60 did not compare the water-air distribution of the three models to imaging data therefore, the accuracy of the water distribution could not be determined from their study. Kim et al. 61 compared LBM and FM simulations for an air-water system in compacted silica sand and compared the predicted wetting and non-wetting phase distributions to X-ray and neutron imaging data. A qualitative comparison of FM and 2D LBM results with X-ray and neutron images showed that the FM model was able to predict similar liquid water distributions in the compacted silica sand as obtained from the X-ray and neutron images while being computationally more efficient than the 2D LBM used in their study. Vogel et al. 60 and Kim et al. 61 however did not study fibrous media, such as gas diffusion electrodes. Agaesse et al. 62 used μ-CT reconstructions of partially saturated GDLs at different capillary pressures to validate numerically simulated liquid water distributions using a pore network model and a full morphology approach. Their work demonstrated the capability of the numerical models to predict capillary pressure-saturation curves and saturation profiles similar to those obtained from the μ-CT reconstructions. However, none of these studies have analyzed the gas transport in the partially saturated media therefore, the effect and importance of the liquid water distribution on predicting gas transport cannot be determined from these studies.
Partially saturated GDL morphologies obtained from PNM, FM and LBM are commonly used in order to estimate effective diffusivity and gas permeability for different GDL microstructures under dry and wet conditions. Therefore, it is also critical to assess the errors on transport property predictions associated with inaccuracies during the water intrusion process. Gas transport in the GDL microstructures has been simulated using PNM, 36,37,63,64 LBM 10,21,30,45,[65][66][67][68][69] and porescale CFD simulations, 23,[25][26][27]70,71 however the obtained results are seldom compared to experimental observation and the sensitivity of transport property predictions to errors in water distribution have not been discussed.
In this article, liquid water intrusion in a dry GDL reconstruction is simulated using FM and compared to experimental μ-CT tomographic reconstructions of the partially saturated GDL sample in order to estimate the validity of this approach. To simulate the liquid water intrusion, a cluster based algorithm based on the full morphology approach is developed. To evaluate the validity of the liquid water distributions obtained from the simulation results, gas transport simulations are performed to estimate the effective diffusivity of dry and partially wetted GDL samples obtained from the full morphology simulations and from tomographic reconstruction using the Fick's second law solver in OpenFCST, 72 an open-source fuel cell simulation package based on the finite element method (FEM). Algorithm and implementation details for the water intrusion and gas transport are discussed in Theory section and the results are discussed in Results and discussion section. This work uses a powerful but simple means of of quantifying the accuracy of the simulation results by performing a voxel-by-voxel comparison of real and simulated water distribution. Effective diffusivity of the dry and partially saturated GDL reconstructions from CFM and tomography are computed and compared to previously reported experimental data from literature. Numerical water intrusion model.-Water intrusion into the dry GDL reconstruction was simulated numerically using a boundary based method. Figure 1a shows a schematic of the boundary based method where water from a circular reservoir region, shown in Figure  1b, representative of the pipe in the experimental setup, is intruded into the sample.
The liquid pressure in the reservoir was increased to intrude the water into the GDL as shown by the arrows in Figure 1a. To simulate the water intrusion and track the progression of the liquid water front into the GDL a quasi-static algorithm based on the full morphology (FM) approach outlined in Reference 53 was developed. The steps of the developed algorithm are described in the paragraphs below and summarized in Algorithm 1. This algorithm uses capillary pressure and the local pore radius of the porous media to predict the movement of the liquid water front.
The inputs to the algorithm are, 1 the GDL microstructure as a 3D array or stack of images ( 0 ), 2 the contact angle (θ), 3 the boundary intrusion points for the liquid water, and 4 the number of discrete steps (n) for incrementing the capillary pressure.
The 3D reconstruction of the dry GDL sample was used as the input to the FM algorithm. The GDL was assumed to be hydrophobic with a contact angle of 110 • (similar to the reported contact angle of PTFE, i.e., 112 • . 36 ) The boundary intrusion points were defined as a circular region with a diameter of 1500 μm at z=0 as shown in Figure  1b (z is the through-plane direction and the direction in which the Algorithm 1. Full morphology algorithm for water intrusion. All capital Greek symbols denote 3D arrays with dimensions of the input microstructure. Read original input microstructure 0 ; Define the intrusion points, contact angle (θ) and number of steps (n); Compute p c for every pore pixel in 0 using Equation 1 Initialize p in as the minimum liquid pressure to start water intrusion Initialize l = intrusion points ( l has a value of 1 at the locations where water is present.) Compute p = max( p c )− p in n while s<1 do procedure IDENTIFY LIQUID WATER CLUSTERS Obtain the set voxels where p c < p in Cluster connected voxels in and store in with each cluster denoted by a unique integer j Identify cluster indices ( j) in connected to l and store in array u end procedure procedure Record l for mesh generation Increment p in = p in + p end procedure end while water was intruded into the sample in the experiments). The number of discrete steps was varied between 8 and 100.
The first step in the algorithm was the computation of pore size distribution (PSD) and capillary pressure ( p c ) required to intrude each pore. The PSD for the GDL microstructure was computed using the sphere fitting algorithm described in Ref. 73. The PSD algorithm first evaluates the distance transform at pore voxels in 3D, i.e., the distance of a pore voxel to the closest solid voxel. Then, the distance transform map is used to find all possible locations where spheres of radius r can be fit (i.e., the distance transform value should be greater than or equal to r ). These locations are treated as sphere centers with the neighboring pore voxels within a distance of r all being assigned a pore radius of r . The process is repeated by incrementing r until the maximum value of distance transform map is reached. Thus, every pore voxel is either the center of a sphere with a radius equal to its distance transform value or part of a larger sphere. This approach is similar to that presented by Agaesse et al. 62 with a difference being that in the present algorithm the radius r is assigned to all pixels within r distance from the center. This eliminates the need to use a morphological dilation 62 to intrude liquid water.
Using the pore radius, the capillary pressure ( p c ) required to intrude each pore was computed using the Young-Laplace equation, where γ is the surface tension of the liquid (water in this case), θ is the contact angle between the solid surface and liquid and r is the pore radius. For implementation purposes, the capillary pressure ( p c ) is defined as, where p l and p g is the liquid and gas pressure respectively. Using this convention for capillary pressure leads to hydrophilic pores having a negative capillary pressure and hydrophobic pores having a positive capillary pressure value. In this work p g is taken as zero.
The propagation and tracking of the liquid water front is the most crucial step in the algorithm. To improve the efficiency of the algorithm a cluster based approach was formulated as described in Algorithm 1. First, all locations with a capillary pressure lower than the specified liquid pressure at a given step were identified and clustered into groups with unique integer values ( j) and stored in a 3D array ( ). This step identified all possible locations where the liquid water could be present at the given liquid pressure p in . The next step was to identify the clusters in which were connected to the existing water filled pores denoted by l . This was performed using a logical AND operation between and l which returned the unique integers (u) for cluster locations where the water could propagate. To intrude liquid water into the identified clusters, the x, y and z co-ordinates of all points in corresponding to the integer values in u were found and the corresponding locations were set to 1 in l .
The saturation (s) at each step is then computed using, where ε( l ) denotes the volume fraction of voxels with liquid water and ε V denotes the porosity of the microstructure. The PSD algorithm combined with the cluster based approach distinguishes it from the FM algorithm previously reported in References 22,46,47,53 which uses a series of image erosion and dilation operations with spherical structuring elements. The water intrusion and PSD algorithms were coded using Python as a part of the open-source package OpenFCST. 72 Gas transport analysis.-Meshing.-Voxel based meshes were generated for the partially saturated GDL microstructures obtained from the microscopy images and numerical water intrusion simulations. The meshes were created by direct conversion of image voxels into mesh cells. Therefore, the mesh cells were uniform and cubic with edge length of 1.33 μm and consisted of three material ids corresponding to solid fibers, void phase and liquid water. For gas transport simulations, a new mesh containing only the cells with void phase material id was used. Due to the large size of the computational domain, i.e., 1.995 × 1.995 × 0.282 mm 3 (nearly 477 million cells at a voxel resolution of 1.33 μm × 1.33 μm × 1.33 μm), performing numerical gas transport simulations on these domains was computationally expensive. To circumvent this problem, the 3D stack of images was scaled using a nearest neighbor interpolation 74 by factors of 2 and 4 to coarsen the voxel resolution to 2.66 μm and 5.32 μm in all directions, respectively. A comparison of the effect of the mesh resolution on the gas transport will be discussed in Dry and wet diffusivity of GDL section.
Gas transport governing equation and solution strategy.-Gas transport in the GDL is assumed to be steady state and analyzed using Fick's second law, where D is the bulk diffusion coefficient of the species i, c tot is the total gas concentration (assumed constant), and x i is the molar fraction of the diffusing species i. For pores filled with liquid water, the diffusion coefficient was assumed to be zero because of the four orders of magnitude difference between the molecular oxygen diffusion coefficient in air and water. 75 The boundary conditions used to simulate the oxygen transport in the GDL microstructure are, and (Dc tot ∇x O 2 ) · n = 0 everywhere else, [5] where x in O 2 and x out O 2 are the oxygen molar fractions at the inlet ( 1 ) and outlet ( 2 ) planes respectively. The gas transport equations were discretized using the finite element method. The resulting system of equations was solved using the conjugate gradient solver in deal.II 76 and integrated in the framework of the open-source package OpenFCST. 72 The effective diffusion coefficient (D eff O 2 ) of the GDL was computed using, where N out O 2 is the total molar flow rate of oxygen at the outlet plane 2 , L is the shortest distance separating inlet and outlet planes, A is the cross-section area and, c in O 2 and c out O 2 are the concentrations of oxygen at the inlet and outlet boundary faces, respectively. The effective oxygen diffusion coefficient is used to compute the formation factor for diffusion using, [7]

Results and Discussion
Water intrusion was simulated on the dry GDL reconstruction obtained from μ-CT images. The reconstructed GDL domain had dimensions of 1995 μm × 1995 μm × 281.96 μm (1500 × 1500 × 212 voxels) with a resolution of 1.33 μm in all directions and a porosity of 65.8%. As shown earlier in Figure 1, water intrusion was simulated from the boundary at z=0 using the cluster based full morphology (CFM) algorithm developed. For this study, the GDL was assumed to be hydrophobic with a contact angle of 110 • . Figure  2 shows the simulated liquid water profiles in the GDL at capillary pressure of 1, 4 and 8 kPa.
Due to the cluster based approach used, water intrusions steps are relatively fast for a GDL domain size with 477 million DOFs. To perform a single water intrusion step at a given capillary pressure the algorithm took an average of 6 minutes on a single core of Intel(R) Xeon(R) CPU E5-2690 v2 with a clock speed of 3.00 GHz. Computation times for water intrusion simulations are seldom reported in literature. Vogel et al. 60 reported that both the computational time and cost (in terms of RAM requirement) for PNM, FM and LBM were of the order of 1, 1.5 and 10 4 -10 5 respectively. Based on the values reported by Vogel et al. 60 the computational time for PNM would be lower than the CFM model used in this study but the computational time for LBM would be on the order of 4-40 weeks on a single core processor. In the section below, three properties are compared between the numerical liquid water intrusion algorithm and μ-CT imaging results, namely: i) capillary pressure-saturation curve, ii) 3D liquid water profiles, iii) gas diffusivity in partially saturated GDLs.
Capillary pressure-saturation curves.-Capillary pressuresaturation curves are commonly used to characterize liquid water transport in porous media. 7,36,46,53,77 Figure 3 shows the saturation as a function of capillary pressure for the simulated water intrusion using the CFM model, the experimental tomographic reconstructions and the ex-situ experimental data in Reference 7 for the same material. For the CFM simulations, a constant contact angle was assigned to all the pore surfaces. Since the local contact angles for the GDL reconstruction were not known a parametric study was performed by varying the contact angle between 95 • -120 • in the Young-Laplace equation. An increase in the contact angle for the simulations resulted in an increase in the capillary pressure required to obtain a given saturation without any change in the shape or slope of the capillary pressure-saturation curve.
Comparison of the simulated capillary pressure-saturation curves obtained using the proposed CFM model to those obtained from μ-CT reconstructions and previously reported experimental data by Gostick et al. 7 for the same material in Figure 3, shows that the CFM model predicted higher saturations for capillary pressures of 4-8 kPa compared to the saturation values from μ-CT reconstructions and previously reported in Reference 7 for all contact angles used. However it must be noted that the sigmoid shape of the simulated capillary pressuresaturation curve resembles the experimental data by Gostick et al. 7 The difference in the saturation values at a given capillary pressure between the simulation and experiments is likely due to the use of a 180 • contact angle for estimating the pore sizes. The different contact angles used in Figure 3 for the CFM simulations are merely a fitting parameter using the Young-Laplace equation. Although the CFM predicted capillary pressure-saturation curve for the sample used in this study is quite different from the μ-CT and literature curves, a FM model has been shown to be able to predict the capillary pressuresaturation curve for SGL 24BA. 62 Therefore, analysis of different samples might be needed to improve the limitations of the CFM model. Schulz et al. 53 proposed an improvement to the FM model by creating fictitious spheres with diameters increased by the cosine of the contact angle to account for the wettability of the walls. Although this is an improvement over the traditionally used FM model, it would be difficult to model the real GDL surfaces where multiple contact angles might exist along the walls of the same pore. The method used in this study to fit the capillary pressure with a contact angle using the Young-Laplace equation is an approximation and a future study would investigate the implementation of a more accurate representation of the local contact angle in the GDL. Figure 4 shows the liquid water distribution at a saturation of 39% at half through-plane depth of the GDL reconstruction using contact angles of 100 • , 110 • and 120 • . It can be seen that for the same saturation (of 39%) the liquid water distribution for the different contact angles is identical in the GDL due to the reasons explained earlier. It must also be noted that for the CFM simulations the effect of the contact angle on the saturation and water distribution can be neglected as long as the wettability (i.e., hydrophobicity of the pores) is maintained and therefore, the contact angle can be used as a fitting parameter. This warrants the comparison of the water distribution between the μ-CT reconstructions and CFM simulations of the partially saturated GDLs at similar saturations rather than capillary pressure which will be discussed in Liquid water distributions section.
Liquid water distributions.-Capillary pressure-saturation curves provide averaged information about the liquid water transport through the whole porous media. However to validate the profile and distribution of liquid water in the porous media, a one-to-one comparison must be made between the numerically simulated and microscopy reconstructed 3D liquid water distributions. To validate the liquid water intrusion algorithm developed in this study, the 3D water distributions from the simulations were compared to the 3D water distributions from μ-CT reconstructions. In order to compare liquid water distributions with the same volume fraction of water, simulations and experiments were compared at similar saturation values. It must be noted that since the simulation uses a discrete pore size distribution of the GDL obtaining identical saturation values as the experiments is difficult. To achieve similar saturation values to the experiments, the number of steps n in Algorithm 1 was increased to 100. Figures 5 and 6 show a comparison of the experimental and simulated liquid water distributions in the x-y plane at 25%, 50% and 75% through-plane (z) depth of the GDL from the inlet at a saturation of nearly 30% and 46% respectively. Comparing the liquid water distribution between the simulations and experiments in Figure 5, it can be seen that the simulated water clusters in Figures 5b, 5d and 5f appear in similar regions as those seen in the μ-CT reconstructions in Figures 5a, 5c and 5e. However the simulated water clusters appear more discrete compared to the liquid water clusters from the μ-CT reconstructions. The discrete nature of the simulated liquid water clusters also exists at the higher saturation of nearly 46% as shown in Figures 6b, 6d and 6f. Additionally, for both the saturations, i.e., 30% and 46%, the simulation predicts excessive liquid water clusters, especially near the outlet, compared to the experimental reconstructions as shown by the comparison of Figures 5e and 5f and Figures 6e  and 6f.
Based on the comparison of the liquid water distributions, two issues exist between the simulated and experimental water distributions. The first one is the presence of discrete capillary fingers of liquid water in the simulations compared to the larger clusters in the μ-CT reconstructions. The discrete liquid water clusters seen in the simulated liquid water distributions in Figures 5 and 6 can be attributed to the algorithm used for estimating the pore sizes for the GDL. Figure  7 shows the pore size distribution and liquid water distribution at s = 32.01% for a GDL slice at 50% through-plane depth. As discussed in Numerical water intrusion model section, the PSD algorithm assigns spherical radii to the pores. Therefore, pore voxels which are closer to the fibers or in between two large pores are assigned a smaller pore radius as shown in Figure 7a. Consequently, the liquid water does not intrude these pores as shown in Figure 7b where the capillary pressure is 3.12 kPa corresponding to all pores with radius bigger than 15.9 μm that are connected to the water network being flooded. Note that this algorithm is similar to that used in previous FM results for fuel cell porous media, e.g., References 22,46,47,62, therefore it is likely that the previous work experienced similar problems. The FM liquid water distributions shown by Agaesse et al. 62 also showed a larger number of clusters compared to the experimental distributions. Since, Agaesse et al. 62 did not show 2D profiles of the liquid water distributions in the GDL slices it was difficult to observe the granularity in the simulated liquid water distribution. Similar discrete capillary fingers can also be seen in the FM simulation results on compacted silica sand by Kim et al. 61 Image processing operations to better discretize the void space into pores using advanced methods analogous to those developed for pore network extraction 49,62 could help improve the granularity of the liquid water distribution. Another improvement to the current model could be the use of Purcell toroid model 78 instead of the Young-Laplace equation to account for the liquid water intrusion. These improvements to the existing CFM model will be investigated in a later study.
The second issue is the prediction of excessive liquid water clusters in the simulations as seen in Figures 5d, 5f, 6d and 6f. This is likely due to the variation of the local contact angle in the GDL pores which has not been accounted for in the simulations.
Although the primary source of error between the simulated and reconstructed liquid water profiles is attributed to the formulation of the numerical model, the uncertainties associated with the experimental setup and μ-CT reconstructions must also be pointed out. It must be noted that the μ-CT reconstructions were obtained using manual segmentation of gray-scale images 30 which might have resulted in an overprediction of the connectivity of the liquid clusters. Details of the experimental methods and image processing procedures can be found in Ref. 30. Inconsistencies in the reconstructions can be seen in Figure 8 where a liquid water cluster (near top edge) is originating from the outlet which is supposed to be in contact with a hydrophobic membrane. This would affect the boundary condition for the CFM simulations where the outlet face in contact with the hydrophobic membrane is assumed to be non-penetrable by the liquid water. Further, the liquid water behavior in the hydrophilic membrane might also influence the boundary condition for the liquid water intrusion at the inlet face. This indicates the need for well controlled experiments together with development of optimal reconstruction algorithms to reduce uncertainties due to image analysis, leaks in the system and evaporation. These would also facilitate the development of more representative numerical models. Figure 9 shows the local average saturation profiles in the throughplane direction for the cases discussed above. Local average saturation profiles are obtained as the ratio of the number of pixels assigned to the water phase to the sum of pixels in the void and water phase in  a 2D cross-section (XY plane) at a given distance from the bottom (inlet plane) of the GDL. The experimental reconstructions exhibit a steeper gradient in the local saturation across 25-75% of throughplane depth. Saturation gradients of 22.17% and 17.97% are observed at s=30.43% and 46.79% in the μ-CT reconstructions, respectively, compared to 11.15% and 6.1% variation in the corresponding simulations. The local saturation profiles from tomography and simulated water distributions exhibit less than 8% discrepancy in saturation between 25%-75% of through-plane depth for both the saturation cases. The tomographic water distributions show a larger variation in the local saturation from inlet to outlet compared to the CFM water distributions. Further, it can be seen that the numerical simulations are affected by the boundaries which result in an extrema near the inlet and outlet. This is likely due to the pore size computation algorithm where the distance transform near the edges is computed by looking In order to quantify the agreement between experimental and numerical results the number of liquid water voxels which are identical in both images as a fraction of the total liquid water voxels is shown in Figure 10. The random distribution line indicates the average agreement that would be obtained if water was assigned randomly to the void space for a given saturation. Since higher saturation means more pores will have water, the probability that a pixel would be assigned to the water phase increases proportionally with saturation. Comparing the agreement of the simulations to the average random line, it can also be seen that the proposed model provides significantly higher accuracy (15-40%) compared to a random water distribution. Figure  10 shows that for saturations between 29-50% the CFM model is able to predict 56-60% of the liquid water voxels observed in the μ-CT  reconstructions. The quantitative one-to-one comparison of the liquid water distributions from the experimental μ-CT reconstructions and numerical simulation using FM model shows that the developed CFM model is able to predict the liquid water intrusion in the GDL for saturations in the range 29-100% with a accuracy of 56% and higher. In the 6-29% saturation region, which is of particular interest for fuel cell simulations, μ-CT data was not available. Therefore, the accuracy of the CFM algorithm cannot be estimated in this region, however it Figure 9. Local saturation profiles in the through-plane direction for the tomographic and simulated liquid water distributions for the cases described in Figures 5 and 6. is likely to be between 40 and 60% based on the results at higher saturations as well as the results from Agaesse et al. 62 where similar agreement is observed using a FM approach.
Although the accuracy of the CFM model for a saturation of 6% is only 22%, this is much higher than a random distribution of liquid water which would have an accuracy of less than 6% (with the value approaching 6% for a large number of realizations). This suggests that the numerical model is able to capture some of the characteristics of the water intrusion even at these low saturations. As discussed earlier, there are several uncertainties associated with the μ-CT reconstructions, primarily the boundary conditions that should be used at the inlet and outlet of the GDL, and simulations parameters such as local wettability and contact angle in the GDL pores. Therefore, the current level of agreement between the simulated and experimental water distributions was deemed acceptable. This figure also provides a quantifiable benchmark to compare future models to study liquid water intrusion.
Gas transport in partially saturated GDLs.-Liquid water distribution predictions were required in order to predict transport in the GDL at varying saturation levels. In order to analyze if the current level of accuracy on water distribution was appropriate to estimate transport properties, the proposed CFM model was used to generate the pore space available to gas transport for partially saturated media and the results of mass transport simulations were compared to the effective diffusivity predictions obtained from simulations on the wet μ-CT reconstructions. This comparison would determine whether the current level of accuracy in the water distribution is appropriate to predict gas transport in the partially saturated GDL.
Dry and wet diffusivity of GDL.-As described in Meshing section, gas transport was simulated on the partially saturated GDL reconstructions by direct conversion of the image voxels to mesh cells. To improve computational efficiency, the cell resolution was coarsened to 2.66 μm in all directions from the original 1.33 μm. Gas transport was simulated using Fick's second law as described in Gas transport governing equation and solution strategy section. The average simulation time for the dry GDL sample with nearly 46 million DOFs (domain size of 1.995 mm × 1.995 mm × 0.282 mm) was 3.6 hours on a single core of Intel(R) Xeon(R) CPU E5-2690 v2 with a clock speed of 3.00GHz. Figure 11 shows the formation factor in the Cartesian directions as a function of the saturation for the partially saturated GDL reconstructions from μ-CT data and numerical simulations using CFM model. Formation factor is defined as the ratio of the effective diffusivity to the bulk diffusivity. At zero saturation, the formation factor in the x, Figure 11. Variation of the formation factor in the x, y and z directions with saturation of the GDL from μ-CT and numerical reconstructions compared to previously reported literature data for partially saturated Toray TGP-H-120 GDLs. 18,20 For the numerical and μ-CT reconstructions, x and y are in-plane directions and z is through-plane direction. Tranter et al. 20 measured the inplane (IP) diffusivity of partially saturated Toray TGP-H-120 GDL with 5% PTFE. Hwang and Weber 18 measured the through-plane (TP) diffusivity of partially saturated Toray TGP-H-120 GDLs with 10% PTFE loadings.
y and z direction is 0.418, 0.404 and 0.224 respectively. The average in-plane (x and y) formation factor is 0.411 for the dry GDL reconstruction which agrees with experimental in-plane formation factor values in the range of 0.31-0.54 17 for the Toray TGP-H-120 sample with 10% PTFE loading and porosity in the range of 0.61-0.73. The through-plane (z) formation factor for the dry GDL reconstruction is 0.224 which is also in agreement with previously reported experimental values in the range of 0.14-0.33 for Toray 120 samples. 11,12,18 The formation factor in all directions decreases with an increase in saturation. At a given saturation, the formation factor for the liquid water distributions from CFM model (solid line) is higher than those obtained from the μ-CT reconstructions (dashed line) in the x and y direction. However the difference in the formation factor values for the x and y direction between the μ-CT and numerical reconstructions of the partially saturated GDLs is less than 20% for saturations below 40%. For higher saturations, the difference in the formation factor values increases to nearly 67% for s=42%. The formation factor in the z direction is nearly identical for both the CFM and μ-CT reconstructions at different saturations.
The higher deviation observed in the in-plane formation factors for saturations above 40% might appear as a surprise because the accuracy of the liquid water distribution for these cases is above 60% whereas for saturations below 40% the accuracy of the liquid water distribution is lower. To analyze the reason for the larger discrepancy in the predicted effective diffusion coefficients at high saturation levels, saturation profiles for GDL slices in the x and y direction are shown in Figure 12 at low and high saturation levels. Saturation profiles for GDL slices in the z direction are given in Figure 9 for the same saturation levels. Figure 12 shows that the μ-CT reconstructions exhibit a much larger variation in the local saturation compared to the CFM simulations at high saturation with local saturation in certain regions as high as 70%. The high saturation local areas lead to a significant drop in the effective porosity in these regions and consequently a sharp decrease in the effective diffusivity. These peaks are not seen in the CFM simulations where the local saturation increases nearly uniformly through the whole domain. This explanation is further strengthened by the loss of a percolating pore path in the x direction in the μ-CT reconstructions at a saturation of 76% which is also not observed in the CFM simulation predictions. Figure 9 shows that large variation in saturation are not observed in the z direction even at high saturation, and in this direction there is a much higher accuracy for the through-plane (z) formation factor. The peaks in μ-CT reconstruction could be due to a higher water affinity toward these regions due to a lower contact angle, or to issues with reconstructed images as shown in Figure 8. Assuming the former, the difference in the growth of the local saturation peaks between the μ-CT reconstruction and numerical simulations is likely due to changes in local contact angles. Experimental knowledge of local contact angles and improvements to the CFM model might help to better capture the gradients in the local saturation profiles. Figure 11 also shows a comparison of the formation factors for partially saturated GDLs from μ-CT and CFM simulations to previously reported literature data for Toray TGP-H-120 GDLs. 18,20 Tranter et al. 20 measured the in-plane diffusivity for partially saturated Toray TGP-H-120 GDL with 5% PTFE. The simulated formation in-plane (x and y) factors for the μ-CT and CFM simulations have a maximum formation factor difference of ± 0.05 compared to the experimental values in Reference 20.
Hwang and Weber 18 measured the through-plane diffusivities for Toray TGP-H-120 GDL with different PTFE loadings. Hwang and Weber 18 data for 10% PTFE had a higher porosity than the μ-CT reconstruction used in this study, 0.73 vs 0.66. The wet TP formation  The present analysis shows that the difference in the IP and TP formation factors between the μ-CT and CFM simulated structures was less than 20% for saturations below 40%. Other algorithms, such as LBM, VOF or Lagrangian-Eulerian formulation, might be able to predict more realistic water distributions but given the good agreement between μ-CT and CFM effective diffusivity predictions this improvement might not justify the increased computational expense.
Effect of voxel size.-The GDL reconstruction used in the current study was extremely large (470 million voxels) and required large computational resources (RAM and CPU time) to perform the numerical simulations. To circumvent this problem coarsening the image resolution to decrease the voxel count (and hence the computational requirement) was explored. Three sections were extracted from the dry GDL reconstruction and coarsened by factors of 2 and 4, corresponding to voxel resolution of 2.66 μm and 5.32 μm in each direction, respectively, to quantify the relative difference in the formation factor for diffusion with respect to the voxel resolution. The three stacks labeled as Stack 1, Stack 2 and Stack 3 had dimensions of 931 μm × 931 μm × 281.96 μm (700 × 700 × 212 voxels), 532 μm × 532 μm × 281.96 μm (400 × 400 × 212 voxels) and 532 μm × 532 μm × 281.96 μm (400 × 400 × 212 voxels), respectively, and were taken from different parts of the dry GDL reconstruction. For the partially saturated conditions, Stacks 1, 2 and 3 were selected from the same location (as the dry GDL reconstruction) from the μ-CT GDL reconstruction at an average saturation of 46.79% and capillary pressure of 6 kPa. Table I shows the formation factors in the in-plane (x and y) and through-plane (z) direction for Stack 1 at different voxel resolutions under dry and partially saturated (s global =46.79%, p c =6 kPa) conditions. It was ensured that coarsening of the voxel resolution did not affect the porosity and percolating pore volume and the change in the porosity, percolating pore volume and saturation was less than 1% of the original value. The error ( ) defined as, The formation factors under dry and partially saturated conditions increase with a decrease in the voxel resolution for all the stacks. Under dry conditions, it was found that applying a coarsening factor of two and four results in an overestimation of the formation factors of at most 4.35% and 22% respectively, for all three cases. Every level of coarsening reduces the computation time by a factor of approximately 10 from the original image computation time of about 25,000 seconds.
For the partially saturated GDL, the local saturation in the three stacks was different from the average saturation of the entire GDL due to the stacks being extracted from different locations of the full GDL microstructure. Stack 1 had the highest local saturation of 62.52% which led to a significant decrease in the formation factor. For the partially saturated Stack 1, the 2X and 4X coarsening led to an average of error of 9% and 33% respectively, in the in-plane formation factor. The through-plane (z) formation factor for Stack 1 was extremely low (order of 10 −4 ) thereby the coarsening led to extremely high errors of 80-360%. For Stack 2 and Stack 3 the local saturations were 17.46% and 44.88% respectively. The coarsening for Stack 2 and Stack 3 showed a similar trend to the dry formation factors. Therefore, the coarsening factor of 2X, used in Dry and wet diffusivity of GDL section, was deemed to be acceptable as it resulted in less than 7% error in the formation factor except for near saturation threshold value. The relatively low error in the effective diffusivity due to coarsening also indicates that the μ-CT reconstructions for the GDL could be imaged with a higher pixel dimension without loss of accuracy for transport studies.

Conclusions
An efficient algorithm for water intrusion based on cluster formation and the full morphology approach was developed and used to simulate water intrusion in a μ-CT reconstruction of a Toray TGP-H-120 GDL with 10% PTFE content. The simulated water distributions were compared to μ-CT reconstructions of the GDL at different saturations using three metrics, namely, capillary pressure-saturation curve, liquid water distributions and effective diffusivity.
Comparison of the capillary pressure-saturation curves between the μ-CT reconstructions and CFM simulations showed that the shape of the capillary pressure-saturation curve was similar to the previously reported literature data for the same sample. A parametric study for the contact angle showed that an increase in the contact angle shifted the curve forward on the capillary pressure axis while the water distributions at the same saturation were unaffected. The low agreement in the capillary pressure-saturation curve was attributed to the limitations of the FM approach which relies on spherical pores. However, a very good agreement between the capillary pressure-saturation curves from FM and μ-CT reconstructions was shown by Agaesse et al. 62 using a similar approach. Therefore, numerical study on multiple datasets is required to identify the reasons for this deviation and improve the CFM model.
A one-to-one comparison of the water distributions at similar saturations between the μ-CT reconstructions and CFM simulations showed good agreement, with more than 56% of the liquid water voxels at identical locations for saturations between 29-100%. Agreement at lower saturations in the GDL would be more beneficial for the prediction of transport parameters during in-situ fuel cell operation 79 because the liquid intrusion would be governed by the underlying catalyst layer or micro-porous layer and proceed until breakthrough of the liquid water into the channels at saturations of 8-25%. [79][80][81] In this study, however, the objective was only to evaluate the accuracy of the CFM model to predict liquid water distributions and transport properties. A major difference between the CFM and μ-CT water distributions, was the appearance of smaller and discrete liquid water clusters compared to the larger water clusters in the μ-CT reconstructions. This was partly due to the sphere fitting algorithm used for the study which fitted smaller spheres near the GDL fibers. Further improvement of the numerically predicted water distributions would require knowledge of the wettability and contact angles for the local pore walls as well as development of better algorithms to characterize the pore sizes. Although this information is difficult to obtain for the commonly used GDL materials, it might be available for radiation grafted GDLs 82 where controlled HI and HO regions can be engineered within GDLs.
Gas transport was studied on the dry and partially saturated GDLs using Fick's second law. The effective diffusivity values computed for the dry GDL sample were found to be in good agreement with previously reported experimental data in literature. For the wet effective diffusivity values, the CFM predicted liquid water distributions were sufficient to predict the in-plane effective diffusivity values with less than 20% error for saturations below 40%. The significant drop in the in-plane effective diffusivity values for the μ-CT reconstructions at saturations greater than 40% was attributed to the high local saturation variations observed in the in-plane direction, which create bottlenecks in the percolating gas network and eventually lead to the loss of a percolating gas network. The reason for high local saturations observed in the μ-CT reconstructions is likely a lower contact angle in some regions which was not accounted for in the CFM simulations where the contact angle was considered uniform. The through-plane effective diffusivities were similar for both the CFM and μ-CT simulations which was likely due to a better agreement in the local saturation profiles in the through-plane direction. The predicted in-plane diffusivities also agreed well with literature data while the predicted through-plane diffusivities were underestimated by a factor of 2 compared to the experimental data. This discrepancy was attributed to the differences in the predicted saturation distributions and material variability which was difficult to analyze in the current study due to only one 3D reconstruction available for the GDL.