Investigating Phase‐Change‐Induced Flow in Gas Diffusion Layers in Fuel Cells with X‐ray Computed Tomography

The performance of polymer-electrolyte fuel cells is heavily dependent on proper management of liquid water. One particular reason is that liquid water can collect in the gas diffusion layers (GDLs) blocking the reactant flow to the catalyst layer. This results in increased mass-transport losses. At higher temperatures, evaporation of water becomes a dominant water-removal mechanism and specifically phase-change-induced (PCI) flow is present due to thermal gradients. This study used synchrotron based micro X-ray computed tomography (CT) to visualize and quantify the water distribution within gas diffusion layers subject to a thermal gradient. Plotting saturation as a function of through-plane distance quantitatively shows water redistribution, where water evaporates at hotter locations and condenses in colder locations. The morphology of the GDLs on the micro-scale, as well as evaporating water clusters, are resolved, indicating that the GDL voids are slightly prolate, whereas water clusters are oblate. From the mean radii of water distributions and visual inspection, it is observed that larger water clusters evaporate faster than smaller ones.


Introduction
The performance of polymer-electrolyte fuel cells (PEFCs) and other multiphase flow technologies is significantly dependent on liquid-water management [1][2][3][4]. This is particularly true for PEFCs at low operating temperatures and during startup operations due to hindered reactant delivery by water in cathode [5][6][7][8][9]. Because of the exothermic oxygen reduction reaction (ORR) at the cathode, a thermal gradient develops during operation in the through-plane direction, with the hottest location in the catalyst layer (CL). At higher temperatures (~80 °C), this thermal gradient, in combination with the dependence of vapor pressure on temperature, promotes removal of water in a vapor form [1,4,5,[10][11][12][13]. Water vapor within the CL travels through the gas diffusion layer (GDL) to the gas channels (GCs) (see Figure 1), where it condenses due to the decrease in temperature. This type of flow, which is due to the evaporation and condensation of water, is known as phase-change-induced (PCI) flow [4,5,14].
Wang and Wang [4] were the first study in the PEFC field to investigate this phenomenon using a multidimensional model that incorporated electrochemical heat generation, phase-change, two-phase flow, non-constant gas phase pressure, and the entire PEFC domains. By doing so, the authors in [4] also clearly identified the importance of thermal gradients as the driving force for PCI flow within the porous media of the fuel cell. Although water is removed in the vapor phase, depending on a PEFC's operating temperature, a fraction of the total water has to still be removed in the liquid phase. Thus PEFCs experience two-phase water flow and, consequently, substantially coupled heat and mass transport. As such, effective water management requires an understanding of the interaction between pressure-driven, capillary-driven, and PCI water transport [4,5]. Phase change is not a drive potential or force like pressure and capillary forces, however, the term "PCI" has become the common name in literature for heat-driven mass transport of water by evaporation and condensation in a temperature gradient.

Figure 1 Goes Here
The GDL is a porous fibrous component of PEFCs responsible for the transport of electrons, water byproduct, gaseous reactants, and heat [15]. It is made from carbon fibers which are assembled to form either nonwoven paper, woven cloth, or felt. With pores on the order of 10 μm, these materials have porosities typically ranging from 65 % to 90 % and thickness around 200 -400 μm [16,17]. Generally, cell compression influences the GDL's structure [15] and performance during operation. Because carbon fibers are naturally hydrophilic, GDLs are typically treated with 5 -20 % of polytetrafluoroethylene (PTFE). Due to non-uniformities in the coating, there is a mix of hydrophilic and hydrophobic pores, which causes the overall structure to possess mixed wettability [18,19]. As with most porous media, heat and mass transport properties depend on local morphology in addition to bulk material properties.
Most scientific work concerning transport in porous media has been conducted in the fields of civil and petroleum engineering [20][21][22][23][24]. Although this provides a starting point, there are a number of notable differences between the systems studied in those fields and thin materials such as GDLs and CLs.
It is necessary to re-examine each of the various transport mechanisms as they pertain to engineered systems [18]. To this end, much has already been accomplished for transport mechanisms guided by capillary, convection, and gravitational forces. Non-isothermal phenomena, on the other hand, remain an area that is not well researched [1,5,18]. Amongst existing non-isothermal studies, most do not address multiphase flow; let alone phase change [1]. Furthermore, those that do address multiphase flow are typically simulation-based [18] due to difficulties with an experimental approach [1].
Previous studies have shown PCI flow to be a significant contributor to overall water transport within PEFCs. For instance, Weber and Newman [1], through use of one-dimensional simulations, showed that non-isothermal effects are significant when feed gas flows are or become saturated. According to their results, net evaporation/condensation accounts for only 2.6 % of overall heat generation within a fuel cell. However, the heat generated/consumed by each individually is approximately 100 times that of the net contribution. Additionally, their work shows that a thermal gradient of only a few degrees is required across the GDL to completely remove product water from the CL, with larger thermal gradients needed at lower temperatures. As noted by Kumbur and Mench [18], the GDL provides one of the largest thermal resistances in a PEFC and therefore may experience a temperature gradient in excess of 5 °C.
Kim and Mench [5] conducted an experimental study of PCI flow in which they tested various membrane-GDL combinations. It was found that PCI flow does dominate net water flux at high temperatures (80 °C). Furthermore, it was shown that incomplete saturation of the porous media is key to determining whether or not PCI flow will occur.
Micro-CT with a resolution of 1.3 μm is well fit to non-destructively visualize three-dimensional GDL structures and water filling of GDL pores [15,27,33]. Recent studies indicate that, during PEFC operation, liquid water occupies less than 50 % [34] of the GDL pore volume because of the GDL's hydrophobic treatments, and, in the absence of temperature gradients, capillary fingering is the predominant liquid-water-transport mechanism [32,35]. Previously, X-ray CT was used to study the evaporation of water within GDLs under constant temperature. It was found that evaporation rates at water saturations higher than 10 % scale with the surface area of water and are diffusion limited [15].
In this study, a novel X-ray CT technique to explore PCI flow within a PEFC is presented. Coupled measurements of temperature, thermal gradients, and thermal conductivity are combined with visualizations of GDL morphology and water distribution. The overall results of this study contribute to the general understanding of evaporation phenomena in porous media pertaining to PEFCs.

Sample Apparatus
A custom apparatus (Figure 2), was designed and fabricated to conduct the experiment at the synchrotron X-ray CT beamlines. The design aimed to control sample compression, temperature, temperature gradient, and water capillary pressure. The apparatus adheres to size restrictions for the two different X-ray CT beamlines where data was collected and achieves a balance between the needs for structural stability and an un-obstructed view of the sample. The sample sits inside a polyetheretherkeytone (PEEK) ring. Figure 2 shows a schematic of the experimental setup, its three-dimensional rendering, and photographs of the apparatus at the two synchrotron beamlines.
The upper portion of the apparatus consists of a PEEK support structure, nitrile-PVC tube insulation, stainless steel compression cap, stainless steel or aluminum piston, and copper water injection tube. The piston serves as a thermal conductor and a water pathway to simultaneously heat and fill the sample with water. The compression cap allows the piston to be pressed against the sample so as to ensure proper contact. The lower portion of the apparatus consists of a PEEK support structure, nitrile-PVC tube insulation, stainless steel or aluminum piston, and copper water cooling coil. The lower piston serves as a thermal conductor to remove heat from the sample. With regards to the upper and lower pistons, stainless steel was used for the initial data set. However, due to the low thermal conductivity of stainless steel, aluminum was used for further experiments. See Table 1     IL, USA. The data presented in this manuscript is a subset of that from ALS. The other data from ALS will be presented with the data from APS in a future manuscript.

ALS
Image acquisition was conducted using a 500 μm LuAG scintillator, 5x lenses, and a sCMOS PCO.
Dimax camera. This resulted in 2.2 μm cubic voxels and a horizontal field of view (FOV) of 4.4 mm. A double-multilayer monochromator was used to select a beam energy of 22 keV. Each scan was performed over a rotation range of 180° with 1025 projections and an exposure time of 40 ms.

APS
Image acquisition was conducted using a 20 μm LuAG scintillator, 5x lenses, and a sCMOS PCO. Edge camera. This resulted in 1.33 μm cubic voxels and a horizontal FOV of 3.3 mm. A double-multilayer monochromator was used to select a beam energy of 25 keV. Each scan was performed over a rotation range of 180° with 1500 projections, an exposure time of 50 ms, and a total scan time of 3 minutes.

Materials and Setup
SGL10BA (SGL CARBON GmbH -Fuel Cell Components, Meitingen, Germany) was used as the GDL sample in this study. This sample was chosen because it was previously well studied with X-ray CT [15,32] and is easy to handle. In order to obtain a larger thermal gradient, the sample consisted of a stack of two GDLs. Water capillary pressure was controlled by attaching a flexible tube to the water injection tube and then adjusting the water column height to 2 -3 cm. Water vapor was exhausted from the apparatus through the unsealed interfaces between the sample ring and the PEEK supports. After water injection and before each tomography scan, temperatures and heat flux were recorded. The scans were timed to have these measurements every 13 minutes. For each heat flux, 4 -5 measurements and scans were done; after which the heat flux was increased by stepping the voltage to the heater by 1 V on the power supply. Data collection for a given case was stopped when water completely evaporated within the GDL (as observed with X-ray CT scans). Each case lasted 50 -200 minutes. Cooling was accomplished by using flexible tubing to connect the cooling coil to a refrigerated water bath. Temperature data was obtained through use of a thermocouple data acquisition board (product number IPDAS TC from CyberResearch Inc., Branford, Connecticut, USA). Temperature drops across the apparatus ranging from 15 °C to 28 °C were generated by adjusting the voltage applied to the resistive heater. The setup was insulated resulting in negligible heat escape through the top and sides of the PEEK, ensuring one-dimensional heat transport through the sample.

Data Reconstruction
For data from ALS, preprocessing was conducting using Fiji/ImageJ [36]. This was followed by the Modified Bronnikov Algorithm (MBA) to retrieve phases and reconstruction using Octopus 8.6 [37]. For APS data, all steps of reconstruction were conducted using the TomoPy package (version 0.1.15) for Anaconda/Python. First the sinograms were normalized to the white and dark field projections. Then they were normalized to the background intensity using a scaling factor based on 10 pixels from the right boundary and another 10 pixels from the left boundary. Horizontal stripes were removed using the sym16 wavelet filter with 10 discrete wavelet transform levels and a Fourier space damping parameter of 1. For the actual reconstruction, the Gridrec algorithm was applied [38].

Segmentation and Results Collection
Post reconstruction processing was conducted using Fiji/ImageJ [36]. All images were cropped to include the region of interest but exclude the significant reconstruction artifacts along the edges. Separation of three phases (voids, fibers, and water) was obtained through use of a dry reference image for each sample. For these reference images, voids and fibers were separated using the Otsu algorithm. For images of saturated samples, manual threshold determination was 9 used to separate voids from fibers and water. Manual thresholding was conducted by setting various thresholds and the visually judging the balance between retaining water and removing fibers. The corresponding reference image was then subtracted from the sample image to isolate water. Next, ImageJ's "Open" operation (erosion followed by dilation) was used to remove noise from the images of the isolated phases. The operation was repeated 8 times for ALS data and 6 times for APS data. These numbers of iterations were determined by testing multiple numbers of iterations and then selecting those that visually provided the best balance between removing fibers and retaining water. In both cases, the minimum neighbor count was 4. See "Supplemental Information" section 4 for a sample image at various stages in this process.
Through-plane porosity was determined by counting the number of background pixels in the reference image. This was done separately for each slice; thus corresponding to depth into the GDL.
Through-plane water volume fraction was determined in a similar manner but using the isolated water images. The through-plane saturation was then calculated from the porosity and water content data.
Both pore size distributions and water cluster size distributions were determined using the Local Thickness plug-in [39]. Local Thickness utilizes sphere fitting and assigns a pixel value equal to the radius of the largest sphere whose domain contains the pixel. Classification of pore/water cluster shape was conducted using the ellipsoid factor (EF) method [40] from the BoneJ plug-in [41]. The EF method utilizes ellipsoid fitting to calculate an index (see "Supplemental Information" section 3.3 for details; Figure S5 provides example ellipsoids) used to classify pores/clusters according to shape.

Thermal Conductivity
Assuming one-dimensional heat flow though the setup, heat flux and temperature are coupled by the following formulation of Fourier's Law.
where q is heat flux per unit area, k is thermal conductivity, ΔT is temperature difference, and Δx is position difference. Thermal resistance, R [K W -1 ], can be used to rearrange Equation (1) as follows were Q is heat flux and A is cross-sectional area: Combining Equation (3) and Equation (5, Equation (6 can be derived and solved for the GDL thermal conductivity: where Δx samp is sample thickness, k samp is sample thermal conductivity, A samp is GDL cross-sectional area, Q tot is total heat flux through the apparatus, Δx PEEK is the PEEK ring thickness (approximately same as GDL thickness), k PEEK is the thermal conductivity of PEEK, and A PEEK is the cross-sectional area of the PEEK ring. By assuming no significant inductive or capacitive effects, use of a resistive heater allows for the following equation to be used. R heater is the resistance of the heater and V heater is the voltage across the heater's leads.
The contact resistances are estimated using data from [42] along with the following equation: where r cont is the value as provided by [42] and A cont is the simple geometric area of the contact surface. The thermal conductivity of the sample can then be determined by applying Equation (7), Equation (8), and Equation (9) using the data summarized in Table 2.

13
In order to compare results to previous literature, it is necessary to first define tortuosity, τ. This study uses the definition put forth in [43]: τ = L e L (10) where L e is the diffusion path length and L is the Euclidean distance of the diffusion path.
Corresponding to this definition of tortuosity, effective diffusivity, D eff , is defined as: where D is diffusivity and ε is porosity [43]. The water vapor flux, J, is calculated using: from [5] where R s is the specific gas constant of water, P sat is saturation pressure, and T is temperature in Kelvin. [5] also provides a means of calculating diffusivity: n P o P (13) where D o is diffusivity at reference absolute temperature, T o , and reference pressure, P o ; n is a fitting parameter; and P is pressure. By combining Equation (11), Equation (12), and 14 Equation (13); maintaining the assumption of one-dimensional heat flow; and discretizing differentials, water vapor flux can be expressed as: (14) where M is the molar mass of water, R u is the universal gas constant, P sat,air,top and P sat,air,bot are saturation pressures in air at the top and bottom of the sample respectively, and temperatures are in Kelvin. According to [44], saturation pressure, over the temperature range of 0 °C to 100 °C, is given by:  (15) where f is the enhancement factor (~1), P sat,vap is saturation pressure in water vapor, and temperature is in Celsius. Various -experiment data in the sample. Since water's thermal conductivity is about 20 times higher than that of air, net thermal conductivity is slightly higher at low heat flux values due to increased saturation. This observation is consistent with previous studies where higher thermal conductivities were observed for saturated samples due to water having a higher thermal conductivity than air and providing better fiber to fiber connectivity [45]. The reason such a small decrease in thermal conductivity was observed in this study is that, at all points, saturation levels were relatively low (< 0.2 compared to 0.4 -0.7).   Figure 4a shows that this occurs even when the heat flux is kept constant. The simplest explanation for this is that the system is in a pseudo-steady state. Figure 5b shows the temperature drop over the sample. The temperature drop remains relatively constant, around 8 °C. Unlike the mean temperature, the temperature drop for a given heat flux fluctuates. More specifically, the time data (Figure 5b) shows that the temperature drop tends to increase directly after an increase in heat flux and then decrease with time. At first, this seems like random fluctuations. Further consideration presents the possibility that this may be the combined result of a thermal "response time" and evaporation/condensation. Figure 5c shows near-zero saturation approaching the moment of water injection. Since evaporation/condensation is the proposed mechanism of heat transport, it is reasonable to assume that low saturation would hinder heat transport and thus prevent equalization of the sample's face temperatures. Consistent with this proposed mechanism, the injection of water at 187 minutes corresponds to a significant increase in saturation   Figure 6b shows each saturation data set normalized to its own volume average (see "Supplemental Information" section 5). Essentially, this adjusts each data set such that they may be compared as if taken at the same overall saturation, i.e., the liquid-water content remains unchanged for all heat flux and time series. If liquid-water content stays the same within the enclosed volume, it can only spatially change and redistribute. This water redistribution is clearly seen in Figure 6b where significant evaporation and a decrease in saturation is observed with time at the GDL's hot side and condensation of water is observed 305 μm into the GDL, close to the cold piston. Since it appears that most condensation has occurred by a distance of 305 μm into the sample, it is reasonable to expect that capillary-driven flow would be the dominant mechanism for water transport from 305 μm to the bottom of the sample (~445 μm). The saturation in the middle of the GDL remains unchanged. This behavior is consistent with what would be expected from PCI flow. Figure 6 Goes Here Figure 7 shows volume-rendered images of the GDLs corresponding to the data in Figure 6. These images aid in visualizing the changes described by numeric data. The time sequence is described from the left (t = 53) to the right (t = 157). The dramatic change between the last two time steps corresponds to the "rapid" evaporation shown in Figure 5c. The shift towards smaller water clusters, as will be explained by Figure 8a, is also noticeable. Furthermore, looking at the gray-scale cross-section tomographs of the GDL near the hot location (Figure 7c), fast evaporation, primarily that of large water domains, is observed. This is perhaps due to larger surface area of water droplets, which is exposed as an evaporation front. From the cross-section tomographs near the cold location of the GDL (Figure 7d), water redistribution is observed with almost no evaporation in the first three images and complete dryout in the last image.

Figure 7 Goes Here
Through-plane saturation data is useful for determining the amount of water present. However, it does not provide any insight with regards to cluster geometry. For this, it is desirable to know both the size and the shape of the water clusters as well as compare the cluster geometry to that of the pores.
Precisely defining the geometry of each pore/cluster is unrealistic and does not add much value. As such, sphere and ellipsoid fitting (described in "Supplemental Information" section 3.3) may be used to group clusters by similar geometric properties. The collection of these groups can then be represented as a discrete probability density function (PDF) relative to the total pore/cluster volume. Figure 8 shows the size and EF distributions corresponding to the data in Figure 6 (see Figure S4 for PDFs as functions of both size and EF). Plotting the distributions with respect to time visually demonstrates the impact of evaporation. As seen from the size distributions, the water clusters start with a distribution similar to that of the pores. However, as evaporation occurs, water distributions shift towards smaller cluster sizes.
This can be due to evaporation of larger clusters or breaking of large clusters into smaller ones. Also note that the initial size distribution is already biased towards smaller clusters. For radii smaller than 20 μm, water clusters favor the smallest size possible while pores have peak radii around 10 μm and no observable peak at a lower radius. Although this is visible in Figure 8a, it is more clearly shown in Figure S3c. During evaporation, this preference for forming the smallest clusters possible persists in addition to an overall distribution shift towards smaller radii (from 19.9 μm to 14.5 μm at 0 minutes and 90 minutes, respectively).
For cluster shape, the EF distributions show a preference for negative EF values. This corresponds to oblate water clusters. In contrast, the pore shapes show a preference for prolate ellipsoids. In either case, the mean EF values are close to zero. A trend worth noting is that the magnitude of the mean EF values for water clusters decrease with time (evaporation). This indicates that water cluster shapes become less oblate as evaporation occurs, even though the slight tendency towards oblate shapes remains.  (Figure 9a), it can be seen that the top (hot) region of the GDL contains noticeably less water than the middle and bottom (cold) regions. Another interesting trend occurs prior to the significant decrease in saturation. In the time between 53 minutes and 143 minutes, which is the "slow" evaporation regime, the top experiences the highest evaporation rate while the bottom experiences the lowest evaporation rate. This is due to the fact that the top is at a higher temperature than the bottom of the GDL and confirms the PCI-flow observation. Figure 9b shows that the behavior of the mean water cluster radius is similar to that of the average saturation. This is what one would expect given that Figure   6a shows a decrease in saturation with time and Figure 8a shows a decrease in mean pore radius with time. Figure 9 Goes Here The saturation data in Figure 9a shows the same type of behavior as the overall saturation data ( Figure 5c); primarily, there is a sudden increase in the saturation reduction rate at 143 minutes. As mentioned while discussing Figure 5c, this suggests two regimes of evaporation. However, identifying the first regime (53 minutes to 143 minutes) as the slow regime is actually incorrect. To understand why, it is important to refer back to the experimental setup. Recall that additional tubing was attached to the water injection tube. This additional tubing was then filled with water to achieve a specific water column height and, in doing so, control capillary pressure. Also, the tomographic scans were taken at each heat flux until the water in the GDL was completely evaporated. This, as observations confirm, means that the entire contents of the water column has also been evaporated as it entered the GDL. The reason is that, while water in the GDL is evaporated, the water column provides a reservoir from which replacement liquid water is taken. Therefore, two regimes do exist; one in which evaporated water is replaced due to the reservoir and another in which it is not because the reservoir has been depleted. These two cases are depicted in Figure 9c-d respectively. Since liquid water is replaced during the first regime, this decouples the evaporation and saturation reduction rates. In order to determine the evaporation rate, the rates of saturation reduction and reduction in water column height must be combined. Table 4 shows the calculated evaporation rates assuming that the water column becomes depleted at 143 minutes.
From this data, the first regime is the one that experiences a higher (2 orders of magnitude) evaporation rate. This is due to large saturation and evaporating surface area of water, as well as a substantial amount of water in contact with the hot piston, thus evaporating at faster rates. Lower evaporation rates in the second regime are due to much lower saturation values and evaporating surface areas of water; moreover at these low saturations water is in contact primarily with the cold piston. Given that liquid water is produced during PEFC operation, and at constant current density is injected into the GDL, the first evaporation regime is an accurate reflection of operating conditions. However, the second regime can be applicable during the cell purge at shutdown, where irreducible water saturation needs to be removed from the cell, there is no water production, and the gas phase is no longer at 100% RH. Figure 10 shows the evaporation rates predicted by Equation (14)   From area-averaged saturation plots in the through-plane GDL direction, it is evident that water near the top of the GDL (hot location) evaporates at a faster rate compared to water at the bottom of the GDL (cold location). Consideration of water cluster size and shape reveal how the geometry of the water clusters change during evaporation. The initial water size distribution closely follows that of the pore size distribution. As time progresses and more water clusters evaporate, the mean radius of the water clusters decreases from 19.9 μm to 14.5 μm; more than 5 μm. Furthermore, a large peak in water distribution is observed for very small water cluster sizes. From the ellipsoid factor (EF) data, voids have a slight tendency to be prolate while water clusters have oblate shapes.
The saturation data for the sample as a whole shows a point at which a dramatic increase in the saturation reduction rate occurs. This suggests that there are two regimes of evaporative water transport; one representing PEFC operation and the other potentially representing cell purge. The first regime actually experiences a higher evaporation rate even if saturation decrease is slow. This is due to the presence of a water reservoir during the first regime. Because water is produced during PEFC operation, it is expected that this first regime reflects operating PEFC conditions. Net water vapor flux per kg of water was calculated and it agrees well with previous experimental studies. The follow-up study will explore PCI flow under a larger window of temperatures and temperature gradients and with various GDLs.       36 Figure 9. a) Volume-average saturation with a fit (blue line) indicating evaporation when water reservoir was connected to the sample and a fit (red line) indicating evaporation without water reservoir. b) Mean water cluster size (radius) as a function of time. Data is split into thirds by equal spacing along the through-plane direction. Dashed lines indicate times at which more water was introduced. Marker shapes group data by heat flux (see Figure 4). c-d) Show (not to scale) the two different evaporation 37 cases; with and without a water reservoir respectively. e-f) Gray-scale images obtained from reconstructing data for dry and wet samples respectively. These images are labeled to show the top, middle, and bottom regions referred to in a-b. Fibers, PTFE, and water are light gray while pores are dark gray. Figure 10. Water vapor flux (left axis) and equivalent current density (right axis) for the entire GDL, as predicted by Equation (14). Circles correspond to the tortuosity-porosity value provided by [5], triangles correspond to the tortuosity and porosity (ε 1 ) values provided by [15], and squares correspond to the tortuosity provided by [15] combined with the porosity (ε 2 ) values determined in this experiment.

Symbols
Dashed lines indicate times at which more water was introduced.