Micro-Scale Analysis of Liquid Water Breakthrough inside Gas Diffusion Layer for PEMFC Using X-ray Computed Tomography and Lattice Boltzmann Method

The main objective of this work is to predict the breakthrough pressure of liquid water transport through the gas diffusion layer (GDL) and/or micro porous layer (MPL) used in polymer electrolyte membrane fuel cells. The integration of structural GDL and MPL with Lattice Boltzmann Method is primary focused. The numerical predictions are also compared with experimental data. The interaction between liquid phase and different surface treatments of solid structures controls the evolution of liquid water and the change of capillary pressure. The geometries of GDLs and MPLs were obtained by three dimensional reconstructed micro-structure images from both nanometer and micrometer-scaled high spatial resolution X-ray computed tomography (CT). The predictions of water breakthrough pressure agree with the data observed in the experiment. They also reveal that the breakthrough pressure and liquid water evolution inside the GDL samples are different when the wetting properties of GDL and/or MPL are changed. The detailed microporous property can be obtained using high spatial resolution image from vision

The enhancement of the transport of liquid water and reactant gases in polymer electrolyte membrane fuel cells (PEMFCs) is a critical subject that has been greatly studied due to its critical importance to fuel cell performance at high current densities. [1][2][3] The liquid transport and the concurrent two-phase flow in the gas diffusion layer (GDL) and micro porous layer (MPL) are the most widely studied. In general GDL materials for PEMFCs are made of carbon fiber based cloths and paper. These materials are typically porous to allow for the transport of reactant gases to the catalyst layer as well as transport of condensed water from the catalyst layer. To accelerate the removal of liquid water, GDLs are usually impregnated with non-wetting polymer such as polytetrafluoroethylene (PTFE) to create hydrophobic characteristics. Further, a thin MPL which consists of carbon powder and polymer binders is applied to the GDL side facing the catalyst layer in the membrane electrode assembly (MEA) to further increase cell performance and mechanical stability. The complex structural and chemical heterogeneity of the GDL and MPL make studying the transport of liquid water and obtaining the solution for mass transport losses substantially complicated. 4,5 Many researchers have studied the transport of liquid water through the GDL and MPL in order to develop an understanding of the resistance of reactant gas transport due to water accumulation.  These included the observation of vapor condensation and liquid breakthrough in a GDL using an environmental scanning electron microscope (SEM). 13,14 They proposed a treelike transport mechanism in which micro-droplets condensed from vapor agglomerate to form macro-droplets which eventually flow preferentially toward larger pores and breakthrough. Pasaogullari and Wang 15 also assumed a treelike water transport behavior in GDLs in their two-phase flow model. Bazylak et al. 18 observed the dynamic changes in water breakthrough locations and they described that phenomena using a dynamic and interconnected network of water pathways within the GDL. Manke et al. 20,21 and Hartnig et al. 22 investigated liquid water evolution and transport during fuel cell operation with synchrotron X-ray radiography. They observed an ''eruptive transport'' mechanism in GDL pores near the channels, which they describe as the quick ejection of droplets from the GDL into the gas channels. Both ex-situ and in-situ experiments have demonstrated the existence of low resistance water transport channels within a GDL and that water transport and breakthrough are dynamic processes. However, the morphology of transport channels and the dynamics of water transport in these channels need to be further investigated. Lu et al. 26 reported the liquid water breakthrough across GDLs with and without a MPL. Their ex-situ experimental setup closely simulated a real fuel cell configuration and operating conditions. Their results reveal that recurrent breakthroughs indicate the presence of an intermittent water drainage mechanism in the GDL. They also concluded that the MPL not only limits the number of water entry locations into the GDL but also stabilizes the water paths. Shimpalee et al. 28 investigated several custom GDLs with two layers of MPL under different fuel cell operating conditions. The two micro layers had different properties and particle size. The pore morphology of those samples including pore size distributions and transport properties such as MacMullin numbers were reported and correlated with their performance data.
Recently, high resolution X-ray computed tomography (CT) was used to generate three dimensional (3D) volume rendering of fuel cell materials. The technique has become a powerful tool to describe morphology and study transport properties inside complex structures like GDLs, MPLs, and catalyst layers. [29][30][31][32][33][34][35][36][37] For example, Epting and Litster 30 used microscale X-ray CT to characterize the morphology of the GDL and provide detail analysis of the local oxygen concentration inside the GDL. Hong et al. 31 combined nanoscale X-ray CT with the data from electrochemical diagnostics to enhance the understanding of materials and electrode designs for fuel cell particular under reversal tolerance anode. Zenyuk et al. 32 reported liquid water evaporation rates in different types of GDL. The evaporation rates were measured using in-situ set-up at synchrotron micro X-ray CT. They also used a simplified 3D mathematical model with liquid water evaporation studies from X-ray CT to determine the transport limitation from water evaporation. This team also used the microscale CT to understand the structural properties of GDL under different compression pressure. 33 Holzer et al. 34 applied their quantitative relationship to describe the effective properties including gas diffusivity, permeability, and electrical conductivity for a dry GDL. They used microstructure from 3D X-ray CT to characterize phase volume fractions, geodesic tortuosity, constrictivity and hydraulic radius at different compression levels. The micro-macro-relationships were then established and reported for the in-plane and through-plane directions. Those images also confirms that all pore space inside GDLs is not spherical shape and therefore using only Young-Laplace equation to analyze the water breakthrough pressure and its phenomena will not be sufficient.
Various types of numerical models have been proposed to solve water management problems in the fuel cell. Conventional computational fluid dynamics (CFD) based fuel cell models have been used on the macro-scale but is not suitable for multi-phase water management in the micro-scale gas diffusion media which has complex boundary geometry. A number of works are based on the continuum two-phase flow model,  which describes the flow and mass transport on the basis of Darcy's law. However, the scarcity of experimental data on many of the necessary relationships and parameters, such as the water saturation, dependent relative permeability, effective diffusivity, and air/liquid water capillary pressure, limited the application of these models to GDL materials. A pore-network model that has a long history in the study of porous media like soils and rocks 62 has recently been used in understanding water transport in GDL materials. [63][64][65][66][67][68][69] The pore network model maps a complex pore space continuum onto a regular and irregular network of pore bodies and pore throats. Several works have shown that invasion-percolation process, which is a strongly capillary-driven process at the limiting case of zero fluid velocity, may be an important mechanism for water transport in GDL. However, most of the pore-network models such as in Refs. 63 to 65 focused on the numerical determination of the macroscopic two-phase properties and the relative phase permeability as a function of water saturation, and few works have been done to clarify the mechanism of water transport through a GDL.
In this work, CFD with Lattice Boltzmann Method (LBM) is introduced as an alternative simulation approach to predict the solution of gas/liquid transport in the porous media. LBM is one class of CFD which solves Boltzmann equations to simulate the fluid with collision models. Due to the particulate nature and local dynamics (discrete lattice mesh), LBM is suitable in dealing with complex boundaries, including microscopic multi-phase problems. Garcia-Salaberri et al. 70 also used this simulation methodology from open source code to predict the effective diffusivity of the GDL samples for dry and water-invaded conditions. However, their GDL samples did not include MPL. In this work, the 3D reconstructed images from X-ray CT showing the micro-structure of GDL and MPL are used as model geometry to exercise the computational model. The main objective of this work is to predict the breakthrough pressure of liquid water transport through the GDL and/or MPL and the numerical predictions are compared with experimental data. The effect of wetting parameter also known as contact angle of liquid water on the solid structures is also taken into account. From 3D simulation, the liquid water evolution during water breakthrough process are reported and discussed.

Model Development
Lattice boltzmann method (LBM).-LBM was originally developed as an improved modification of the Lattice Gas Cellular Automata (LGCA) to remove statistical noise and achieve better Galilean invariance for fast flows. 71,72 LBM is one of the most powerful techniques for computational fluid dynamics for a wide variety of complex fluid flow problems including single phase, free surface, and multiphase flow model in complex geometries. The methods concept of streaming and collisions of particles that incorporate the essential physics of microscopic and mesoscopic processes so that the macroscopic averaged properties obey the desired macroscopic equations. 73 Boltzmann's transport equation is defined as follows: where f i is the particle distribution function in direction i, e i is the particle discrete velocity and i is the collision operator. LBM makes use of a statistical distribution function with real variables, preserving by construction the conservation of mass, momentum, and energy. 73 In the most common approach, the collision operator can be approximated by the Bhatnagar-Gross-Krook (BGK) single relaxation time (SRT) from: The Boltzmann's transport equation with a single relaxation time in the Lattice Bhatnagar-Gross-Krook (LBGK) model for collision operator can be written as: is the equilibrium distribution function and τ is the relaxation time which is related to the macroscopic velocity. Usually, the equilibrium distribution function adopts the following expression: where c s is the sound speed (c s = c/ √ 3), ρ is the macroscopic density, u is the macroscopic velocity, δ is the Kronecker delta, α and β subindexes denote the different spatial components of the vectors appearing in the equation, and w i is the weight factor. The multiscale Chapman-Enskog expansion gives us the relation between the macroscopic viscosity (ν) and the relaxation parameter: where x is the lattice space and t is the time step. For the positive viscosity, τ > t 2 is required as a stability condition in addition to the relaxation time around 0.5. In the multiphase fluid flow and multicomponent, red and blue particle distribution function were represent two different fluids (e.g., gas phase and liquid phase, liquid phase and solid phase, and gas phase and solid phase). 73 The total particle distribution function is defined as: where k denotes either the red or blue fluids. The collision operator for Boltzmann's transport equation with a single relaxation time model can be written as: where ( k i ) A is the collision operator similar to the Lattice Bhatnagar-Gross-Krook model in Equation 2, and ( k i ) B is the additional collision operator contributes to the dynamics in the interface an generates a surface tension defined as: [8] where A k is a parameter which determines the surface tension by using the definition of mechanical surface tension (σ), F is the local color gradient, θ i is the angle of lattice direction, and θ f is the angle of local color gradient (e.g., for fluid-solid interface, θ f is contact angle).
The macroscopic fluid pressures are calculated from the equation of stage: LBM schemes are classified as a function of the spatial dimensions m and the number of distribution functions n, resulting in the notation DmQn. The most common schemes in two dimensions are the D2Q7 and D2Q9, while in three dimensions the most used schemes are the D3Q13, D3Q15, D3Q19 and D3Q27. In this work the commercial LBM solver, XFlow 2016, was chosen to perform the calculation. This solver uses the twenty seven velocity D3Q27 lattice model as shown in Figure 2.
The advanced collision operators were approximated by Multiple Relaxation Time (MRT) scheme. In general, collision operators in MRT model schemes are formally defined by: where the collision matrixŜ i j is a b × b diagonal relaxation matrix, m eq i is the equilibrium value of the moment m i and M i j is a b × b matrix which transforms the distribution function to macroscopic moment. 74,75 The collision operator is based on a multiple relaxation time scheme. However, as opposed to standard MRT, the scattering operator is implemented in central moment space. The relaxation process is performed in a moving reference frame by shifting the discrete particle velocities with the local macroscopic velocity, naturally improving the Galilean invariance and the numerical stability for a given velocity set. 76,77 Raw moments can be defined as: [13] and the central moments can be defined as: The computational lattice size was varied from 0.1 μm to 1 μm giving the number of total lattice elements to be 202,312 to 12,234,030 elements, respectively. They were dependent on the size and the resolution of computed geometry. For nanoscale geometry, the domain sizes were in the range of between 1.50 × 10 −5 mm 3 and 1.66 × 10 −4 mm 3 . For microscale geometry, the domain sizes was around 0.42 mm 3 . The liquid water was fed from the top to the bottom of the geometry with the volume flow rate of 43 μl/min for all cases. All sides of the computational geometries were set as wall boundary. According to the nanoscale-to-microscale size of geometry, the time step was also set to be as small as 0.001 to 1.0 microsecond per time step. The maximum computational time was 120 hours on 8 cores in a single node of Intel Xeon 2.8 GHz -256 GB RAM.

Micro X-ray computed tomography (CT).-Collection of micro
X-ray computed tomography data was performed at Beamline 8.3.2 of the Advanced Light Source (ALS) at Lawrence Berkeley National Laboratory (LBNL). 14 keV X-rays were selected using a doublemultilayer monochromator while a 0.5 mm LuAG scintillator, 5× lenses, and a sCMOS PCO. Edge camera provided 1.33 μm cubic voxels and a horizontal field of view (FOV) of 3.3 mm. An exposure time of 300 ms yielded 8000 counts on the 16-bit camera. For each tomographic scan, 1025 projections were collected over a rotation of 180 • . The raw data files were reconstructed into image stacks using Octopus 8.6 and phase retrieval was performed with the Modified Bronnikov Algorithm (MBA). A sample 3.2 mm in diameter was cut from the GDL and mounted between a flat stamp and aluminum stage. All images used in this work were uncompressed. Further details on the sample holder and reconstruction algorithms can be found in our previous publication. 32 For image processing, the FOV was cropped to about 2 × 2 mm. Using Fiji/ImageJ, three phases were identified and separated; gas diffusion layer (GDL) (carbon fibers, binder, and PTFE), micro-porous layer (MPL), and pores. This was done by using Otsu algorithm and by manual determination of the gray-scale threshold values associated with each phase through visual comparison. Once the phases are separated, the STL data was generated for solid and void phases using BoneJ plugin within ImageJ. The information of characterization of samples can also be found in our previous publication. 33 Nano X-Ray computed tomography (CT).-The high resolution 3D images of the entire GDL thickness and the MPL were obtained using a laboratory-scale, nano-scale X-ray computed tomography (nano-CT) instrument (UltraXRM-L200, Xradia, Inc., Pleasanton, CA, USA). The nano-CT utilizes an X-ray source with a rotating copper anode providing an 8 keV X-ray beam. Prior to the imaging process, the samples were reduced in size in radial (in-plane) direction to reduce the beam path length across the sample, thus increasing the signal-to-noise ratio during the imaging. Ideally, all sample are reduced to the same size of that of the field of view. However, for GDL carbon papers, reducing the size impacts the internal structure of these sparse, fibrous materials and internal tomography within so a larger sample is necessary. In this work, a laser micro-scale milling machine (ESI, QuikLaze 50ST2, Portland, OR, USA) was used to mill out a pillar of the GDL sample with a diameter of approximately 200 μm across the entire GDL thickness. Small gold particles were placed on the top of the GDL sample for metrology and sample drift correction during reconstruction. A set of tomographic images were obtained through Zernike phase contrast imaging with large field of view optics with 2 × 2 pixel binning on the X-ray detector, yielding 128 nm pixels and 150 nm optical resolution. The Zernike phase contrast is necessary for imaging of low atomic number materials or low density materials, such as carbon and polymers, which typically present little absorption contrast at 8 keV. The tomographic data was reconstructed with Zeiss's filtered back projection software, resulting in the voxels of approximately 128 nm in a cylindrical volume with a 65 μm diameter and height. The phase contrast imaged was corrected from interface halos characteristic of phase contrast to volumetric images using the algorithm developed in Kumar et al., 78 which corrects phase contrast artifacts in the reconstructed data. Binary thresholding was applied to the 3D image to obtain the segmented structures of the fibers, binders, and PTFE. For the imaging of the MPL volume, the laser micro-mill was used to mill out a pillar that was 20 μm diameter. It is noted that we observed damage on the outer 1 μm of materials, but did not observe any heat damage with in the sample. The geometry and mesh were generated from within the internals of the sample and thus were not influenced by the surface damage. The Zernike phase contrast images were obtained using the high resolution optics with 16 nm pixels and 50 nm optical resolution. The projection images were reconstructed into a cylindrical 3D CT image with a 16 μm diameter and height. By thresholding, the solid structures of the MPL were obtained without phase contrast correction since it is not necessary on fine-scale features where halo artifacts provide sufficient volumetric contrast. 79 Finally, the reconstructed data was visualized in Avizo (FEI Visualization Sciences Group, Hillsboro, OR, USA). Avizo was also used to generate the 3D surface mesh and tetrahedral mesh used in the simulations.
Water breakthrough pressure measurement.-The water breakthrough pressures of the commercially available GDLs were measured in an ex-situ setup, as shown in Fig. 3. In this setup, the GDL was sandwiched between two stainless steel plates. Water was injected into the GDL (or GDL/MPL) through the bottom plate via uniformly drilled holes of 1 mm diameter with 1 mm interspaces in a φ = 25 mm circle. The actual water/GDL contact area was limited to a diameter of 22 mm by using a 125 μm thick adhesive-backed PTFE film (McMaster-Carr), which prevents lateral water leakage from occurring. To ensure the proper GDL compression, a 125 μm thick hydrophilic PVDF membrane (GVWP02500, 0.22 μm, Millipore) of a diameter of 22 mm was placed in the center hole of the PTFE film. This hydrophilic membrane has a negative capillary pressure and is capable of high flow rate of water; therefore, it does not affect the water breakthrough behavior, but acts to distribute water uniformly over the GDL surface. The top stainless steel plate is a mirror image of the bottom plate and was used as the water outlet. GDL samples with φ = 34 mm were die cut and tested in the experiment. The compression of GDL was monitored with a load cell and kept at 200 psi, which is a typical GDL compression in a fuel cell stack.
The water injection rate was controlled by a syringe pump (Harvard Apparatus 11 Elite). After an initial testing of a wide range of injection rates, the water injection rate for the GDL water breakthrough studies was fixed at 43 μL/min, which corresponded to an equivalent water production rate at current density of 2.0 A/cm 2 and to a capillary number on the order of 10 −6 . The liquid pressure, as referenced to the atmospheric pressure, was measured with a differential pressure transducer (Model FDW2AR, Honeywell) and was recorded with a National Instruments DAQ system at 100 Hz. Assuming that the air pressure in the GDL is equilibrated with the ambient pressure, the pressure transducer directly measures the capillary pressure (i.e., the pressure difference between the water and the air across the meniscus within GDL.) Thermal gradients within the GDL were avoided by controlling the temperatures of the top and bottom plates to a common temperature with two electrical heaters. The bottom plate was designed with a water reservoir of about 5 mL. This reservoir volume was much larger than the amount of water injected to GDL (typically less than 0.5 mL); it is therefore reasonable to assume that the water was maintained at a constant temperature during the injection even though the water pipeline itself was not temperature controlled in the experiment.
The GDL materials investigated in the work and the measured water breakthrough pressures are listed in Table I.

Results and Discussion
Water breakthrough in gas diffusion layer without micro porous layer.- Figure 4 shows a single piece from a GDL sample (SGL 24BA) reconstructed by Nano X-ray CT. Figure 4a show raw data of selected unprocessed virtual slices from reconstructed scans. Figure 4b presents an image of 3D rendering with the size of 65 μm diameter. This image also provides the separation of 3 carbon fibers and PTFE/binder structure. It is noted that the PTFE/binder structure shown in this image is noisy and some binders could not be clearly detected. These undetected binders found in this work is around 10% of the total volume of the sample. Figure 4c is the computational geometry converted from 3D rendered image and it is also split into 2 solid parts, fibers and PTFE/binder. Therefore the properties of both solid parts can be individually input (e.g., contact angle of liquid  water to the solid surface). In this section of study, the contact angle of PTFE/binder structure was varied at 40 • , 90 • , and 140 • and the contact angle of carbon fibers were fixed at 90 • . It is noted the change in wettability of binder/PTFE can control the different surface treatment during the manufacturing process. The liquid water was fed from the top to the bottom of the geometry with the same condition as given in the experimental setup (43 μl/min). Figure 5 shows the prediction of liquid water evolution from LBM with different wettability of PTFE/binder structure. The maximum time for these calculations was 23 hrs. The shape of liquid water and its growth are different with the change of contact angle of PTFE/binder structure. The high hydrophobicity case (140 • contact angle) presents significant different appearance compared to the others (90 • and 40 • contact angles) as it needs more energy to push the water through the GDL. The evolution of liquid water is high sphericity for higher hydrophobicity (non-wetting) of PTFE/binder as liquid water has high contact angle on the structure. The pressures when the liquid water breaks through the GDL from all three predictions are given in Fig. 6. This graph shows pressure versus dimensionless time (tD) for three different wettability of PTFE structure. This pressure was numerically measured by the difference between the pressure of liquid water at the location where the liquid water is fed and pressure at the outlet. So when the liquid water is moving through the GDL and the exit location is occupied by air, this pressure is also known as capillary pressure. Once the liquid water passes through the GDL and continuously flows inside the pore-network of GDL toward the exit, the  pressure difference is stagnation pressure (sum of dynamic pressure and static pressure). Also in this work, dimensionless time is defined as the ratio of local time to the time approaching steady state. All plots show similar profiles as the capillary pressure rises to the maximum when the liquid water passes through the GDL then it decreases to the steady state value as capillary pressure no longer holds and it changes to stagnation pressure as liquid water flows through the pore structure of GDL. The capillary breakthrough pressure was chosen at the maximum pressure for each case. The case of high contact angle or non-wetting of PTFE structure gives the highest breakthrough pressure of 4000 Pa. The breakthrough pressure decreases as the PTFE structure changes to the more wetting material. However, the case of 90 • contact angle shows slightly higher water breakthrough pressure than 40 • contact angle.
The next simulation was computed using the same GDL sample but the image was taken and rendered from the Micro X-ray CT. Figure 7a the raw data of cross-section tomograph for grey-scale and thresholded data of SGL 24BA. 33 Mostly, Otsu algorithm was used to segment the data. Although the binder shows some microporosity, its nano-porosity is not resolved with micro-CT. 33 Figure 7b shows the detailed structure of GDL and it has a much larger size than the image taken by Nano X-ray CT but this image has a lower resolution (i.e, 2.58 × 1.20 × 0.24 mm 3 ). As shown in the figure, both fibers, PTFE structure, and binders are combined to one solid piece. In order to compare the prediction of breakthrough pressure with the previous results that were calculated using the smaller geometry from Nano X-ray CT, three random locations of the GDL in Fig. 7b were chosen and they all have the same volume as the volume of geometry imaged from Nano X-ray CT. Figure 7c shows the three GDL structures that were randomly picked from different positions shown Fig. 7b (Positions 1, 2, and 3). These GDLs have different structure and porosity. Positions 1, 2, and 3 have porosity of 0.7, 0.5, and 0.6, respectively. As mentioned above that the image taken by Micro X-Ray CT, the solid parts of GDL cannot be separated due to its limited resolution. Therefore only the wetting property of the single solid part can be input into the model. Figure 8 shows the prediction of pressure profile against dimensionless time for three different positions as given in Fig. 7. Once again, the maximum pressure of each profile represents the water breakthrough pressure. This figure also provides the snapshots of water evolution inside the GDLs that have the same contact angle of 90 • . Note that the longest time of calculation in these studies was 12 hrs. The overall predictions give the same summary as given in the simulation using higher resolution image from Nano X-ray CT. The maximum breakthrough pressure is shown at the high contact angle (i.e., 140 • ) and both contact angle of 90 • and 40 • show similar pressure profile. The GDL at Position# 1 where the solid structure is loose and has the highest porosity gives the maximum breakthrough pressure of 1750 Pa with stagnation pressure of 1450 Pa for contact angle of 140 • and the minimum of 1100 Pa with stagnation pressure of around 1020 Pa for the contact angle of 40 • . The water evolution also confirms that the liquid water simply move from the top to the bottom of GDL. At Position # 2 where the GDL has the densest of solid structure and the lowest porosity, the highest breakthrough pressure rises to 3650 Pa with stagnation pressure of 1750 Pa and the lowest breakthrough pressure increases to around 1600 Pa with stagnation pressure of around 1370 Pa. With the lower wettability of solid structure, the breakthrough pressure seems to be more sensitive to the porosity of the GDL. The snapshots of the water movement reveal that it requires more pressure than other positions to push liquid water through the GDL. The liquid water also needs to fill up the top before moving down to the bottom due to the high density of solid material located at the top of the GDL. The GDL at Position# 3 has the porosity value of 0.6 which it is not too coarse or too dense. It has the dense structure at the top and bottom but has high porosity in the middle. The maximum pressure at contact angle of 140 • is 2100 Pa with stagnation pressure of approximately 1200 Pa and the minimum pressure at contact angle of 40 • is 1200 Pa with stagnation pressure of around 820 Pa. The water transport behavior is similar to the case of GDL at Position# 2. But once it passes through the dense area of solid structure at the top, the liquid water easily move to the bottom and then additional pressure is created to pass through it. The average maximum breakthrough pressure from those three GDL at the contact angle of 140 • is around 2500 Pa and the average minimum breakthrough pressure at the contact angle of 40 • is around 1320 Pa. Moreover the average maximum stagnation pressure from those three GDL at the contact angle of 140 • is around 1467 Pa and the average minimum breakthrough pressure at the contact angle of 40 • is around 1070 Pa. From the experiment described above (i.e., Water Breakthrough Measurement), the breakthrough pressure of SGL 24BA GDL is 1770 Pa and the stagnation pressure (steady state pressure) is around 1120 Pa which are in the range of model predictions.
The model prediction of water breakthrough pressure was extended to the larger size of SGL 24BA GDL taken from Micro X-ray CT as shown in Figure 9. This GDL is around two orders of magnitude larger than the size taken from Nano X-ray CT. The pressure profiles of the three cases at different contact angles vary from the smaller size in Fig. 8  shows liquid water spreading layer by layer of the GDL from top to bottom.
Breakthrough pressure of gas diffusion layer with micro porous layer.-In order to predict water breakthrough pressure of GDL with MPL, the structure of MPL needs to be included into the model. Figure 10 shows the 3D reconstructed GDL, MPL, and PTEF/Binders used in this simulation. Figure 10a show the structure of MPL of SGL 25BC from Nano X-ray CT with the dimension only 3.375 cubic micron. But the fiber and PTFE shown in Fig. 10b has the dimension of 15,625 cubic and it has the area of 277 times larger than the geometry of MPL shown in Fig. 10a. Moreover, there is an evident crack on MPL as provided in Fig. 10c. Therefore, before creating the complete geometry of GDL and MPL, the surface area of MPL from Fig. 10a must be increased by manually performing the periodic repetition of the Nano X-ray CT data. So that the MPL has the same surface area as GDL as shown in Fig. 10b. Then the space was also manually made to represent the crack shown in Fig. 10c. Figure 10d presents the complete geometry taken by Nano X-ray CT that was used for water breakthrough simulation using LBM. Due to the complexity of MPL micro-structure, the maximum computational time was increased to 100 hrs. Figure 11 presents the snapshots of water evolution inside the GDL and MPL. As the MPL is located on only one side of the GDL (i.e., between GDL and catalyst layer), the water breakthrough pressure prediction was performed at both sides where the liquid water was pushed from MPL to GDL (Figure 11a) and from GDL to MPL (Figure 11b). With the 3D structure taken from Nano X-ray CT, the fiber, PTFE/binder, and MPL structure can be individually separated.
In this study, only wettability of PTFE/binder structure was varied and the rests were kept constant at 90 • contact angle. Both Figures 11a and  11b show different progression of liquid water during breakthrough process. When the liquid water was pushed from MPL to GDL, liquid water spread all over the MPL before it moves toward the GDL. This is because the pressure required to push the liquid water thru the MPL structure or the crack in MPL is higher than the pressure used for liquid water spreads on the top of MPL. But when the liquid water moved from GDL to MPL, liquid water easily transported and grew through the fiber and PTEF/binder structure inside the GDL then liquid water accumulated over the MPL layer and moved past MPL once it reached the breakthrough pressure. Figure 11c shows the close-up views of liquid water passes through MPL and crack from both in Figure 11a and Figure 11b. Figure 12 presents the pressure profiles during the movement of liquid water when it was pushed from MPL to GDL (Figure 12a) and from GDL to MPL (Figure 12b). The contact angle of solid fiber of GDL was fixed at 90 whereas PTFE/binder of GDL and solid structure of MPL were varied at 40 • , 90 • , and 140 • . The overall profiles are similar to the predictions of GDL without MPL from Nano X-ray CT as the pressure increases until it reaches the breakthrough pressure then it decreases to the steady state value. The contact angle of 40 • gives the lowest pressure profile following by the pressure profiles at 90 • and 140 • contact angle. When the liquid water was pushed from the MPL to the GDL as shown in Fig. 12a, the increase in pressure with dimensionless time for 140 • and 90 • contact angle values are similar as liquid water was spread all over the top of MPL before it moves to the GDL. Because the required pressure to push the liquid water completely through the MPL and the crack is higher than the pressure needed to spread liquid water on top of the MPL. After the liquid is fully spread out these two cases give different pressure profiles. For the case at a 40 • contact angle, the pressure profile is different compared to the other two cases. This is due to the hydrophilic nature of the solid surface which allows liquid water to flow through the MPL and crack with less pressure than non-wetting surface. It gives the maximum pressure of 12.5 kPa before it drops to around 10.0 kPa at stable pressure. The case at a 90 • contact angle shows the maximum pressure of 21.0 kPa and then the pressure decreases to a steady state value of 19.5 kPa. For the case of a 140 • contact angle, the pressure increases to reach the maximum at 26.5 kPa then it slightly reduces to reach the steady value at 24.0 kPa. Figure 12b shows the pressure profiles with different contact angles when the liquid water was pushed from GDL. The profiles of all three cases are similar but with different values. This is because the liquid water transports from the coarse structure of GDL through the fine/dense structure of MPL and therefore these pressure profiles are proportional to the wettability of porous material. The breakthrough capillary pressure is 26.0 kPa, 20.0 kPa, and 17.0 kPa for the contact angle of 140 • , 90 • , and 40 • , respectively. It is noted that the overall values of pressure profiles and breakthrough pressure are much higher than the experimental data (i.e., For SGL 25BC, the breakthrough pressure is 7600 Pa and the stagnation pressure is 4750 Pa). This could be due to 1) the sample size used in experiment is much larger than it is used in the LBM model and 2) the MPL used in experiment has multiple cracks with different sizes whereas there is only one artificial crack in the MPL geometry and it much smaller than what was observed due to the geometry size used in the model. It is noted that Shimpalee et al. 80 demonstrated the water evolution and its breakthrough pressure inside the GDL with no-crack MPL. Their predictions give very high breakthrough pressure as around 10 times of experimental value. The reconstructed micro-structure of SGL 25BC GDL from Micro X-ray CT geometry was also used to exercise the model. Figure  13a shows cross-section tomograph, thresholded data with blue represents fiber and red denotes MPL, and 3D rendered image of SGL 25BC (GDL and MPL). 33 As mention earlier that the detailed structure of GDL/MPL from Micro X-ray CT gives only one single solid domain. Therefore, the MPL from this geometry was converted to solid layer with several cracks also shown in Fig. 13a. In order to predict the liquid water transports through both MPL and cracks, the cracks were preserved while the solid part of MPL was replaced by porous model and its properties were taken from the measurement in 3D reconstructed micro-structure of MPL of SGL 25BC from Nano X-ray CT as shown in Fig. 10a. In the LBM used in this work, the porous model applies the Darcy-Ergun formulation to calculate pressure drop ( p) across the media: 73 where k is the permeability, v is the velocity normal to the porous surface, c E is the Ergun coefficient of 0.2, and th is the thickness. The porosity of the MPL is 52% and the permeability is around 1.0 × 10 −15 m 2 for all directions. Figures 13b and 13c present the snapshots of water transport during the water breakthrough progression. Figure  13b shows the snapshots when the liquid water was pushed from MPL to GDL and Figure 13c shows when the liquid water was pushed from GDL to MPL. The contact angle of solid structure was set at 90 • for both figures. The behavior of liquid water is similar to what it was shown in Figs. 11a and 11b using the geometry from Nano X-ray CT. Figure 13b shows the liquid water spreads over the top of MPL before it breaks through the MPL and then GDL. In Figure 13c, the liquid water grows throughout the GDL structure and stays on top of MPL before it breaks through the cracks and MPL. With this approach, the maximum calculation time was reduced to 36 hrs. Figure 14 presents the pressure profiles during water evolution starting from MPL to GDL (Figure 14a) and from GDL to MPL (Figure 14b). The contact angle of solid structure was varied over the contact angle of 40 • , 90 • , and 140 • . The overall profiles are similar to the predictions of GDL without MPL (i.e., SGL 24BA) as the pressure increase to reach the breakthrough pressure then it decreases to the stable value. The wettability of 140 • gives the highest pressure profile following by the pressure profiles of 90 • and 40 • contact angle conditions. When the water was pushed from the top of MPL to the GDL as shown in Fig. 14a, the increasing in pressure with dimensionless time for all three contact angle values are similar as liquid water was spread all over the top of MPL. Then all three cases give the pressure profiles differently. The case of 40 • contact angle gives the maximum pressure of 4500 Pa before it drops to 4000 Pa at the stagnation pressure (stable pressure). The case of 90 • contact angle shows the maximum pressure of 6000 Pa and then the pressure decreases to 4500 Pa. For the case of 140 • contact angle, the pressure increases to reach the maximum at 7800 Pa then it slightly reduces to reach the steady value at 6900 Pa. As already mention, the maximum pressure is the breakthrough pressure or capillary breakthrough pressure. Once liquid water pass through and occupy all pore networks in the GDL/MPL, the capillary pressure will no longer hold thus, reducing to static liquid water pressure. Figure 14b shows the capillary pressure profiles with different contact angles when the liquid water was push from GDL to MPL. Once again, the pressure profiles are similar to the predictions from Nano X-ray CT structure. The breakthrough capillary pressure is 8000 Pa, 5700 Pa, and 5000 Pa and the stagnation pressure is 6000 Pa, 4600 Pa, and 3630 Pa for the contact angle of 140 • , 90 • , and 40 • , respectively.
The values of breakthrough pressure from both orientations shown in Figure 14 are close to the value observed in the experiment especially at the contact angle of 140 • . With the range of breakthrough pressure with contact angle, the contact angle of 135 • from Fig. 14a and 131 • from Fig. 14b represent the angle that are needed to represent the breakthrough pressure measurement. The predictions of breakthrough pressure are also much lower than the predictions from geometry taken by Nano X-ray CT. This is because the size of SGL 25BC used in the model is compatible to the one used in the experiment even though the detailed PTFE/binder has been ignored. Further, the number of cracks and their size are same as the real sample as the resolution and the scale from Micro X-ray CT can capture those details.

Conclusions
In this work, a three dimensional numerical simulation using the Lattice Boltzmann Method has been introduced to study the water transport inside the porous material used in polymer electrolyte membrane fuel cells. The liquid water breakthrough dynamics across the GDL with and without a MPL were predicted and compared with an ex-situ experimental data. The geometries of 3D reconstructed microstructure GDL and MPL from both microscale and nanoscale X-ray CT were used to exercise the model. From the images taken by both resolutions, the overall predictions are in the range of data observed in experiment. The motion of liquid water in the pore structure MPL can only be predicted using the geometry taken by Nano X-ray CT. However, the size of 3D reconstructed MPL requires tremendously scale up in order to create a complete GDL with MPL model geometry. The simulation technique by applying the transport properties of MPL from Nano X-ray CT to porous layer of MPL from Micro X-ray CT appears to predict well and it will require less computational resources. The breakthrough capillary pressure and the evolution of liquid water inside the GDL samples are different with the change of wetting properties of GDL and/or MPL. Due to the irregular and inconsistent shape of pore space in porous material used in PEMFC, it is impossible to only solve Young-Laplace equation and provide the qualitative conclusion. The LBM technique can successfully simulate to two-phase transport inside the porous media using the geometries from both Nano and Micro X-ray CT depending on the purpose of the research. The properties of carbon fiber and other components such as binders and PTFE can be separated using high resolution image from Nano X-ray CT but the sample will be very small whereas, Micro X-ray CT give much larger sample and wide-ranging vision of porous material but the detailed structure inside cannot be observed and separated.