Mass Transport Limitations of Water Evaporation in Polymer Electrolyte Fuel Cell Gas Diffusion Layers

Facilitating the proper handling of water is one of the main challenges to overcome when trying to improve fuel cell performance. Specifically, enhanced removal of liquid water from the porous gas diffusion layers (GDLs) holds a lot of potential, but has proven to be non-trivial. A main contributor to this removal process is the gaseous transport of water following evaporation inside the GDL or catalyst layer domain. Vapor transport is desired over liquid removal, as the liquid water takes up pore space otherwise available for reactant gas supply to the catalytically active sites and opens up the possibility to remove the waste heat of the cell by evaporative cooling concepts. To better understand evaporative water removal from fuel cells and facilitate the evaporative cooling concept developed at the Paul Scherrer Institute, the effect of gas speed (0.5–10 m/s), temperature (30–60 °C), and evaporation domain (0.8–10 mm) on the evaporation rate of water from a GDL (TGP-H-120, 10 wt% PTFE) has been investigated using an ex situ approach, combined with X-ray tomographic microscopy. An along-the-channel model showed good agreement with the measured values and was used to extrapolate the differential approach to larger domains and to investigate parameter variations that were not covered experimentally.


Introduction
During polymer electrolyte fuel cell (PEFC) operation, water is generated as a byproduct of the electrochemical reaction. This water, besides being absorbed by the proton conductive membrane and evaporating on the anode [1,2], commonly leaves the PEFC cathode either in liquid or gaseous form through the gas channels in the cathode side of the bipolar plates [3,4]. The transport path for either gas or liquid water removal is the same, with both having to progress from the cathode catalyst layer (CL) through a fibrous gas diffusion layer (GDL) to the gas channels. In PEFCs, besides distributing the reactant gases, the GDL provides electrical conductivity between the CL and the bipolar plates, conducts heat away from the CL, and serves as mechanical support to the membrane [5][6][7][8]. In the case of gaseous removal of water, evaporation might occur in the CL or the GDL structure and is followed by Fickian diffusion through the GDL's pore spaces to the gas channel [9]. For liquid removal, product water accumulations in the GDL eventually grow so large in size that they break through into the gas channel where the formed droplets are carried away by the gas flow [10][11][12][13][14][15]. This mechanism of water accumulation and subsequent growth of clusters has the undesired side effect of blocking the pore space, which should be avoided in order to maintain undisturbed supply of reactant gas to the active sites in the CL [16,17]. It is thus desirable to transport as much water as possible out of the cell in vapor form, minimizing the accumulation of liquid water.
A novel cooling strategy developed at the Paul Scherrer Institute aims to use patterned wettability GDLs to remove heat through the phase change from liquid to vapor [18][19][20]. The artificial distribution of hydrophobic and hydrophilic GDL domains enables tailormade patterns of water filled lines with high levels of evaporation rates and necessitates a good understanding of the vapor removal mechanisms, ultimately limiting the rate of evaporation. Understanding the characteristic mechanisms and limitations to evaporative water removal in the framework of PEFCs and quantifying them has been the aim of many scientific works, both experimentally [21][22][23][24][25] and numerically [26][27][28][29][30][31][32]. However, the reported evaporation values show a wide spread [31,33] and some measurements performed on differential cell scales leave open questions of how evaporative flux scales with GDL thickness and gas type are used in the evaporation rate measurements, as reported by Lal et al. [21]. They observed that doubling the diffusive distance did not halve the diffusive flux, as was expected under differential conditions, while the evaporation increase due to changing from nitrogen to hydrogen in the gas channel fell short of the theoretical effect the increased diffusivity should have. Boundary effects, caused by the differential cell approach with very small evaporative domains, could be a potential origin of these inconsistencies. They can lead to large deviations when upscaling to larger domains or even make it impossible to do so, also because accumulations of water vapor in the gas stream can become decisive or invalidate the assumption of differential operation.
Hence, in this study, we perform direct evaporation rate measurements from a differential ex situ cell, imitating a fuel cell GDL and gas channel. To control all relevant boundary conditions during the experiment and ensure comparability between multiple runs, liquid water was prevented from entering the GDL domain, resulting in pure diffusive transport and avoiding sample specific water cluster formations. The goal of this study is to quantify evaporative water removal for differently sized evaporative domains and help unify the available literature results. In addition to the presented experimental results, a numerical 2D model was applied. It accounts for diffusion in the GDL domain and combined convective and diffusive transport in the gas channel to provide further insights into the governing processes and phenomena. Spatially resolved information of the distribution of the water vapor concentration profiles deepens the understanding of water evaporation and helps explain the experimentally observed trends. Moreover, the modelling work allows extrapolations to larger domains and to compare the differential cell results to evaporation rates determined in larger cells.

Materials and Methods
The experiments were designed to emulate evaporative water removal from a saturated catalyst layer through the GDL into a gas channel and investigate the effects of temperature, gas flow, and evaporation domain size on the evaporation rates. An ex situ cell design [11,34] (Figure 1a) was chosen to assure controllability of the boundary conditions and compatibility with X-ray tomographic microscopy (XTM), which is commonly used to determine the distribution of water inside fuel cell GDLs [22,23,[34][35][36][37][38][39]. The configuration of water and GDL during the evaporation measurement can be seen in the XTM images in Figure 1b. (a) Experimental setup consisting of syringe, injection cell, and laser levels (L1 and L2), as well as the cell used to measure evaporation rates. The close-up of the injection cell shows the flow field with a single gas channel (gas path highlighted in green), the injector used to guide the water to the base of the gas diffusion layer (GDL), the pool gasket with an exemplary evaporation domain geometry, and the seal used to prevent leakage and control the compression of the GDL. (b) Through-plane X-ray tomographic microscopy (XTM) image along the gas channel of the cell showing the arrangement of the cell components as well as the waterfront established under the GDL. The close-up of the injection cell shows the flow field with a single gas channel (gas path highlighted in green), the injector used to guide the water to the base of the GDL, the pool gasket with an exemplary evaporation domain geometry, and the seal used to prevent leakage and control the compression of the GDL. (b) Reconstructed tomography images perpendicular to and along the gas channel of the cell, showing the arrangement of the cell components as well as the waterfront maintained under the GDL.

Materials and Setup
The GDL, a TGP-H 120 10 wt% PTFE (Toray Industries, Inc., Tokyo, Japan) with a nominal thickness of 370 µm and an average porosity of 78%, was placed in the cell, consisting of a flow field and an injector, both made from graphite (BMA5, Eisenhuth GmbH as well as the cell used to measure evaporation rates. The close-up of the injection cell shows the flow field with a single gas channel (gas path highlighted in green), the injector used to guide the water to the base of the gas diffusion layer (GDL), the pool gasket with an exemplary evaporation domain geometry, and the seal used to prevent leakage and control the compression of the GDL. (b) Through-plane X-ray tomographic microscopy (XTM) image along the gas channel of the cell showing the arrangement of the cell components as well as the waterfront established under the GDL. The close-up of the injection cell shows the flow field with a single gas channel (gas path highlighted in green), the injector used to guide the water to the base of the GDL, the pool gasket with an exemplary evaporation domain geometry, and the seal used to prevent leakage and control the compression of the GDL. (b) Reconstructed tomography images perpendicular to and along the gas channel of the cell, showing the arrangement of the cell components as well as the waterfront maintained under the GDL.

Materials and Setup
The GDL, a TGP-H 120 10 wt% PTFE (Toray Industries, Inc., Tokyo, Japan) with a nominal thickness of 370 µm and an average porosity of 78%, was placed in the cell, consisting of a flow field and an injector, both made from graphite (BMA5, Eisenhuth GmbH & Co. KG, Osterode am Harz, Germany). This GDL was chosen to keep overall evaporative fluxes moderate to avoid thermal gradients as much as possible. The onset of thermal limitations was observed for thinner GDLs at 80 • C in a similar setup [21]. A seal was placed around the GDL, preventing leakage of vapor and controlling the GDL compression. The seal is made up of a 160 µm PTFE and a 25 µm FEP foil. In combination with a 125 µm groove in the flow field, this resulted in a compressed GDL thickness of 310 µm and a compression of 16%. The flow field simulates a gas channel (width 800 µm, depth 300 µm) of a fuel cell, while the injector is used to supply the water to the base of the GDL via a small hole. Both the injector and the flow field can be heated using resistive heating pouches and their temperature monitored with temperature sensors (PT100).
To define the evaporation domain, a hydrophobic FEP foil with a rectangular cutout was used. The foil contained water in the rectangular cutout and prevented uncontrolled spreading of the water underneath the GDL. This way, a well-defined evaporation domain was created with a predictable geometrical water surface area. The evaporation domains realized in this study all had a width of 1.6 mm and varying length along the channel between 0.8 and 10 mm. Note that, to change the size of the evaporation domain, a new cell had to be assembled with a new GDL, seal, and gaskets.
The evaporation domain was filled with deionized water using an external syringe pump (Legato ® 110, KD Scientific Inc., Holliston, MA, USA) to a pressure p w , which was controlled using a laser level sensor (L1) placed on a riser column. Every centimeter of water in the column above the evaporation domain in the cell resulted in an increase in water pressure of 1 mbar. For the experiments in this study, L1 was set to 5 mbar (5 cm) to obtain good filling of the evaporation domain without penetration of water into the bulk of the GDL. Preventing the water from advancing into the GDL ensured that the geometrical distance from the waterfront to the gas channel was even, reproducible, and defined only by the compressed GDL thickness.
X-ray tomographic microscopy in a CT scanner (Nanotom, GE Sensing & Inspection, Wunstorf, Germany) was employed to verify proper filling of the evaporation domain with water. Additionally, it enabled the validation of the level of compression of the GDLs in the cells and confirmed that the water did not extend into the GDL domain. An example of this can be seen in Figure 1b, showing the water in the evaporation domain of the pool gasket and following the contour of the GDL without penetration. The average vertical distance from the waterfront to the gas channel was confirmed using this method and found to be 310 µm.
The gas flowing through the channel in the flow field was dry nitrogen at the respective cell temperature with velocities between 0.5 m/s and 10 m/s. These velocities span flow rates commonly encountered in air fed cathode gas channels of large cells as well as high velocities of interest for evaporative cooling. The flow rate was controlled by a mass flow controller (F-201CV-500-AAD-33-V, Bronkhorst, AK Ruurlo, The Netherlands) at 25 • C in accordance with the experimental parameters in Table 1. For comparability, all flow rates reported in this study refer to the dry gas flow rate at 25 • C, before being heated and guided into the cell.

Measurement Procedure
At the start of each measurement, water was injected into the system until the water level in the riser reached L1 (Figure 1a), at which point the injection was stopped. Owing to the presence of evaporation in the cell, the water volume in the system then started to decrease, resulting in a drop of the water level in the riser. The time (t L1−L2 ) it took for the water in the riser column to drop to a second laser level L2, which was located one centimeter below L1, was then measured. With the measured time and the known volume of water stored between L1 and L2 (∆V L1−L2 = 1.96 µL), the evaporation rate q was determined. Note again, that to ensure a constant evaporation front in the cell, the laser levels L1 and L2 were both positioned a few centimeters above the cell, resulting in a slight overpressure in the water of~500 Pa. This pressure, while not enough to push any water into the GDL domain, acts as a buffer to prevent movement of the water front in the cell during evaporation rate determinations and ensured that the change in water volume manifested in the riser exclusively. All measured evaporation rates were corrected for a system intrinsic loss rate (1.33 nL/s), which originated from diffusion through the tubes, imperfections in the fittings, as well as the riser column being open to the ambient atmosphere at the top. The latter was kept low by the riser column extending roughly one meter beyond L1, limiting the diffusion of water from this location while keeping the system open to ambient pressure at the water air interface in the riser. This type of measurement was repeated three times for the temperatures and gas speeds investigated. A full list of the experimental parameter variations is available in Table 1.

Numeric Analysis
A two-dimensional, isothermal, macro-homogeneous along-the-channel model was implemented in MATLAB ® (R2019b, © The MathWorks, Inc., Natick, MA, USA), taking into account Fickian diffusion in both the GDL and gas channels, as well as convection in the gas channel. Figure 2 gives an overview of the two-dimensional modelling domain along the center line of the gas channel. Note that, with this 2D modeling approach, contributions from the third dimension (perpendicular to the gas channel) are not captured. The simulation domain also includes the inlet and outlet domains of 1.3 mm in addition to the evaporative domain length (L). the GDL domain, acts as a buffer to prevent movement of the water front in the cell during evaporation rate determinations and ensured that the change in water volume manifested in the riser exclusively. All measured evaporation rates were corrected for a system intrinsic loss rate (1.33 nL/s), which originated from diffusion through the tubes, imperfections in the fittings, as well as the riser column being open to the ambient atmosphere at the top. The latter was kept low by the riser column extending roughly one meter beyond L1, limiting the diffusion of water from this location while keeping the system open to ambient pressure at the water air interface in the riser. This type of measurement was repeated three times for the temperatures and gas speeds investigated. A full list of the experimental parameter variations is available in Table 1.

Numeric Analysis
A two-dimensional, isothermal, macro-homogeneous along-the-channel model was implemented in MATLAB® (R2019b, © The MathWorks, Inc., Natick, MA, USA), taking into account Fickian diffusion in both the GDL and gas channels, as well as convection in the gas channel. Figure 2 gives an overview of the two-dimensional modelling domain along the center line of the gas channel. Note that, with this 2D modeling approach, contributions from the third dimension (perpendicular to the gas channel) are not captured. The simulation domain also includes the inlet and outlet domains of 1.3 mm in addition to the evaporative domain length (L). While a more in-depth description of the model implementation is included in the Appendix A, the main parameters are given in Table 2.  While a more in-depth description of the model implementation is included in the Appendix A, the main parameters are given in Table 2.

Results and Discussion
As the experiment was designed to only evaporate water from a predefined geometry underneath the GDL without liquid water penetration into the GDL domain, water transport to the gas channel occurs primarily in the form of diffusion. Once in the gas channel, the gas stream carries the vapor out of the cell. The dependence of the evaporation rate q on the gas flow rate v Gas,in can be seen in Figure 3a for a temperature of 60 • C.  1 (0.176 + 0.0023 T) × 10 -4 m 2 /s Diffusivity (H2O in H2) 1 (0.8912 + 0.0039 T) × 10 -4 m 2 /s TP coefficient ( ) 2 0.28 -IP coefficient ( ) 3 1.9 × -

Results and Discussion
As the experiment was designed to only evaporate water from a predefined geometry underneath the GDL without liquid water penetration into the GDL domain, water transport to the gas channel occurs primarily in the form of diffusion. Once in the gas channel, the gas stream carries the vapor out of the cell. The dependence of the evaporation rate q on the gas flow rate , can be seen in Figure 3a for a temperature of 60 °C.  Especially for the longer evaporation domains, the evaporation rate first increases drastically with increasing gas flow rates, following the limit of maximum moisture saturation in the gas channel as indicated by the dotted line. For shorter domains, this effect is not so pronounced, as almost all investigated gas flow rates offer sufficient convective vapor removal to maintain low saturation in the gas channel. For lower velocities, a decrease of the evaporation rate was observed, indicating a significant increase of saturation in the carrier gas over the evaporation domain. At higher gas flow rates, for all L, a plateau was reached, after which an increase of gas flow rate showed only a minor or no increase in evaporation. This indicates that the mechanism limiting vapor transfer has shifted away from the gas channel to the GDL domain, and diffusion driven transport is now dominating the evaporation characteristics. It is assumed that, under the given conditions, diffusion is the only transport mode in the GDL. The vertical error bars represent the maximum error observed when propagating the measurement uncertainties of ∆t L1−L2 and ∆V L1−L2 .
The values obtained from the model show very good agreement with the measured values, indicating that the effects taken into consideration describe the vapor transport processes under the investigated conditions accurately enough to allow extrapolations exceeding the experimentally investigated parameter range.
For different temperatures, at L = 10 mm, the effect of gas flow rate is shown in Figure 3b. An initial increase at low gas velocities can be seen as well as a plateauing at higher gas velocities. The general behavior appears very similar in terms of the average velocity required to reach the plateau being comparable for all temperatures. It was found that, at sufficiently high velocities, the plateau values of the different temperatures relate to each other in accordance with their respective saturation pressures: where q is the evaporation rate at a respective temperature T, D the diffusivity, p sat is the water saturation pressure, and v Gas is the average gas velocity. This highlights again the transition from convective (saturation) limiting conditions at low velocities to diffusion limited vapor transport at higher velocities. As the evaporation domain did not extend underneath the entire GDL along the gas channel, the diffusive transport can be expected to have a different behavior at the edges of the domain than over the bulk. The contribution of this edge region to the overall evaporation rate as function of L was quantified, considering only results for gas velocities v Gas,in ≥ 4 m/s to ensure diffusion limited conditions (Figure 4a).
It was found that the evaporation rates follow a linear increase with L for all temperatures indicated by the linear fits (dashed lines, see Figure 4a). These fits have a non-zero ordinate value that can be interpreted as the edge contribution to the overall evaporation rates and enables the quantification of the two diffusion domains: (i) Through-plane (TP) diffusion over the bulk of the evaporation domain; (ii) In-plane (IP) and through-plane (TP) diffusion on the up and downstream edges of the evaporation domain.
These are indicated schematically in Figure 4b with the bulk and edge contributions in green and red, respectively. Note that this illustration aims just to indicate the locations of these regions, not to define their exact dimensions. This interpretation holds only under differential conditions, where a change in L only affects the contribution of the TP diffusion component to the overall evaporation rate, while the edge contributions remain the same. This effect could be very well observed in this study when varying L between 0.8 and 10 mm, as shown in Figure 4a for high gas velocities (≥4 m/s). The experimentally obtained  Table 3. These slope values represent the upper limit of diffusive vapor transport through the GDL domain under optimal, purely diffusion limited conditions at the respective temperatures. differential conditions, where a change in L only affects the contribution of the TP diffusion component to the overall evaporation rate, while the edge contributions remain the same. This effect could be very well observed in this study when varying L between 0.8 and 10 mm, as shown in Figure 4a for high gas velocities (≥4 m/s). The experimentally obtained values for both the slopes of the linear fits and their ordinate values are reported in Table 3. These slope values represent the upper limit of diffusive vapor transport through the GDL domain under optimal, purely diffusion limited conditions at the respective temperatures.  The same averaging and linear fitting has been applied to the evaporation rates obtained by the model at same values for L as used in the measurements and the resulting  The same averaging and linear fitting has been applied to the evaporation rates obtained by the model at same values for L as used in the measurements and the resulting values for the slopes and ordinates are included in Table 3. The slope values show good agreement, while the ordinate values differ slightly. From this, it can be seen that, at 30 • C, an evaporation domain length of L = 0.5 mm contributes almost equally as much to the total evaporation as its edge domain. At 60 • C, this contribution is even larger, comparable to an Energies 2021, 14, 2967 9 of 21 evaporation domain of~1 mm in length. With increasing L, the contribution of this edge domain decreases, reducing to 10% at~5 mm and~8.5 mm at 30 • C and 60 • C, respectively.
More detailed insight into the edge effects is obtained when looking at the spatial distribution of water vapor inside the GDL and the gas channels. Figure 5 shows the relative humidity (RH) obtained from the model for L = 0.8 and 2.4 mm as well as for v Gas,in = 0.5 and 6 m/s at a temperature of 60 • C inside the modeling domain in accordance with Figure 2. In this representation, the difference in water concentration between low and high gas velocities over the length of the evaporation domain becomes apparent. While for v Gas,in = 0.5 m/s (Figure 5a,c), the relative humidity in the gas channel increases substantially over the evaporation domain (up to~20 and~40% for L = 0.8 mm and 2.4 mm, respectively), it remains low (<10%) in the case of v Gas,in = 6 m/s (Figure 5b,d). In the GDL, clear differences can be observed especially around the downstream edge of the evaporation domain. For the case of low gas velocities, the humidity gradients between the bottom of the GDL (line D) and the GDL/gas channel interface (line B) become smaller towards the end of the simulation domain, while they remain almost the same for the case of high gas velocity.
values for the slopes and ordinates are included in Table 3. The slope values show good agreement, while the ordinate values differ slightly. From this, it can be seen that, at 30 °C, an evaporation domain length of L = 0.5 mm contributes almost equally as much to the total evaporation as its edge domain. At 60 °C, this contribution is even larger, comparable to an evaporation domain of ~1 mm in length. With increasing L, the contribution of this edge domain decreases, reducing to 10% at ~5 mm and ~8.5 mm at 30 °C and 60 °C, respectively.
More detailed insight into the edge effects is obtained when looking at the spatial distribution of water vapor inside the GDL and the gas channels. Figure 5 shows the relative humidity (RH) obtained from the model for L = 0.8 and 2.4 mm as well as for , = 0.5 and 6 m/s at a temperature of 60 °C inside the modeling domain in accordance with Figure 2. In this representation, the difference in water concentration between low and high gas velocities over the length of the evaporation domain becomes apparent. While for , = 0.5 m/s (Figure 5a,c), the relative humidity in the gas channel increases substantially over the evaporation domain (up to ~20 and ~40% for L = 0.8 mm and 2.4 mm, respectively), it remains low (<10%) in the case of , = 6 m/s (Figure 5b,d). In the GDL, clear differences can be observed especially around the downstream edge of the evaporation domain. For the case of low gas velocities, the humidity gradients between the bottom of the GDL (line D) and the GDL/gas channel interface (line B) become smaller towards the end of the simulation domain, while they remain almost the same for the case of high gas velocity. For a more detailed comparison, the relative humidity values shown in Figure 5 were extracted at lines A-D for , = 0.5 m/s and 6 m/s, respectively. It can be seen, especially For a more detailed comparison, the relative humidity values shown in Figure 5 were extracted at lines A-D for v Gas,in = 0.5 m/s and 6 m/s, respectively. It can be seen, especially for L = 2.4 mm, that the humidity extracted at the GDL/gas channel interface (line B) shows a distinct increase along the evaporation domain (starting at x = 1.3 mm) along the x-axis, thus reducing the gradients through the GDL, resulting in reduced evaporation along the channel ( Figure 6). Contrarily, for v Gas,in = 6 m/s (Figure 6b), the humidity in the gas Energies 2021, 14, 2967 10 of 21 channel remains rather constant over the evaporation domain (line A). A slight humidity increase can be observed at the GDL/gas channel interface (line B), which drops back down to~10% downstream of the evaporation domain as the vapor distributes in the gas channel. With a decreasing humidity gradient, the diffusive flux decreases as well, resulting in a transition away from diffusion-limited to convection limited mass transport when the RH in the channel is dominating. It is apparent that, owing to this saturation increase, the evaporation rate will decrease along the channel until, eventually (for a large enough evaporation domain), the relative humidity in the gas channel will reach unity and no more evaporation can occur. The closer the saturation in the gas channel gets to RH = 1, the more the overall evaporation rate depends on the convective removal via the gas channel. for L = 2.4 mm, that the humidity extracted at the GDL/gas channel interface (line B) shows a distinct increase along the evaporation domain (starting at x = 1.3 mm) along the x-axis, thus reducing the gradients through the GDL, resulting in reduced evaporation along the channel ( Figure 6). Contrarily, for , = 6 m/s (Figure 6b), the humidity in the gas channel remains rather constant over the evaporation domain (line A). A slight humidity increase can be observed at the GDL/gas channel interface (line B), which drops back down to ~10% downstream of the evaporation domain as the vapor distributes in the gas channel. With a decreasing humidity gradient, the diffusive flux decreases as well, resulting in a transition away from diffusion-limited to convection limited mass transport when the RH in the channel is dominating. It is apparent that, owing to this saturation increase, the evaporation rate will decrease along the channel until, eventually (for a large enough evaporation domain), the relative humidity in the gas channel will reach unity and no more evaporation can occur. The closer the saturation in the gas channel gets to RH = 1, the more the overall evaporation rate depends on the convective removal via the gas channel. The model was used to extrapolate to larger evaporation domains and higher temperatures that were not covered experimentally. The reader should keep in mind that the resulting evaporation rates and humidity distributions only hold true for evaporation experiments or cell designs with dedicated liquid feed for evaporative cooling. In conventional PEFCs, the evaporation rates at static conditions will be limited additionally by the local current density and the respective water production. Figure 7a shows the evaporation rates at 60 °C for inlet gas velocities between 0.5 and 10 m/s and L up to 200 mm. The model was used to extrapolate to larger evaporation domains and higher temperatures that were not covered experimentally. The reader should keep in mind that the resulting evaporation rates and humidity distributions only hold true for evaporation experiments or cell designs with dedicated liquid feed for evaporative cooling. In conventional PEFCs, the evaporation rates at static conditions will be limited additionally by the local current density and the respective water production. Figure 7a shows the evaporation rates at 60 • C for inlet gas velocities between 0.5 and 10 m/s and L up to 200 mm.
As can be seen, the evaporation rates increase with the length of the evaporation domain as well as the gas velocity. As L increases, the evaporation rates reach a limit; for lower gas velocities, this limit is reached after shorter L compared with higher gas velocities. The dotted lines indicate 25% increments in average relative humidity in the gas channel outlet. The lines indicating 20-100% were extracted by determining the humidity in the gas channel outlet at the respective conditions and extrapolated for visibility. The~0% line is tangential to the 10 m/s curve at L → 0 mm, where the average humidity in the gas channel was closest to 0%. Initially, for small values of L and thus close to 0% average humidity in the channel outlet, the slope of q with L coincides with the value obtained at differential conditions ( Table 3). As L gets larger, the transition from diffusion to saturation limit vapor transport can be observed. For the lower gas velocities, this occurs already at low values of L, while for the higher gas velocities, L can be larger until, even at v Gas,in =10 m/s, a significant deviation from pure diffusion limitation can be observed around L = 20 mm as saturation in the gas channel starts to increase. It can be seen that, at an average humidity of 25% in the gas channel, the diffusive transport, while no longer under differential conditions, still plays a dominant role in the vapor transport mechanism. As the average gas channel humidity increases, the convective transport limitation starts to dominate more and more while local evaporation reduces. As an example, for v Gas,in = 2 m/s, 75% average humidity in the gas channel is already reached at L = 40 mm, where after it would require another 160 mm of evaporation domain to reach full saturation at L = 200 mm. As can be seen, the evaporation rates increase with the length of the evaporation domain as well as the gas velocity. As L increases, the evaporation rates reach a limit; for lower gas velocities, this limit is reached after shorter L compared with higher gas veloc-  The temperature dependence of the evaporation rates for L up to 200 mm is shown in Figure 7b at a constant gas inlet velocity of 6 m/s. With increasing L and temperature, the evaporation rates increase. For low values of L, the evaporation rates follow different initial slopes as diffusive transport is limiting in this regime and, with increasing temperature, the diffusivity as well as the concentration gradient through the GDL (driven by the vapor pressure at the water interface) increase. As the average humidity in the gas channel increases for longer evaporation domains, we see a deviation from the respective initial slopes. Overall, it can also be observed that, with increasing temperature, the gas channel saturates faster, reducing the evaporation domain length required for 75% humidity in the channel from L = 130 mm at 30 • C to 95 mm at 80 • C. A humidity of 100% in the gas channel is not reached in this example as the gas velocity of v Gas,in = 6 m/s requires an evaporation domain length greater than 200 mm to be fully saturated under the given conditions.
A variation of GDL thickness has been investigated numerically by reducing the simulated GDL thickness by 50% while leaving all other parameters the same to understand the influence of diffusive path length. To better visualize this effect, Figure 8a shows the results in the form of the normalized evaporation rate, q norm , as a function of L for different gas inlet velocities. This normalization is done by dividing the obtained evaporation rates at any given condition with the maximal amount of vapor removable by the fully saturated gas flow at the respective conditions. Note that the normalization to the limiting value at full gas saturation is equivalent to the calculation of the average relative humidity. In an attempt to avoid excessive overlap, only four velocities are shown. In this representation, a thinning of the GDL and the accompanying reduction in diffusive distance results in a higher rate of evaporation as indicated by the steeper increase of evaporation rate with L for thinner GDLs (Figure 8a). This also results in faster saturation of the gas in the channel and, as such, the convective limitation is reached for lower values of L. Ideally, if evaporation is to occur evenly everywhere along the channel, the gas velocity would have to be higher for thinner GDLs than for thicker GDLs.
When only considering through-plane diffusion in the GDL, according to Fick's law, the length contributes inversely to the diffusive flux: From the ratios shown in Figure 8b, we see that, in the present system, a reduction of the diffusive pathway by 50% does not result in a doubling of the evaporation rate. This is because, for short evaporation domains, where diffusion is limiting, the contribution of edge diffusion is not negligible and, because it does not scale linearly with the GDL thickness, we do not see the expected rate doubling. As the contribution of this edge diffusion to the overall vapor removal decreases with increasing L, a relative evaporation rate increase can be observed until the rising saturation in the gas channel starts to push the system towards a convective limitation, reducing the effect from the GDL thickness reduction. Thus, for any gas velocity, a GDL thickness reduction only increases the total evaporative water removal if L is small enough to avoid fully saturating the gas channel. This effect of increased evaporation rates has been observed experimentally by Lal et al. [21] who, with a comparable setup, observed that a thickness reduction from TGP-H 120 to TGP-H 060 resulted in an evaporation rate increase that was lower than expected from the reduction in diffusive path length.
While water is generated on the cathode side, evaporation can also occur on the anode. In regular fuel cell operation, this is can occur in specific cases owing to water migration through the membrane [9], but, in recent years, concepts for active water injection on the anode side have also been investigated for cooling purposes [18]. With hydrogen as gas on the anode, the diffusivity of water vapor compared with nitrogen is increased. The effect of this increased diffusivity on evaporation is shown in Figure 9a in the form of q norm at 60 • C for different gas velocities. to TGP-H 060 resulted in an evaporation rate increase that was lower than expected from the reduction in diffusive path length. While water is generated on the cathode side, evaporation can also occur on the anode. In regular fuel cell operation, this is can occur in specific cases owing to water migration through the membrane [9], but, in recent years, concepts for active water injection on the anode side have also been investigated for cooling purposes [18]. With hydrogen as gas on the anode, the diffusivity of water vapor compared with nitrogen is increased. The effect of this increased diffusivity on evaporation is shown in Figure 9a in the form of at 60 °C for different gas velocities. Similar to increasing the diffusive vapor transport by reducing the diffusive distance, changing the carrier gas from nitrogen to hydrogen results in an increase in evaporation rates due to the higher diffusivity of water vapor in hydrogen. At this temperature, the diffusivity ratio between hydrogen and nitrogen is~3.58 [40], resulting in a more drastic change of q with regards to L in hydrogen for small values of L, and thus also a faster buildup of saturation in the gas channel. Contrarily to the reduction of diffusive distance, increasing the diffusivity affects all diffusive transport equally. Thus, for very short evaporation domains, the evaporation rate increases by a factor of 3.58, as shown in Figure 9b, obtained by dividing the evaporation rate in hydrogen by the evaporation rate in nitrogen. As before, the evaporation rate gain reduces with increasing values of L as the saturation in the gas channel slows down the rate of evaporation until, eventually, the gas in the gas channel is fully saturated. For this reason, an increased diffusivity due to a change of gas type does not result in the same increase occurring in the evaporation rates, unless L is unreasonably short. Lal et al. have also observed this effect experimentally for the change from nitrogen to hydrogen as carrier gas with the measured evaporation rates falling short of what would be expected from the change in diffusivity [21]. Note that the effects of changing gas density due to the large difference in molar masses between water vapor and carrier gas are not captured by this model, as Equation (4) assumes that the total gas density stays constant. Similar to increasing the diffusive vapor transport by reducing the diffusive distance, changing the carrier gas from nitrogen to hydrogen results in an increase in evaporation rates due to the higher diffusivity of water vapor in hydrogen. At this temperature, the diffusivity ratio between hydrogen and nitrogen is ~3.58 [40], resulting in a more drastic change of q with regards to L in hydrogen for small values of L, and thus also a faster buildup of saturation in the gas channel. Contrarily to the reduction of diffusive distance, increasing the diffusivity affects all diffusive transport equally. Thus, for very short evaporation domains, the evaporation rate increases by a factor of 3.58, as shown in Figure 9b, obtained by dividing the evaporation rate in hydrogen by the evaporation rate in nitrogen. As before, the evaporation rate gain reduces with increasing values of L as the saturation in the gas channel slows down the rate of evaporation until, eventually, the gas in the gas From Figures 8 and 9, it becomes especially apparent how important it is to maintain differential conditions in the gas channel when investigating diffusion limited evaporation rates from fuel cell GDLs. This requirement becomes increasingly harder to maintain experimentally for cases with high evaporation rates (high temperatures, thin GDLs, or gases with high vapor diffusivity) and would necessitate very high gas velocities. Furthermore, evaporation rate normalizations with evaporative area have to be carefully examined for the conditions under which they have been acquired, because, for small domains, a normalization with area suggests very high rates of water removal. In reality, such high rates of evaporation do not scale to larger evaporation domains. This is exemplarily shown in Figure 10a,b, which are obtained by dividing the evaporation rates from Figure 7

Conclusions
Diffusive and convective limitations to the evaporation of water from fuel cell GDLs were investigated using an ex situ cell capable of measuring accurate evaporation rates at various fuel cell operation conditions. Furthermore, a 2D along-the-channel model was applied, enabling insights into spatially resolved humidity distributions. The model results showed good agreement with experimental data without the need for any fitting parameters, using only literature values for the structural bulk properties of the GDL together with information about the geometry of the experiment. While not the subject of the current study, it is important to note that heat transfer limitations may also play a role if the local evaporation rates become too high or thermal conductivities of the measurement setup materials are too low.

Conclusions
Diffusive and convective limitations to the evaporation of water from fuel cell GDLs were investigated using an ex situ cell capable of measuring accurate evaporation rates at various fuel cell operation conditions. Furthermore, a 2D along-the-channel model was applied, enabling insights into spatially resolved humidity distributions. The model results showed good agreement with experimental data without the need for any fitting parameters, using only literature values for the structural bulk properties of the GDL together with information about the geometry of the experiment.
Increasing the gas velocity results in an increase of the evaporation rate until eventually a plateau was reached, after which further increases of the gas velocity show almost no further evaporation rate increase. A clear transition from convection limited evaporation, at low gas velocities, to diffusion limited evaporation, at high gas velocities, was observed and the average humidity in the gas channel was identified as the ruling property for the transition. Up to an average humidity of 25% in the gas channel, diffusion through the GDL was the main limitation for the evaporation of water. With increasing saturation in the gas channel, the diffusive vapor transport through the GDL decreased and ceased once the channel saturation reached 100%, as is the case for low gas flow velocities or long domains. In general, an uneven evaporation distribution along the channel was observed where the majority of the water evaporates over a short initial distance, saturating the gas channel and decreasing evaporation further downstream.
From both the experiments and the model, the contributions of through-plane and edge diffusion in the GDL could be identified for diffusion-limited conditions. The contribution of edge diffusion can be almost 50% for a 1 mm evaporative domain, but levels of 10% at L =~5 mm and~8.5 mm at 30 • C and 60 • C, respectively. From the through-plane contribution, an upper limit for diffusive vapor transport through the GDL was determined to be 2.44 nL/s/mm at 60 • C.
Reducing the GDL thickness by a factor of 2 did not result in doubling the observed evaporation rates, which was found to be due to the contributions of edge regions of the evaporation domain as well as the accumulation of water vapor along the gas channel. With hydrogen as the carrier gas, and for very short domains, an increase in evaporation was found, which coincided with the increase of diffusivity compared with nitrogen gas. For larger evaporation domains, the saturation buildup in the gas channel reduces the diffusion flux through the GDL, lessening the impact of the increased diffusivity.
These factors highlight the importance of careful analysis of the boundary conditions under which evaporating rates are determined experimentally. Especially values normalized by area run the risk of unintentional bias, as experiments with short domains overestimate the rates of evaporation when scaled to larger domains. Contrarily, measurements performed with large domains do not resolve how unevenly the evaporation is distributed across the evaporative domain.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Acknowledgments:
The authors would like thank Yu Kakizawa for his contribution of analyzing the velocity profiles in the gas channel and the GDL using Geodict simulations. Furthermore, Thomas Gloor and Marcel Hottiger are acknowledged for their support towards setup development and software implementation.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. D and v are considered constant along the channel and stepwise constant (over one spatial increment) in the through-plane (TP) direction. As laminar flow is expected in the gas channels (Re < 200), v is assumed to show only a contribution in the x direction. Moreover, the evaporation process at the water surface is modelled as a fully saturated air boundary condition, thus S = 0 inside the computational domain. Based on these assumptions, Equation (A1) can be simplified for two-dimensional diffusion and onedimensional convection without a source term as follows:

Glossary
with D i as the diffusion coefficient and v i as the gas velocity in the respective direction i.

Appendix A.2. Model Implementation
In order to solve Equation (A2), the finite difference method was applied and implemented in MATLAB®. Equation (A3) shows the implicit discretization of the partial differential equation (Equation (A2)), using the central difference scheme for the diffusion process and the upwind scheme for the convection process (as the Peclet number is sufficiently large (Pe > |2|) [43]). The implicit implementation is unconditionally stable with respect to the time step ∆t, and thus allows for rapid and accurate calculation of the steady-state solution.
where c t x,z is the water vapor concentration at position x, z at time t and ∆ is the step size of the respective variable.
The implemented boundary conditions are depicted in Figure 2. Gas with no humidity and a defined velocity profile v x (z) enters the channel at the inlet and leaves the channel on the right side at the outlet, leading to convective fluxes across the system boundaries. In order to model the evaporation at the interface between the water and the GDL, fully saturated gas is assumed and, therefore, a Dirichlet boundary condition (full saturation at the interface) is implemented. The saturation concentration is calculated with the Antoine equation as a function of temperature and ideal gas behavior is assumed. Moreover, no-flux conditions are assigned to all other boundaries and the diffusive flux across the channel inlet and outlet is assumed to be negligible. The velocity profile is obtained by first determining v x (z, y) using the following equation from Shah et al. [44]: where H ch and W ch are the height and width of the channel, respectively, with n = 2.0125 and m = 3.6739 (according to Shah et al. [44] for the given channel geometry). This profile is then averaged over the y-direction to obtain v x (z) as depicted in Figure A1 (solid lines).
To assess the validity of this approach, the obtained profiles were compared to 3D velocity profiles from GeoDict (FlowDict, Navier-Stokes(-Brinkman) LIR from Math2Market GmbH, Kaiserslautern, Germany) simulations of gas flowing through an equivalent channel domain over a Toray GDL. For comparison, the simulated velocity profile was averaged over the width of the channel domain as well ( Figure A1, hollow lines). The simulated profiles show good agreement with the analytical velocities in the investigated velocity range and only differ slightly, in that the simulated velocity profiles exhibits some penetration into the upper regions of the GDL domain and a slightly shifted peak position in the channel. by first determining ( , ) using the following equation from Shah et al. [44]: where and are the height and width of the channel, respectively, with n = 2.0125 and m = 3.6739 (according to Shah et al. [44] for the given channel geometry). This profile is then averaged over the y-direction to obtain ( ) as depicted in Figure A1 (solid lines).  To capture the inhomogeneous effect of the GDL on the diffusion process, the diffusivity in the GDL domain D GDL was derived from the binary diffusivity D with an in-plane (IP) and TP relative diffusivity coefficient, d IP and d TP , respectively, so that D GDL x = d IP * D and D GDL z = d TP * D. These coefficients were obtained from the literature for bulk GDL properties of Toray GDLs [41,42] and are listed in Table 2.
The formulation of Equation (A3) for each position (x, z) in the model domain with implemented boundary conditions yields a system of linear equations that is solved numerically at each time step for c t x,z . The convergence criteria are set to a relative change in concentration lower than 10 −10 between two iterations. To obtain accurate results, the spatial resolution was refined in both dimensions until a further duplication of resolution resulted in less than 0.5% change of the result, yielding a spatial resolution of ∆x = ∆z = 1 µm.
In order to account for in-plane diffusion effects, the calculation domain in the x direction is required to be larger than L. Therefore, an additional 100 µm was stepwise added on both sides of the evaporation domain, until a further increase of the GDL length led to no significant changes (<0.5%) in the calculated concentration profile. This was achieved at 1.3 mm added length. As the maximum possible gradient on both sides of the evaporation domain is identical and nearly independent of L, this additional channel length is added to both sides of all evaporation domains investigated (see Figure 2).
The total evaporation rate is computed as the stationary flow rate of water exiting the model domain at the channel outlet, as follows: