Lattice Boltzmann simulation of liquid water transport in gas diffusion layers of proton exchange membrane fuel cells: Parametric studies on capillary hysteresis

simulated using a 3D Color-Gradient Lattice Boltzmann model. Simulations were carried out on microstructures of plain and impregnated fiber substrates of a Freudenberg H14. The GDL microstructures were reconstructed from high-resolution X-ray micro-computed tomography ( μ -CT). For the distinction of carbon fibers and polytetrafluoroethylene (PTFE) in the binarized microstructures an in-house algorithm was developed. The additive was specified heterogeneously in the GDL through-plane direction employing a PTFE loading profile as derived based on μ -CT image data. In the in-plane direction the additive was furthermore defined in a realistic fashion near carbon fiber intersections. Prior to parametric studies on capillary behavior a sophisticated modeling approach for semipermeable membranes had to be developed to account for experimental boundary conditions. Capillary hysteresis was then investigated by simulation of intrusion and drainage curves and subsequent comparison to testbench data.


Introduction
In the wake of proceeding climate change, reduction of greenhouse gas emissions becomes an increasingly urgent task. As the transportation sector is responsible for more than 20% of global emissions [1,2], the automotive industry is recently undergoing a disruptive technological transition from conventional combustion engines towards emission-free electrical motors. In the e-Mobility sector clean drivetrain solutions are found in battery-electric (BEV) and fuel cell-electric vehicles (FCEV). The latter is especially compelling for heavy-duty P. Sarkezi-Selsky et al.
typically a carbon-based cloth, paper or felt, which is commonly treated with a hydrophobic agent such as polytetrafluoroethylene (PTFE) to promote its water removal efficiency [7].
Liquid water transport through these GDL microstructures is complex and experimental investigations on the underlying mechanisms are costly and time-consuming. Modeling and simulation of multiphase flows have gained therefore increasing attention in the scientific field within the past decades. The GDL microstructures for such simulations are generally reconstructed from μ-CT imaging data or generated by stochastic methods [8][9][10]. In recent years also modeling of additive distributions has increasingly approached researchers' focus, as the hydrophobicity of PTFE is majorly affecting water transport through GDLs. According to Mathias et al. [7] PTFE is distributed rather homogeneously or heterogeneously along the material thickness, depending on the GDL manufacturing procedure. In previous works such through-plane additive profiles were approximated using either artificial functions [11,12] or experimental data based on SEM-EDX analysis [13,14]. In the in-plane direction PTFE is furthermore predominantly observed in the vicinity of carbon fiber intersections, as it was also visualized in SEM images by Zamel et al. [15]. Most modeling works followed a phenomenological approach, where additive is distributed within these intersecting areas by utilizing morphological image opening [16][17][18]. Other researchers defined distance-based probability functions [12,19] or employed simulated annealing to distribute additives accordingly [20].
In recent years researchers have developed powerful simulation techniques to predict multiphase flows in reconstructed microstructures. Pore scale modeling approaches such as Full Morphology (FM), Lattice Boltzmann method (LBM), Pore Network Modeling (PNM) or Volume-of-Fluid (VOF) method have gained increasing popularity in the scientific field. For more information on different pore scale modeling approaches the reader is referred to [21][22][23][24][25]. Among these methods LBM is especially compelling with respect to numerical stability, suitability for massive parallelization and applicability to arbitrary geometries [26,27]. The latter is particularly advantageous with regard to parametric uncertainty, as simulations can be conducted on lattice representations of real geometries. Modeling of structural effects, e.g., the capillary valve effect [28], thus does not necessarily require additional geometric parameters. To this day a multitude of numerical studies has been conducted, in which liquid water transport through the porous GDL was simulated using the Lattice Boltzmann method. Most of these works focus on investigating liquid water distributions for the intrusion up to breakthrough under variation of parameters such as compression rate [29][30][31], wettability [11,14,29,[32][33][34][35][36][37][38][39] or inlet conditions [40][41][42]. In other studies LBM is utilized to derive effective transport properties based on simulated [37,43,44] or experimentally observed [45,46] saturation profiles. Capillary pressure-saturation characteristics of GDLs, on the other hand, are only very rarely investigated with LBM. Hao et al. [47] employed a free-energy LB model to simulate liquid water transport through stochastic reconstructions of Toray TGP-H-090 carbon papers. The GDL microstructures were modeled with 10 and 30 wt% PTFE, which was randomly distributed on the carbon fiber surfaces. Simulated intrusion and drainage curves showed a fair agreement with experimental data in the intermediate saturation range. Towards low and high saturations the deviations from the measurement increased, which was accredited to the effect of surface pores and deficiencies in the stochastic reconstruction. The simulated curves were furthermore fitted and a modified Leverett function was presented. In addition, a grid resolution of 2.5 μm per voxel was found to be sufficient. This conclusion was, however, drawn solely from a comparison of simulated primary drainage curves for different lattice resolutions.
Yang et al. [36] investigated capillary behavior of Toray TGP-H-090 as well and utilized a pseudopotential LB model. The GDL microstructure was reconstructed using stochastic methods and PTFE contents of 10 and 20 wt% were distributed randomly on the solid surfaces. Depending on the additive load, the GDL reconstruction exhibited a lattice resolution of 1.8 and 2.0 μm per voxel respectively. For the given additive loadings transport of liquid water was then simulated. However, a complete hysteresis was not recovered, as only the intrusion process was investigated. In comparison to experimental data the simulation results showed for both, 10 and 20 wt% PTFE, only a rough qualitative agreement.
Lately Zhu et al. [48] employed a diffuse-interface Lattice Boltzmann model to investigate compression effects on liquid water intrusion. For this, a GDL-type microstructure was randomly generated with stochastic methods and compressed by ratios of 10 and 20% using solid mechanics simulations. The lattice resolution was 2 μm per voxel, an additive was not considered. Comparison of the capillary-pressure saturation characteristics upon intrusion showed an impeding effect of compression on liquid water transport, owing to pore space constriction. The results were, however, not compared to neither experimental nor other simulation data. Capillary hysteresis was furthermore not recovered, as the drainage process was not investigated.
Most recently Zhang et al. [49] studied the intrusion process for a Toray TGP-H-060 carbon paper and a Freudenberg H2315 carbon felt. In the former case the GDL microstructure was generated using stochastic methods whereas in the latter case it was reconstructed from μ-CT image data with a resolution of 1.86 μm per pixel. A PTFE loading of 7 wt% was realized by morphological opening, for the TGP-H-060 furthermore a binder phase was added. Liquid water transport was then simulated for both, in-plane and through-plane direction, using a diffuse-interface LB model. For the Toray carbon paper capillary resistance to water flow was observed to be higher in the throughplane as compared to the in-plane direction. That trend appeared to be reversed for the Freudenberg carbon felt. These findings were explained with different pore size and tortuosity distributions for the two GDLs under investigation. Capillary hysteresis was, again, not part of the study and the simulated behavior was also not compared to any reference.
This research work presents a novel contribution to the field of multiphase flow simulations through porous carbon felt GDLs in following new modeling approaches and elucidating new parametric influences on capillary behavior. For this, a numerical framework is developed to simulate liquid water transport through carbon fiber felts utilizing a 3D Color-Gradient Lattice Boltzmann model and real GDL microstructural data. The carbon felt GDLs are reconstructed from μ-CT image data by manual thresholding and with a lattice resolution of 0.96 μm per voxel, which is at least by a factor of two finer as compared to related previous works [36,[47][48][49]. For the simulations GDL subsets are selected and domain sizes are chosen according to a new approach based on the analysis of porosity profiles for the binarized GDL microstructures. Distinction of PTFE from carbon fibers is also included in this study by development of an algorithm specifying the additive in a realistic fashion near fiber intersections. The differentiation of PTFE and carbon fibers is, in addition, carried out according to a through-plane loading profile, as derived based on μ-CT image data. This combined approach for the distinction of PTFE and carbon fibers in GDL microstructure reconstructions from μ-CT presents a novelty, to our knowledge a comparable route has yet only been taken for stochastic reconstructions [12]. Inspite being usually neglected, are surface effects addressed in this study, as well, by development of two different approaches to model semipermeable membranes according to an experimental setup. Capillary behavior is then thoroughly investigated by simulating both intrusion and drainage of liquid water to recover the complete hysteresis. Relations for the capillary pressure-saturation characteristics are furthermore analyzed upon variation of parameters such as lattice resolution and contact angle. The simulation results are lastly not just validated against theory, but to a meaningful reference of experimental capillary-pressure saturation curves as well.

-CT image acquisition
As the basis for microstructure reconstruction μ-CT imaging is performed on plain and impregnated fibers of a H14 carbon felt GDL by Freudenberg. The scanning is conducted using a SkyScan 1172 desktop μ-CT by Bruker MicroCT (Belgium) and applying the experimental conditions as listed in Table 1. From the raw data single grayscale images with a high resolution of 0.96 μm/pixel are reconstructed using the software NRecon by Bruker MicroCT (Belgium). All reconstruction parameters and correction factors (gray scale, beam hardening correction, ring artifact correction) are kept identical for the different samples as listed in Table 1.

Measurement of capillary pressure-saturation relations
Capillary pressure-saturation (pcS) curves for the carbon fiber felts are measured with an in-house apparatus (pcS testbench) as shown in Fig. 1(a). For this, the GDL sample is sandwiched between a hydrophilic membrane (Merck Durapore Membrane Filter 0.45 μm HV) and a hydrophobic membrane (Whatman Membrane Filter 0.2 μm WTP Type) on the water and gas side respectively. After sealing the micro-fluidic device with tightening screws gas pressure is then controlled using a syringe (Hamilton 1005TLL 5.0 ml) mounted onto a syringe pump (NE-500). A pressure transducer (DMP 311 from BD sensors, ±0.5 bar, ±0.1%) is used with a data logger (Agilent 34970 A Data Acquisition) to record the gas phase pressure during the measurement. Homogeneity in the gas phase pressure across the sample surface is established using a gas distributor on top of the hydrophobic membrane. On the water side pressure is kept constant at 1 atm by providing the liquid via a basin, which is sealed to prevent evaporation. This container is placed on a balance (OHAUS Explorer EX225) to monitor the amount of water intruding into the porous sample. By changing the volume inside the syringe with a constant rate of 50 μl/min, gas pressure is then ramped between +25 kPa and −25 kPa to impose a capillary pressure equivalent to cap = gas − liq across the GDL sample and enforce subsequent intrusion and drainage of liquid water. In preliminary tests it was ensured, that the gas pressure ramp rate was sufficiently low to obtain a quasi-stationary saturation for any given capillary pressure. With this experimental setup capillary pressure-saturation curves are recorded, whereat each measurement consists of 5 subsequent intrusion and drainage cycles to assure reproducibility. For the comparison with simulation results in the further sections of this work only the last cycles are considered.

Binarization of carbon felt GDLs
Microstructures of plain and impregnated fibers are reconstructed from μ-CT image data derived as described in prior subsection. For this, 2D grayscale image stacks from μ-CT ( Fig. 2(a)) are imported into GeoDict [50] and further post-processed. At first, the quality of the raw image data is enhanced by using a sharpening [51] and a median filter [52], each with a radius of one voxel. The filtered data is afterwards binarized using manual thresholding. For the binarization of the fiber substrate a threshold value is chosen to achieve a porosity of 75%. This value is estimated by utilizing manufacturer specifications on an areal weight of 65 g/m 2 and thickness of 150 μm for the H14 carbon fiber felt [53] and an estimated carbon fiber density of 1.75 g/cm 3 based on the measurements of Rashapov et al. [54]. The binarized 3D microstructure of the GDL fibers is visualized in Fig. 2(b). This binarization workflow is the same for the plain and impregnated fibers except for the manual thresholding. Due to lack of a porosity value from literature for the impregnated fibers it is assumed that the additive component (PTFE, conducting soot) is mainly located towards the GDL surface areas, as observed in other works [12,55,56]. The threshold for binarization is hence chosen such as to achieve for the GDL core area a local porosity similar to the unimpregnated fibers. Following this approach resulted in an overall porosity of 71.41% for the impregnated fibers. A comparison of the porosity profiles for the plain and impregnated fiber substrate is shown in Fig. 5

Domain sizes and porosity profiles
In the field of microstructure-resolved simulations it is an everpresent challenge to find a representative elementary volume (REV) being the smallest microstructure subset still recovering the overall macroscopic structural properties [57,58]. In this work a smaller cutout is selected as REV based on an analysis of the porosity profiles for the whole microstructure ( Fig. 2(b)) of the binarized fiber substrate. This approach of REV selection presents, to the best of the authors' knowledge, a novel contribution to the scientific field. Fig. 3 shows the porosity profiles along the two in-plane as well as the throughplane direction of the GDL. The first in-plane direction ( Fig. 3(a)) exhibits a porosity profile with periodic fluctuations. In accordance with the observations on a Freudenberg H2315 by Fishman et al. [59], we also detect locally alternating minima and maxima in the porosity for roughly every 500 μm. This characteristic can be related to the carbon felt manufacturing technique of hydro-entangling, where loose carbon fibers are mechanically intertwined by arrays of water jets [7,12,55,60]. The second in-plane direction ( Fig. 3(b)) is parallel to the production direction and exhibits no comparable periodicity, resulting in a rather uniform porosity profile. The porosity profile through the GDL thickness ( Fig. 3(c)) exhibits a U-type shape typical for carbon felt GDLs [59,61]. Based on this porosity analysis a GDL subset with a domain size of 500 × 200 × 156 voxels is selected as REV. The rationale for this choice is on one hand to cover at least one period of the porosity fluctuations in the in-plane X direction. Due to the uniformity of the porosity in the in-plane Y direction, on the other hand, the system size in the second dimension is kept smaller to reduce computational expense in the simulations. Exploiting the periodicity of the porosity in the in-plane X direction a subset is chosen to encompass a local minimum in the domain center. For the in-plane Y direction the cutout is selected from the center of the full binarized GDL microstructure ( Fig. 2(b)). The cutout positions for the GDL subset are indicated as blue dashed lines in Figs. 3(a) and 3(b). In the through-plane direction the whole GDL thickness is considered. As discussed in Section 5.2 and in the Supplement material, provides this novel approach a GDL subset, which is representative for the twophase transport characteristics of the overall GDL. Fig. 3(c) shows the porosity profiles for the selected subset as blue solid line, the increase in fluctuations arises from the domain size reduction in the two in-plane directions. For the impregnated fibers an appropriate subset ( Fig. 5(a)) was selected as well, following the same approach as described above.

Modeling of surface effects
According to sketch Fig. 1(b) GDL samples in the pcS testbench are sandwiched between a hydrophilic and a hydrophobic semipermeable membrane on the liquid and gas inlet side respectively. To achieve a simulation setup similar to the experimental conditions, these membranes are considered in our modeling as well. In a simple approach these membranes are modeled by adding two monolayers to the GDL microstructure surfaces ( Fig. 4(b)). Assuming, however, that inspite the absence of a nominal compression pressure the membranes face a certain bending force from the air-tight testbench sealing, they are possibly squeezed into the largest pores of the GDL sample surface. In such a scenario a more sophisticated membrane model is required and hence developed in this work. For the binarized GDL subset median and maximum pore diameters are determined by calculating 2D pore size distributions (PSDs) with the Granulometry module in GeoDict for each individual layer along the through-plane direction. At this point it is noted, that this tool calculates geometric pore size distributions (gPSDs) by finding for each individual pore a sphere of maximum diameter, which still can be fitted into the pore volume. When the investigated porous structure, however, is reduced to a thickness of only one single layer, the pore volume essentially becomes a pore area and the sphere of maximum diameter effectively is a circular disk. For anisotropic GDL microstructures it is hence possible that the 2D-gPSDs of individual single slices (2D) are vastly different to the 3D-gPSD of the whole structure (3D). As the GDL porosity increases strongly towards the structure's surface it is furthermore possible that the 2D pore sizes for the individual GDL surface layers increase by orders of magnitude as compared to the GDL core area and might even exceed the thickness of the whole 3D microstructure. Fig. 4(a) shows the 2D pore diameter profiles as well as the porosity profile along the thickness of the GDL subset for the fiber substrate ( Fig. 2(c)). In the core area (enclosed by dashed vertical lines) values fluctuate around 17-20 μm for the median and 40-80 μm for the maximum in the 2D-gPSDs. Towards the GDL surface, however, the 2D diameters increase strongly and concomitantly with the porosity gradient. Based on the through-plane profile of the 2D maximum pore diameter a threshold value of 75 μm is then defined above which pores should be occupied with the respective semipermeable membrane. In applying this rule to the simple membrane monolayer approach (Fig. 4(b)) results in the more sophisticated membrane model in Fig. 4(c). Even though the difference in the two geometries might look small, it plays a crucial role for the hysteresis behavior as will be discussed in Section 5.1. For the 3D-gPSD of the whole fiber subset (Fig. 2(c)) the median is furthermore determined around phys pore ∼ 21 μm which is in range of common literature values [62,63].

Derivation of PTFE volume fraction profiles
Depending on the manufacturing procedure, PTFE is observed to be distributed rather homogeneous or inhomogeneous along the GDL thickness [7,64,65]. In this work such a through-plane loading profile is derived based on μ-CT image data. This is accomplished by comparing the through-plane porosity profiles of the binarizations for the plain and impregnated fibers. The PTFE volume fraction profile is then derived as the through-plane porosity difference for the plain and impregnated fibers (Fig. 5(b)). This analysis yields a volume fraction of 3.68% for the additive. With carbon fiber and PTFE densities of CF = 1.75 g/cm 3 [54,66] and PTFE = 2.15 g/cm 3 [67], this corresponds to an additive loading of 16.36 wt%. This complies to the range of 5-30 wt% commonly denoted in literature [7,9,68]. However, it is noted, that the impregnation solution is assumed to consist solely of PTFE and additional components such as conducting soot are neglected in being out of the scope of this study. For improved readability the terms PTFE and additive are from here on used interchangeably.

Distinction of PTFE and carbon fibers in binarized microstructures
In μ-CT imaging it is usually not possible to differentiate carbon fibers and PTFE due to similar gray-scale values in the images [69,70]. For impregnated fibers a standalone segmentation procedure thus leads to a binarized microstructure of a a void phase and a solid phase mixture of indistinguishable carbon fibers and additive. To alleviate this issue, an in-house algorithm is developed in this work to distinguish carbon fibers and PTFE in binarized microstructures of impregnated fibers. The algorithm mimics realistic additive distributions by specifying the solid material near fiber intersections as PTFE, which is according to experimental surface analysis a preferred location for additive deposition [12,15,55,71,72]. In a first step identification of these fiber intersections is achieved by application of a morphological opening operator [73] on the binarized microstructure of the impregnated fibers ( Fig. 5(a)). For this, a sphere is used as a structuring element with a radius of 6 voxels. This one-step approach, however, would lead to the definition of PTFE as bulky clusters falsely encompassing the whole fiber intersections, including not just additive but also carbon fibers. In a second step PTFE is hence not specified in the bulk of the fiber intersections, but on the solid surfaces in their vicinity. This procedure is carried out under two conditions: • In the through-plane direction additive is assigned with a probability according to the derived volume fraction profile (Fig. 5(b)). • PTFE is only specified on solid surface sites, which are either part of a fiber intersection cluster or neighboring an already present additive surface site. In this way a branch-like propagation of surface PTFE emanating from the fiber intersections is achieved.
As above algorithm only differentiates additive from the overall solid, the porosity of the binarized microstructure is kept constant. Consideration of PTFE as a second solid phase does thus not alter the pore morphology from μ-CT and additional geometric effects such as pore constriction do not occur. Under these conditions PTFE is then differentiated from carbon fibers until the additive covers a desired percentage of the overall solid surface on the impregnated fibers. Experimental investigation of spatial distributions of PTFE, however, is still difficult. SEM-EDX provides in general the technical capability to distinguish PTFE and carbon fibers in GDL samples [7,64,65], but the analysis is oftentimes restricted to higher additive contents due to deteriorating surface effects [69]. In front of these experimental limitations the PTFE surface coverage is therefore estimated in this work with 50%, a value which is amongst others also used by Yu et al. [34]. The final binarized microstructure of the impregnated fiber substrate with separated fibers and PTFE is visualized in Fig. 5(d). As a result of the porosity and additive loading profile ( Fig. 5(b)) the local PTFE surface coverage is heterogeneous (Fig. 5(e)). With regard to additive modeling and simulative studies it is noted, that the PTFE surface coverage provides a more meaningful parameter describing mixed wettability as compared to its loading. This is given by the fact, that hydrophobic additive only affects capillary behavior when it is in contact with the pore space.
according to a single-phase collision term . In this work a single relaxation time operator (SRT) [74] is employed with the relaxation rate = 2∕(6 + 1) setting the time = 1∕ of approaching local equilibrium depending on the fluid's kinematic viscosity . The equilibrium distributions eq are derived from kinetic theory [75] as with a lattice-specific set of discrete velocities Macroscopic local density and velocity entering Eq. (3) are recovered by zeroth and first order distribution moments.
The fluid phase pressure is linked to the density via with a speed of sound of = √ 1∕3 for the D3Q19 lattice. After the collision step the evolved fluid distributions are then propagated in a streaming step to neighboring lattice sites in the directions . (

Multiphase LBM
Multiphase flow is modeled using the Color-Gradient Model (CGM) by Rothman & Keller [76]. In this phenomenological ansatz two phases are described as a 'red (r)' and a 'blue (b)' fluid.
In the CGM the location of the fluid-fluid interface is described by the phase-field parameter , The orientation of the interface is approximated by the color-gradient which is calculated using a fourth order isotropic discretization scheme [77] ∇ = 3 As the kinematic viscosities of the respective fluid phases ∈ { , } are in general unequal, an effective relaxation parameter (Eq. (2)) has to be defined with̄as an averaged kinematic viscosity. In the bulk of the respective phase the effective relaxation parameter is simply given as ef f = . Throughout the multiphase region near the fluid-fluid interface, however, an interpolation scheme has to be employed to ensure a smooth transition from to and vice-versa. In this work ef f is determined by interpolation of the relaxation time ef f = 1∕ ef f along the interface position given by according to the formula of Grunau et al. [78].
Surface tension is generated by utilizing a perturbation operator [79,80] after the single-phase collision step (Eq. (1)) * * ( The perturbation step generates surface tension, but does not provide phase immiscibility. Thus a recoloring operator [82,83] has to be employed in succession , * * * ( , * * * ) The parameter controls the thickness of the fluid-fluid interface and is set to 0.85 in this work. In the streaming step the new distributions of both fluid phases are lastly assigned to the lattice neighbors, similar to Eq. (9) Fluid-solid interactions are primarily governed by the definition of propagation rules near solid boundaries. At ordinary solid surfaces the full-way bounceback rule [75] is employed on both fluid phases. In the case of semipermeable membranes, however, it is only applied to one of the respective phases, allowing only the other fluid phase to permeate the membrane. Wettability is furthermore modeled utilizing the wetting boundary condition by Leclaire et al. [84]. In this formulation the contact angle is imposed locally on the fluid sites neighboring the solid surface by adjusting the angle between the fluid-fluid interface normal = ∕ | | and the wall normal of the solid surface. The adjustment is carried out every time step by minimizing the objective function for the current orientation of the color-gradient , wall normal and the desired contact angle . Eq. (21) is solved iteratively utilizing the secant method ) .
Following Leclaire et al. [84] we also set = 1∕2 and = 2. Eq. (23), however, is different to (1) in [84], as a slightly improved convergence is observed. For the calculation of , fluid densities have to be provided on the solid surface. This is accomplished with the following scheme [85] ( ) = extrapolating the fluid density of the th phase from adjacent fluid sites (index ) to a solid surface site . The weights ( ) correspond to the standard lattice weights for the respective lattice direction to a neighboring fluid surface site.
The wall normal is calculated a priori by applying the isotropic gradient operator (Eq. (13) using = 3 iterations in accordance to [84] and the standard D3Q19 lattice weights for ( 2 + 2 + 2 ) . Above presented LB model is validated against the theoretical test cases of Jurin's [86] and Washburn's [87] law, further details are provided in the Supplement material.

Computational domain and boundary conditions
The simulation setup is generated with the aim to be as close to the experiment as possible, hence applying similar boundary conditions as in the pcS testbench ( Fig. 1(b)). The reconstructed GDL microstructure is sandwiched between a hydrophilic and a hydrophobic semipermeable membrane on the water and gas side, as described in Section 3.3. In the default setup the domain of the porous GDL sample consists of 500 × 200 × 156 voxels with a lattice resolution of = 0.96 μm per voxel. The membranes consist of single layers in the case of the simplified model (Fig. 4(b)). In the more sophisticated membrane model they also partly occupy the pore space in the GDL surface areas with variable penetration depth (Fig. 4(c)).
At the bottom and top boundaries of the computational domain liquid and gas phase pressures are imposed using the pressure boundary condition originally proposed in 2D by Zou and He [88] and implemented by Hecht and Harting [89] for the D3Q19 lattice. In the two in-plane directions the domain boundaries are governed by periodic boundary conditions. Both fluid phases are initialized with a density ratio = 0 liq ∕ 0 gas of unity for 0 liq = 0 gas = 1. Kinematic viscosities are chosen as liq = 1∕6 and gas = 1∕12 to set a dynamic viscosity ratio of = liq liq ∕( gas gas ) = 2. The real ratios in the physical air-water system are with phys ∼ 1000 and phys ∼ 17.5 different, but are not applicable in the LBM simulations due to numerical instability. Application of these real ratios in density and viscosity is, however, not strictly necessary, since water transport within the GDL is known to be a capillary-dominated process and hence inertial and viscous effects are negligible [90,91]. In the simulations this dominance of capillary forces is ensured by setting a high Laplace number of La = pore ∕(2 liq 2 liq ) ≈ 49.5 with a surface tension of = 0.125 and median pore diameter of pore = phys pore ∕ ≈ 22. In preliminary numerical studies it is furthermore observed that viscosity ratios of > 2 do not affect the capillary behavior within this numerical setup significantly.
With above parameters subsequent intrusion and drainage of liquid water is simulated, starting from an initially dry GDL microstructure reconstruction. In accordance with the experiment, liquid water transport is established by imposing a capillary pressure ( cap = gas − liq ) across the GDL via the pressure boundary conditions. Whereas liquid phase pressure is kept constant at liq = 1 atm, gas phase pressure is varied over time via an adaptive pressure ramp.

Adaptive pressure ramp
In the pcS testbench gas pressure is controlled by changing the gas volume inside the syringe. Implementation of such a boundary condition in a numerical setup, however, would be time consuming and with questionable benefit. In addition is an experimental timescale of hours not accessible to LBM simulations, due to the high computational expense. An unreasonable choice of a significantly higher ramp rate, however, is prohibited since quasi-stationarity has to be ensured to prevent dynamic effects in the ramping procedure. Another route is therefore taken, in which the capillary pressure is ramped adaptively as a function of the saturation change in a given time interval. This approach is advantageous with respect to the computational expense, as low ramp rates would be only employed where the saturation responds very sensitive to changes in the capillary pressure, which is usually in the intermediate saturation range. For low and high saturations both, intrusion and drainage curve, generally exhibit flat profiles for which the ramp rate can be increased significantly to reduce the simulation duration. The adaptive ramp for the gas pressure is implemented with the following formula: BC gas ( + ) = BC gas ( ) + BC gas ( ) The gas pressure for the pressure boundary condition at the gas side is hence updated every time step with a ramp rate BC . For intervals of sat = 1000 time steps this rate is adjusted BC gas ( + ) = BC gas ( ) ⋅ sat ( ) ∕ sat (29) by rescaling it with the ratio of the desired saturation change sat and the average saturation change per time step ( ) ∕ sat determined for the preceding interval.

Effect of membrane modeling
In a first computational study the effect of the different membrane models on the capillary behavior of the GDL subset for the plain fibers (Fig. 2(c)) is investigated. Liquid water intrusion and drainage are simulated using both numerical setups as presented in Section 3.3 and compared to experimental data. The simulated hysteresis curves are shown in Fig. 6(a) and major differences become apparent. As indicated by dashed colored lines deviates the first intrusion into the initially dry GDL from the subsequent intrusion for both simulation setups, similar observations are made for the measurement. As this work focuses on the investigation of stationary transport properties, this different capillary behavior upon first intrusion is from here on neglected.
When considering the semipermeable membranes in the GDL surface pores a good agreement between simulation and experiment is achieved. The width of the hysteretic area enclosed by the intrusion and drainage curves also matches well with the experiment. Only minor deviations in the low and high saturation range are found, which could be due to unknown experimental conditions or insufficient lattice resolution in the reconstructed microstructure.
For the simplified membrane model, however, the simulation results deviate significantly from the experiment. As the hysteresis width appears to be similar, both intrusion and drainage curve are shifted towards lower capillary pressures. More crucially, the drainage curve ends at a residual saturation of ∼30%, which is three times higher as compared to the other numerical setup.
To understand this different capillary behavior, it is useful to visualize the liquid water distributions (Figs. 6(b), 6(c)) near the end of the respective drainage curves (red circles in Fig. 6(a)). When the membranes are present in the surface pores, liquid water is drained uniformly from the GDL for increasing capillary pressures and only a few water clusters remain in narrowings inside the porous microstructure. In the simplified membrane model, however, major portions of water are entrapped in the GDL core area, having lost capillary contact to the liquid phase inlet. This observation can be further rationalized with the 2D pore diameter profile in Fig. 4(a) and Laplace's law with the surface tension , carbon fiber contact angle CF and pore diameter pore . During drainage, the non-wetting gas phase finds its way through the GDL along preferential paths with large pores exhibiting the lowest capillary barrier. Starting from the gas inlet side, gas phase first invades large surface pores with little effort. Towards the GDL core area gas phase then progressively has to infiltrate smaller pores, and hence faces increasing capillary resistance. After having passed the GDL core area with the smallest pores representing a capillary bottleneck, gas phase then quickly invades increasingly larger surface pores towards the water inlet side. Once having reached the hydrophilic semipermeable membrane at the water inlet side via the largest available pores, gas phase is pushed back and redistributes in the in-plane direction. Lateral expulsion of liquid water from the surface pores then eventually causes a loss of capillary contact between liquid water at the inlet and inside the pores, resulting in a liquid phase entrapment inside the GDL core area. This scenario is found in Fig. 6(b), the preferential pathway for the gas phase in the midsection of the GDL can be furthermore related to the high porosity in the in-plane porosity profile in Fig. 3(a). For the sophisticated membrane model this behavior is not observed, since the large surface pores are covered by the respective membranes. From this first study it is evident that the residual saturation upon drainage is mainly influenced by gradients in the porosity and pore size profiles and hence surface effects of the GDL. It is furthermore observed, that the sophisticated model of flexible semipermeable membranes being bent into the GDL surface pores provides a more realistic representation of the experimental boundary conditions, leading to a better agreement with the experimental data. For all further studies it is therefore set as default.

Influence of lattice resolution
In microstructure-reconstructed simulations it is of high importance to account for all relevant structural details and thus a high resolution is desirable. Concomitantly the aspected physical domain should be large enough to be representative for the macroscopic sample properties. Delimited by the computational expense of simulation techniques such P. Sarkezi-Selsky et al. Fig. 6. Capillary hysteresis for the plain fibers as simulated with LBM for two different membrane model and compared with experimental data (a). Dashed colored lines show the first intrusion for the respective simulation setups. In both cases the simulations are carried out on a GDL subset for the plain fibers with 200 × 500 × 156 voxels and = 0.96 μm (Fig. 2(c)). The carbon fiber contact angle is set to CF = 65 • . For the data points marked with red circles liquid water distribution is visualized in (b) and (c). The color code for the respective phases is as follows: Semipermeable membrane (yellow), carbon fiber (gray), gas (light-gray), liquid water (blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) as LBM, a compromise often has to be found between richness in detail and system size. In a second study the lattice resolution of the default numerical setup is halved from 0.96 to 1.92 μm per voxel by employing the Rescale module of GeoDict. This leads to a reduced domain size of 250 × 100 × 78 voxels for the plain fiber GDL subset (Fig. 2(c)). In analyzing the influence of the lattice resolution on the simulated capillary behavior it is investigated, if a lower resolution would suffice to account for all relevant details. In that case computational expense could be reduced and a larger physical domain size for the GDL subset would be accessible to LBM simulations. Fig. 7(a) shows simulated intrusion and drainage curves for the GDL subset using a coarse lattice in comparison to the results for the fine resolution and the experimental data. In the simulations on the coarse lattice the surface tension is set to = 0.25 in order to maintain the same Laplace number as on the fine lattice. While the intrusion curves are similar for saturations below 50%, the drainage curves appear comparable in the saturation range above 50%. For the coarse lattice, the curves then deviate increasingly towards higher saturations for the intrusion, and towards lower saturations for the drainage process respectively. These trends render the overall hysteresis width too narrow for the lower resolution. In addition is the residual saturation with ∼5% almost halved, as compared to ∼10% for the fine lattice. Such differences can be related to the lack of structural detail in the staircase approximation of the pore-throat geometry within the coarse lattice. The fact that the deviations in the intrusion and drainage curves do not appear equally along the whole saturation range can be linked to different complexity of local pore morphology, i.e., not all pores and throats require a high-resolution representation. This analysis reveals that a coarse lattice with a resolution of ∼2 μm per voxel is insufficient to recover relevant structural details for an appropriate hysteretic behavior. It is thus decided to use a fine lattice resolution of 0.96 μm per voxel for all subsequent studies. For the fine lattice resolution the plain fiber GDL subset is furthermore found to be representative for the overall two-phase transport characteristics of the GDL. This is confirmed by a similarity in the simulated capillary behavior for different GDL cutouts. Further information is provided in the Supplement material. At this point it is noted, that this second parametric study did not ensure for a resolution of ∼0.96 μm to be already sufficiently high. A grid refinement beyond the μ-CT resolution is technically possible, but would lead to a massive increase in computational expense without providing more structural information of the GDL microstructure.

Contact angle variation
Literature provides a bandwidth of contact angles for water on carbon fibers ranging from mostly 65 • up to 90 • [92][93][94][95]. This might be on the one hand due to a variety in material properties of carbon fibers per se. On the other hand are the derived values also depending on the utilized measurement technique such as, e.g., sessile drop method or Wilhelmy-Plate method. Lastly are measured contact angles usually effective values, since they not just include wetting but structural properties such as surface roughness as well. For an overview on aforementioned aspects the reader is referred to the reviews of Arvay et al. [96] and Ozden et al. [97]. In LBM, however, the contact angle is imposed locally and it should not account for structural properties as they result naturally from the microstructure geometry in the lattice representation. In front of all these parametric insecurities a third numerical study is conducted. Capillary hysteresis is simulated for the plain fiber GDL subset using a carbon fiber contact angle CF of 65 • , 75 • and 90 • respectively. Goal of this investigation is to cover a range of possible hysteretic behavior given by the bandwidth of literature values for CF . Furthermore, it has to be verified, that the preliminary choice CF = 65 • provides the best agreement in simulated and measured capillary hysteresis curve. Fig. 7(b) shows simulated hysteresis curves for the different contact angles. In comparing the simulation data with the experiment it becomes evident, that the best agreement is found for a contact angle of 65 • . Secondly, it is observed that for increasing contact angles the simulation results are striding further away from the measurement, and towards lower capillary pressures. This behavior is as expected, since a lower hydrophilicity leads to hindered intrusion and equivalently promoted drainage of the wetting fluid. For a contact angle of 90 • the carbon fibers become neutrally wetting, which is also reflected in the simulated hysteresis curve being symmetrically centered around 0 kPa, i.e., intrusion and drainage require the same net effort to establish transport through the porous structure. In achieving the best agreement with the experiment for a carbon fiber contact angle of CF = 65 • , this value is set as default for all subsequent studies.

Impact of heterogeneous PTFE distribution
In a final investigation capillary hysteresis is simulated for impregnated fibers to analyze the effect of a hydrophobic additive. For this, Fig. 7. Parametric LBM studies on capillary pressure-saturation hysteresis for the plain fibers and comparison to experimental data. In all simulations the sophisticated membrane model is employed, considering membrane material in the GDL surface pores. The LBM results for the default numerical setup of a GDL subset (Fig. 2(c)) with a domain size of 500 × 200 × 156 voxels, a lattice resolution of 0.96 μm and CF = 65 • , serves as a reference and is depicted as solid line in dark-violet. For each of the investigations only one simulation parameter is changed with respect to the reference setup. In (a) lattice resolution is halved to 1.92 μm, leading to a reduced domain size of 250 × 100 × 78 voxels. In (b) the carbon fiber contact angle CF is varied.
intrusion and drainage simulations are carried out on a GDL subset of segmented fibers and PTFE (Fig. 5(d)). Similar to the default numerical setup for the unimpregnated fibers, the domain size is 500 × 200 × 156 voxels with a lattice resolution of 0.96 μm per voxel. Different wetting behavior is furthermore modeled with contact angles of CF = 65 • and PTFE = 115 • for carbon fibers and PTFE respectively. In Fig. 8(a) capillary hysteresis is compared for plain and impregnated fibers. In both, simulation and experiment, leads PTFE inclusion to a shift of the hysteresis towards lower capillary pressures, owing to the hydrophobic effect of the additive. Another indication for the impeding effect of PTFE on liquid water transport is found in the intrusion curve never reaching a maximum saturation of 100% in neither experiment nor simulation. This can be explained with gas phase entrapment in hydrophobic pores. The hindrance of the wetting fluid's transport through the porous GDL becomes also crucially apparent in visualizing two simulation snapshots marked with red circles in Fig. 8(a). For a capillary pressure of ∼1.37 kPa the plain fibers are already saturated by almost two thirds upon intrusion ( Fig. 8(b)), whereas for a similar capillary pressure of ∼1.25 kPa only ∼13% of the pores are filled in the impregnated fibers (Fig. 8(c)).
With regard to the hysteresis width, the LBM results differ significantly from the measurement for the impregnated fibers. In the simulations the span between intrusion and drainage curve remains similar, despite consideration of PTFE. This is in accordance with observations in previous numerical [47] and experimental studies [98]. In contrast to the simulations appears the hysteresis width to have increased significantly in the measurement. Parametric simulative studies on the effect of additive distributions with different loading profiles and surface coverages of PTFE (not shown) revealed, however, only a minor influence on the hysteresis. The additive distribution can be therefore excluded as a reason for this significant deviation between simulation and experiment. The reliability of the experimental data for the impregnated fibers, however, is to some extend questionable, as the intrusion curve exhibits an intermediate change in curvature between ∼30 and 40% saturation. This characteristic has been also observed in other research works, where it is identified as an artifact from contacting issues between the hydrophilic membrane and the GDL sample. In referring to it as 'shoulder', Gostick et al. [98] accounted this artifact to finite size effects only appearing for thin GDL samples of ≤200 μm thickness. Tranter et al. [99], on the other hand, stated that the 'shoulder' would be mitigated upon increasing compression of the GDL sample. An (unknown) sample compression leading to a narrowing of the pore space could eventually be also responsible for an increase in the hysteresis, as this was also observed in other experimental works [98]. The pcS testbench (Fig. 1) utilized for this work was designed to exert no compressive forces on the GDL sample under investigation. However, this is difficult to warrant, since at the same time the assembly has to be air-tight. Due to the non-linear stress strain relation of the GDL [100] relatively small pressures might already lead to noticeable compression effects.

Conclusion
In this work a new numerical framework was developed to simulate capillary transport of liquid water through porous GDL microstructures using a 3D Color-Gradient Lattice Boltzmann model. Microstructures of plain and impregnated fibers of Freudenberg carbon felt GDLs were reconstructed via binarization from μ-CT images with a high-resolution of 0.96 μm per pixel. For the binarized microstructure of impregnated fibers PTFE was differentiated from carbon fibers with an in-house algorithm specifying the additive near fiber intersections. This was done according to a heterogeneous additive loading profile as derived based on μ-CT image data. Such combined approach for the distinction of PTFE and carbon fibers in GDL microstructure reconstructions from μ-CT presents a novelty. For the simulations GDL cutouts were selected and domain sizes were chosen according to a new approach based on the analysis of porosity profiles for the binarized GDL microstructures. These subsets proved to be representative for the overall two-phase transport characteristics of the GDL. Capillary behavior was then investigated in simulating intrusion and drainage using an adaptive pressure ramp in the numerical setup. The results were compared to experimental data. Moreover, also the width of the hysteresis area spanned by the intrusion and drainage curves was thoroughly analyzed. In several numerical studies different simulation parameters were varied to investigate their impact on the capillary behavior and following key observations were made: • The GDL surface area has a strong influence on the capillary pressure-saturation characteristics owing to a steep pore diameter gradient in the through-plane direction. To recover the same capillary behavior as determined experimentally it is thus important to model surface effects appropriately. In the present work capillary hysteresis was simulated for GDL samples sandwiched P. Sarkezi-Selsky et al. Fig. 8. Comparison of simulated capillary pressure-saturation curves for plain and PTFE-impregnated fibers with experimental data (a). The simulations are carried out on GDL subsets for the plain (Fig. 2(c)) and impregnated fibers (Fig. 5(d)). In both setups equal domain sizes of 500 × 200 × 156 voxels are used, with a lattice resolution of 0.96 μm per voxel. The contact angles for carbon fibers and additive are set to CF = 65 • and PTFE = 115 • respectively. For the data points circled in red, liquid water distribution is visualized for the plain (b) and impregnated fibers (c). The color code for the respective phases is as follows: Carbon fiber (gray), PTFE (green), liquid water (blue). For better visualization of the liquid water distribution both, gas phase and semipermeable membranes, are not displayed. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) between two semipermeable membranes according to the boundary conditions in the pcS testbench. Modeling these membranes as monolayers on the respective sample surface, however, led to massive liquid phase entrapment in the GDL core area upon drainage, overestimating the residual saturation. A more sophisticated approach was necessary, in which the membranes were modeled as being partly bent into the GDL surface area occupying the largest GDL pores. • Lattice resolution is a key parameter for deriving valuable results in microstructure-resolved simulations. In the scope of this work a lattice resolution of ≥ 1.92 μm was found to be too coarse to account for all relevant structural details, resulting in a too narrow hysteresis area between intrusion and drainage curve as compared to the experiment. Since even for a fine lattice with = 0.96 μm (μ-CT resolution), the simulation results still exhibited minor deviations from the measurement, it is to this point unclear if this is already a sufficiently high resolution.
• The contact angle is a fundamental LBM parameter for simulating multiphase flows in porous structures. Experimental assessment of reliable values, however, is difficult owing to the complexity of isolating the wetting properties from structural effects. In a sensitivity study on the carbon fiber contact angle a value of 65 • proved to achieve the best agreement of simulated and measured hysteresis curve. Higher contact angles (75 • , 90 • ) led to increasing deviation from experiment, with a shift towards lower capillary pressures due to increased hydrophobicity of the solid material. The different contact angle values did not affect the hysteresis width.
• A hydrophobic additive in the impregnated fibers led in both, experiment and simulation, to a shift of the capillary hysteresis towards lower capillary pressures indicating a hindrance of water transport. In the numerical setup this impeding effect was realized in modeling mixed wettability of the solid material by assigning different contact angles of CF = 65 • and PTFE = 115 • for the carbon fibers and PTFE respectively. For the experiment the hysteresis width increased significantly upon additive inclusion, whereas this could not be observed in the simulation results.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.