Lattice Boltzmann simulation of liquid water transport in gas diffusion layers of proton exchange membrane fuel cells: Impact of gas diffusion layer and microporous layer degradation on effective transport properties

high durability of each cell component. In this work, transport of liquid water through pristine and degraded gas diffusion layers (GDL) was simulated with a 3D Color-Gradient Lattice Boltzmann model. The GDL microstructure was reconstructed from high-resolution X-ray micro-computed tomography ( μ -CT) of an impregnated Freudenberg H14. The effect of a microporous layer (MPL) was considered by reconstruction of an impregnated and MPL-coated H14. Aged microstructures were generated artificially, assuming loss of polytetrafluoroethylene (PTFE) within the GDL and increase of MPL macroporosity as main degradation mechanisms. Liquid water transport within aged microstructures was simulated by imposing a liquid phase flow rate until breakthrough was reached. Subsequently, the GDL microstructures were analyzed for their breakthrough characteristics by means of saturation and effective gas transport properties. When the MPL was pristine, no distinct GDL degradation effect was observable, this was attributed to the MPL dominating capillary transport. MPL aging, however, led to increased saturations and thus to a deterioration of the effective gas transport. With a partially degraded MPL, aging of the GDL then appeared to affect the breakthrough characteristics.


Introduction
The increasing occurrence of natural disasters highlights the need to take action against climate change more than ever.With the reduction of greenhouse gas emissions as one of the key mitigation strategies, the automotive industry is urged to accomplish a transition from fossil fuelbased combustion engines towards clean powertrain solutions.Among these, battery-electric (BEVs) and fuel cell-electric vehicles (FCEVs) are the most promising candidates for emission-free technologies within i.e., bipolar plates, gas diffusion layers (GDLs), catalyst layers (CLs) and membrane, affecting their functionality [4,5].In order to achieve lifetime prolongation and thus widespread commercialization, underlying degradation mechanisms have to be identified and their effects mitigated.However, this endeavor is complicated by the fact that several components exhibit similar aging behavior leading to superposing degradation effects.
With regard to stable long-term FC operation, proper water management is crucial to prevent flooding of the CL and eventual failure of the fuel cell [6,7].In this context, the GDL plays a pivotal role, as it has to remove water from the CL while transporting reactant gases to it.The GDL is made of a porous carbon-based cloth, felt or paper with a thickness of 100-300 μm [8] and pore sizes around 10-30 μm [9].In order to improve its water removal capacity, the GDL is usually impregnated with a hydrophobic agent such as PTFE [10,11].In recent years, addition of a microporous layer (MPL) was observed to improve the water management significantly in restricting water transport paths and thus preventing flooding [12][13][14][15].The MPL typically consists of a carbon powder mixed with PTFE and is coated onto the GDL [16,17].Common MPL thicknesses are around 20-50 μm [18] and the pore sizes reside in the range from 100 to 500 nm [10,19].At this point it should be noted that other works refer to carbon fiber substrates as gas diffusion backings (GDBs) or macroporous substrates (MPSs), whereas MPL-coated substrates are denoted as GDLs.In the frame of this work, however, GDL and MPL microstructures are treated as separate entities and hence carbon fiber substrates are referred to as GDLs.
In long-term FC operation, both GDL and MPL experience degradation phenomena, which deteriorate the water management capabilities.This leads to increased accumulation of liquid water, reducing the effective gas transport and thus the performance [19,20].Eventually, flooding of the porous microstructures can even cause a complete failure of the fuel cell in operation.According to Park et al. [21], GDL aging occurs mainly due to mechanical or (electro-)chemical degradation.In the former case, the underlying mechanisms can be further subdivided into effects caused by compression, dissolution, erosion or freeze/thaw cycling.(Electro-)chemical degradation, on the other hand, is found to be induced by carbon corrosion or radical attack.Depending on the FC operating conditions, these degradation routes result in structural changes by carbon fiber breakage or loss of carbon material and wettability alterations due to solid surface oxidation or loss of PTFE.For an in-depth description of the aforementioned GDL degradation mechanisms the interested reader is referred to the comprehensive review by Pan et al. [22].Aging of the MPL microstructure is in general observed to occur via the same routes, however, its impact is only rarely differentiated from the degradation effects of the GDL [23,24].In addition, also the membrane and catalyst layer exhibit similar aging behavior and thus superpose the effects of GDL degradation [5,25].In order to isolate the different aging effects, most degradation studies are therefore carried out ex situ, e.g., by immersion in oxidizing baths [26][27][28][29][30][31][32], compression [33,34] or freeze-thaw cycling [35][36][37].Such accelerated stress tests (ASTs) provide helpful insights on possible degradation mechanisms for structural and wettability changes in the GDL/MPL.Being conducted under ex situ conditions, however, they also might have limited transferability to real FC operation.In situ ASTs [38][39][40][41], on the other hand, are costly and time-consuming and come along with the difficulty to differentiate between superposing degradation phenomena.
Given the above difficulties in aging experiments, modeling and simulation of degradation effects have gained popularity in the last decade.Moreover, as the water management within the porous microstructures is governed by complex multiphase transport phenomena, especially pore-scale models (PSMs) have attracted researchers' interests.For an overview on commonly used PSMs we refer to the review of Golparvar et al. [42].
In order to investigate the impact of PTFE loss on liquid water transport in aged GDLs, researchers have increasingly conducted simulative studies within recent years.Seidenberger et al. [43] conducted 3D Monte Carlo simulations of liquid water transport through a stochastic GDL reconstruction.GDL aging was modeled by decreasing the PTFE surface coverage.As a result of reduced surface hydrophobicity, a nonlinear saturation increase was observed for proceeding GDL degradation.In addition, liquid water clusters were found to be discontinuous down to a surface coverage of 65%, beyond which flooding occurred.However, the impact on effective gas transport has not been studied.Pauchet et al. [44] utilized pore network modeling (PNM) to study water transport for GDL degradation due to PTFE loss.Wettability alteration was modeled by increasing the fraction  of hydrophilic pores.Water intrusion into the GDL was then simulated and stopped when the liquid phase reached the outlet, i.e., at breakthrough.Below the percolation threshold   no significant impact of  on the liquid phase invasion patterns was detected.However, for  ≥   breakthrough saturation increased suddenly and nonlinearly.This led to severe deterioration of the effective gas transport.Bosomoiu et al. [45] investigated degradation impacts on the effective gas transport in aged and partially saturated SGL 34BC GDLs.The microstructures were reconstructed from μ-CT, however, carbon fibers and PTFE were not distinguished and the study was restricted to the GDL bulk region.Loss of surface hydrophobicity was assumed as main degradation path and was modeled by a decrease in the surface contact angle.Liquid water intrusion was then simulated with a pore morphology-based model in GeoDict [46].For declining contact angles an increasing flooding tendency, i.e., increasing liquid water accumulation in the GDL, was detected.Consequently, a significant reduction in the effective gas phase diffusivity and permeability was observed.Lately, Yu et al. [47] studied the effect of PTFE content and distribution on liquid water flow as simulated with a 3D Lattice Boltzmann model (LBM).The GDL microstructure was a stochastic reconstruction of a Toray TGP090.Based on their results for varying PTFE surface coverages, the authors concluded that a higher PTFE content leads to shorter breakthrough times and an increasing number of breakthrough locations.However, studies on the effective gas transport were again not conducted in this work.
In the past, PSM has been also employed to investigate the impact of MPL defects on liquid water transport.Wu et al. [48] studied the influence of MPL cracks on liquid water invasion using a cubic PNM.The MPL was considered in a quasi-continuum approach by allowing water permeation only through the cracks of the microporous layer.The GDL saturation was found to be strongly reduced by the MPL and to be dependent on the crack locations.Ma et al. [49] also employed PNM and extracted the MPL pore network from FIB-SEM images.The authors concluded that a fresh MPL reduces liquid water accumulation in the GDL significantly, whereas cracks increased liquid water permeability notably in providing preferential water transport paths.Cetinbas et al. [50] studied liquid water transport in SGL24BA and SGL25BC microstructures with direct numerical simulations (DNS).The geometries were reconstructed from μ-CT for the GDL and nano-CT for the MPL microstructure.The simulation results indicated that a crack-free MPL notably decreases the liquid water content within the GDL, owing to very high breakthrough pressures.A cracked MPL was again found to facilitate liquid water transport.Recently, Hou et al. [51] came to similar conclusions for their 3D Lattice Boltzmann (LB) simulations on liquid water flow in stochastic reconstructions of a GDL and a MPL.While for the crack-free MPL a very high capillary pressure was required, randomly generated cracks were found to significantly promote liquid water transport.As pointed out by Pan et al. [22], above studies provide an insight into the general effect of MPLs on water management.However, these works focus primarily on the difference in liquid water transport characteristics arising from the presence of MPL cracks.Therefore, further investigations on the effects of proceeding MPL aging are yet to be conducted.
In this work, we investigate how wettability alteration and structural changes due to GDL and MPL aging by means of PTFE loss and MPL macroporosity increase affect the liquid water and gas transport through the GDL.At this point it should be noted that a change in wettability can occur either due to solid surface oxidation or loss of the hydrophobic additive.This study focuses on the latter case.However, since the GDL porosity is only slightly altered, loss of PTFE is expected to result in similar aging effects as surface oxidation.Hence, the presented research can be considered as a general investigation on the impact of wettability alteration on multiphase transport within GDLs.While the underlying degradation mechanisms are not subject to this work, studying their effect on GDL transport properties gives important insights on the repercussions of GDL aging on the water management and flooding behavior of the fuel cell.Employing a 3D Color-Gradient Lattice Boltzmann model, we simulate the flow of liquid water through a carbon felt GDL of a PTFE-impregnated Freudenberg H14.The GDL microstructure is reconstructed from μ-CT with a high resolution of 0.96 μm per pixel and distinguishing carbon and PTFE subsequent to binarization.Furthermore, a slice of a microporous layer is reconstructed from μ-CT images of an impregnated and MPL-coated Freudenberg H14 (impregnated H14+MPL), representing the MPL macropores.For the simulations the GDL and MPL microstructures are sandwiched between a liquid and gas phase buffer zone.Mimicking the boundary conditions of real fuel cell operation, liquid water intrusion is then enforced by imposition of a liquid phase flow rate at the GDL inlet side and a constant gas phase pressure at the GDL outlet.For the degradation studies aged GDL microstructures are generated in assuming loss of hydrophobic PTFE as main degradation path for the GDL.Proceeding stages of aging are realized by voxel-wise removal of PTFE from the GDL microstructure with a probability according to the through-plane PTFE surface coverage profile.With a high μ-CT resolution of 0.96 μm per voxel, we are capable of accounting local degradation effects by wettability changes on the microscale.To the best of our knowledge, this modeling of local aging effects by PTFE loss presents a novelty.Moreover, MPL degradation is considered in this work as well.As Freudenberg MPLs are generally known to exhibit rather macropores than cracks [52], we assume here that MPL aging results in changes of the macroporosity.Aged MPL microstructures are therefore reconstructed with increasing macroporosities for proceeding stages of aging.As far as we know, a comparable study on modeling and simulation of MPL degradation effects has not been conducted yet.Degradation studies are then conducted by simulation of liquid water transport through GDL and MPL microstructures at different stages of aging and halted upon breakthrough, i.e., when the liquid phase reaches the GDL outlet side.By variation of one degradation effect at a time, we are able to isolate the impacts of GDL and MPL aging on the water management, a differentiation hard to accomplish in experiments.In light of existing work on modeling and simulation, we also believe that this is the first time a joint study of the effects of both, GDL and MPL degradation is reported.The resulting partially saturated and aged GDL microstructures are then analyzed on their local and overall saturation at breakthrough.Furthermore, repercussions on the effective gas transport are studied by calculation of effective tortuosities and permeabilities in the presence of accumulated water.As a result, we derive dependencies for liquid water content and effective transport parameters as a function of PTFE surface coverage and MPL macroporosity, representing different stages of GDL and MPL aging, respectively.

Experimental
μ-CT imaging is performed for an impregnated H14 and an impregnated and MPL-coated H14 GDL by Freudenberg.Raw data with a high resolution of 0.96 μm per pixel are recorded with a SkyScan 1172 desktop μ-CT by Bruker MicroCT (Belgium) and employing the experimental conditions as listed in Table 1.Single grayscale images are then reconstructed from the raw image data with the software NRecon by Bruker MicroCT (Belgium) and utilizing the reconstruction parameters shown in Table 1.The experimental setup for the μ-CT image acquisition is thus the same as in our previous open-access publication [53].

GDL microstructure reconstruction
The microstructure of a PTFE-impregnated carbon felt GDL is reconstructed based on high-resolution μ-CT data for an impregnated Freudenberg H14 (Section 2).First, the raw structural data is segmented into solid and void phase by binarization using manual thresholding.However, direct distinction of the solid phase into support and additive material is not possible due to similar X-ray attenuation coefficients of carbon fibers and PTFE [54].In succession, an in-house algorithm is therefore employed to distinguish carbon fibers and additive in a realistic fashion.Along the GDL thickness the PTFE content is varied to obtain a heterogeneous loading profile characteristic for certain GDLs [55,56].Furthermore, the additive is primarily specified in the vicinity of fiber intersections, where it is predominantly observed in experimental investigations [57,58].An in-depth description of all the post-processing steps on the reconstruction procedure is provided in our previous open-access publication [53].The binarized GDL microstructure comprises a domain size of 500 × 200 × 156 voxels with a resolution of 0.96 μm per voxel.The porosity is given with 73.22% and the PTFE loading equals 16.36 wt% ( CF = 1.75 g∕cm 3 ,  PTFE = 2.15 g∕cm 3 ).The PTFE surface coverage  equals 50% and provides a more reasonable measure for mixed wettability than the additive loading, as PTFE only affects two-phase transport when in contact with the pore space.Thus, GDL degradation is hereinafter mainly discussed by means of PTFE surface coverage reduction.
This study aims at mimicking in operando conditions.As fuel cells typically operate under certain compression pressures, the binarized GDL geometry is compacted artificially beforehand further investigations.The authors are aware of the fact that μ-CT imaging and reconstruction of compressed GDL samples presents an alternative approach.However, as a compression unit introduces a source of interference for X-ray interaction, sample compression generally causes a deterioration of the μ-CT image quality.For this work an artificial compression with a numerical procedure is therefore favored.With the GeoDict module Compress structures are compacted geometrically by reduction of their thickness according to a preset compression ratio  = 1 −  1  0 of reduced  1 to initial thickness  0 .Here, we assume  ∼ 0.22 for a compression load of 2 MPa and manufacturer specifications for a plain H14 [59] and an impregnated H14+MPL [60].This is also in line with a comparable work of Chuang et al. [61] for a Freudenberg H23C8.With an uncompressed thickness of  0 = 150 μm (156 voxels) this results in a reduced GDL thickness of  1 = 117 μm (122 voxels).At this point we acknowledge that a physics-based compression routine such as Elasto-Dict provides in general more accurate results.However, a preliminary study on a smaller GDL sub-domain showed for lower compression rates similar structural and multiphase transport characteristics for the two aforementioned compression approaches using Compress and ElastoDict.
In favor of computational efficiency the former is selected as preferred.After compression the GDL microstructure exhibits a reduced porosity of 65.8% and a PTFE surface coverage of  = 47.0%.The resulting geometry is illustrated in Fig. 2(c), for the remainder of this work we will refer to it as the pristine GDL.

MPL microstructure reconstruction
The macropores/cracks of a MPL play a crucial role for the liquid water invasion into the GDL backing.Therefore, a single layer describing the MPL macroporosity is reconstructed based on high-resolution μ-CT imaging data of an impregnated H14+MPL (Section 2), which acts as an inlet layer for the water intrusion into the GDL.Analogous to the GDL microstructure reconstruction, the raw image data exhibits a resolution of 0.96 μm per pixel and is filtered prior to binarization (Section 3.1 and [53]).For a slice halfway along the MPL thickness the microstructure is then reconstructed by manual thresholding, assuming a macroporosity of  macro MPL = 20 % for the pristine MPL (Fig. 2(f)).With this choice for  macro MPL we orient ourselves to the work of Muirhead et al. [62] for a SGL 25 BC.As will be shown in a subsequent part of this work (Section 6.3), however, the impact of  macro MPL on the breakthrough characteristics appears to become insensitive below a certain threshold.Parametric uncertainty from estimation of the pristine MPL's macroporosity is therefore not regarded as an issue of relevance.Furthermore, it is noted that the reconstruction is restricted to the macropores of the MPL structure, as nanopores are beyond the μ-CT resolution.A visualization of these different pore sizes is provided by scanning electron microscopy images of the impregnated H14+MPL in the Supporting material.Here, negligence of the microporosity is justified by the fact that liquid water is mainly transported through the largest MPL openings with the lowest capillary resistance [13,[63][64][65].Partial reconstruction limited to the MPL macropores is thus expected to be sufficient for the liquid water transport characteristics.Furthermore, the binarized solid material is treated as one component and no distinction is made between carbon powder and PTFE.

Multiphase Lattice Boltzmann -The Color-Gradient Model
In this work, liquid water transport through porous GDL microstructure reconstructions is simulated with the Lattice Boltzmann method (LBM) and the D3Q19 velocity set.Multiphase interactions are modeled based on the Color-Gradient Model (CGM) by Rothman & Keller [66].An in-depth description of the implemented model and all governing equations can be found in a previous open-access publication of Sarkezi-Selsky et al. [53] in which the LB model is also validated against theory and experimental data.

Computational domain and boundary conditions
The numerical setup aims at mimicking the real boundary conditions during fuel cell operation, the GDL microstructures are therefore considered to be in the compressed state (500 × 200 × 122 voxels).In order to investigate possible effects of MPL degradation on liquid water transport characteristics, a MPL is considered as well.However, due to the coating procedure the MPL is found to penetrate into the large pores of the GDL surface regions with high porosity, resulting in a transition zone between GDL and MPL microstructure [52,[67][68][69].In a straightforward approach one would resolve this region by a full reconstruction of the impregnated H14+MPL.This task is, however, beyond the scope of this work due to three reasons: 1. Delimited by the μ-CT resolution, MPL reconstruction is restricted to its macropores.2. The porosity of the impregnated H14+MPL is denoted by the manufacturer [60].However, usage of this information for a full binarization of the impregnated H14+MPL would lead to unknown reconstruction errors, since the MPL microporosity cannot be resolved.3. Due to different grayscale values, joint reconstruction of GDL and MPL microstructures requires trinarization (distinction of three phases), which is a disproportionately more challenging task in comparison to binarization.In front of the aforementioned circumstances another approach is taken in this work.
Here, the microstructures of GDL and MPL are reconstructed separately as described in Sections 3.1 and 3.2.Furthermore, the penetration depth of the MPL is estimated by comparison of the porosity profiles for the impregnated H14 and the impregnated H14+MPL.Based on these profiles the transition zone GDL/MPL is then identified for the uncompressed microstructure reconstructions as obtained by binarization of μ-CT images.Since both samples contain the same GDL microstructure, the uncompressed impregnated H14+MPL is binarized using the same grayscale threshold value as for the uncompressed impregnated H14.With respect to the MPL microstructure and its porosity, we note that above binarization approach only allows for a qualitative analysis of the porosity profiles.This is, however, sufficient for the identification of the GDL/MPL transition zone.Fig. 1(a) shows the derived porosity profiles for the binarized and uncompressed microstructures of the impregnated H14 and the impregnated H14+MPL.As a result of the same underlying GDL microstructure, both profiles are in major parts similar.Smaller deviations in the GDL bulk region are observable, but can be related to fluctuations in the GDL manufacturing process.In the GDL surface region towards the MPL, however, the porosity of the impregnated H14+MPL is observed to be increasingly different.Based on the change in curvature it is then assumed that the transition zone GDL/MPL resides in the GDL surface region as shown in Fig. 1(a).Our approach is therefore comparable to the one by Chevalier et al. [68].We transfer above findings to the compressed GDL microstructure reconstruction (cf.Section 3.1, Fig. 2(c)) and identify the respective surface region in the porosity profile (see Fig. 2(b)) as transition zone GDL/MPL.
For the LB simulations the identified transition zone is cut off, reducing the domain size of the compressed GDL microstructure (Fig. 2(c)) to 500 × 200 × 101 voxels.Subsequently, a single layer of reconstructed MPL macropores (Fig. 2(f)) is added to the bottom of the cut and compressed GDL.The combined GDL/MPL domain is then set to be initially dry, i.e., filled with gas phase and modeled as being embedded between the CL and gas channel (GC) at the bottom and top of the GDL microstructure, respectively.
Buffer zones are defined with heights of 10 voxels for the liquid phase on the CL side and 30 voxels for the gas phase on the GC side.Thus, the overall domain size of the numerical setup comprises 500 × 200 × 142 voxels.Intrusion of liquid product water from the CL into the GDL is then modeled via imposition of a liquid phase flow rate at the bottom domain boundary.At the top domain boundary on the GC side, gas phase is provided with a constant pressure of  gas = 1 atm.In the simulations this is realized with a velocity boundary condition for the liquid phase at the lowermost boundary layer and a pressure boundary condition for the gas phase at the topmost domain layer.The formulations for these external boundary conditions are adopted from Hecht & Harting [70], based on the original works of Zou & He [71].In the two in-plane directions the domain boundaries are governed by periodic boundary conditions.Fig. 1(b) illustrates the final computational domain of the GDL sample embedded between the CL and GC buffer zones for the respective fluid phases.
Fluid phases are initialized with  0 liq =  0 gas = 1, resulting in a density ratio  =  0 liq ∕ 0 gas of unity.Kinematic viscosities are set to  liq = 1∕6 and  gas = 1∕12, leading to a ratio in the dynamic viscosities of  =  liq  liq ∕( gas  gas ) = 2.The real ratios of the physical airwater system are given with  phys ∼ 1000 and  phys ∼ 17.5, but are not accessible in LB simulations due to limited numerical stability.However, since liquid water transport through the GDL is a capillarydominated process, this is not strictly necessary as inertial and viscous effects are negligible [72,73].In order to assure dominance of capillary  Starting from initially dry GDL samples, liquid water intrusion is then simulated with above parameters and stopped upon breakthrough, i.e., when the liquid phase reaches the top domain end of the GDL.

Calculation of effective transport properties for partially wetted GDLs
The LBM simulations are stopped at breakthrough and the partially saturated GDL domains are extracted as simulation snapshots for further evaluation.Subsequently, the liquid water distribution within the porous microstructures is analyzed.First, the voxel-wise saturation is calculated as for each pore space voxel at position (, , )  .Saturation profiles along the GDL thickness are then determined as the average over all  p  and  p  pore space voxels of the respective in-plane direction.The total GDL saturation is eventually derived as the global average over all  p  ,  p  and  p  pore space voxels along the three spatial dimensions.Liquid water within the partially saturated GDL/MPL samples occupies available pore space.This has repercussions on the gas transport and in real fuel cell operation eventually on the performance as well.In order to investigate this degradation impact, we determine tortuosities and permeabilities for the aged and partially saturated GDL/MPL samples using GeoDict simulations.Prior to these studies, the two-phase regions within the partially wetted microstructures are binarized into single-phase regions of the respective fluids.This distinction is made for every pore space voxel by thresholding via the phase-field parameter   (, , ) =  liq (, , ) −  gas (, , )  liq (, , ) +  gas (, , ) ( denoting the position of the fluid-fluid interface.Moreover, a link to the saturation is established with  = 2  − 1 by comparison of Eqs. ( 4) and (1).With the interface mid at  = 0, the two-phase region is differentiated into liquid phase for  ≥ 0 and gas phase for  < 0. In the GeoDict simulations the liquid phase is then considered as an immobile phase blocking progression of the gas phase.In addition, the solid material of the MPL is assumed to allow gas transport on the unresolved nanoscale.In order to investigate how the diffusive gas transport is affected by the pore blockage due to presence of liquid water, the tortuosity  within the partially wetted GDLs is then determined. is a measure for the diffusion resistance in a porous medium and is defined as with  as the porosity of the porous medium. * denotes the relative diffusivity given as the ratio of the effective diffusivity  ef f and the self-diffusion coefficient  0 .The characteristic length of the compressed GDL microstructure (Fig. 2(c)) is with a median pore diameter of around  phys pore ∼ 18 μm much larger than the molecular mean free path  air = 68 nm of air [76].For a Knudsen number of Kn =  phys pore ∕ air ∼ 3.5e −3 the diffusive transport in the GDL is therefore mainly governed by molecular diffusion, while Knudsen diffusion can be neglected.In this work, we determine the effective diffusivity  ef f with the GeoDict tool DiffuDict by solving Fick's law for a given diffusive flux  and concentration gradient ∇ across the porous GDL microstructure.The tortuosity is then eventually obtained via Eqs.( 5) and ( 6).
In addition to the diffusive transport, effective gas permeabilities  gas for the partially saturated GDLs are investigated with GeoDict simulations as well.Under the assumption of creeping (Stokes) flow  gas is determined with the GeoDict module FlowDict by solving Darcy's law  gas = −  gas  gas ∇ gas (8) for a given gas velocity  gas , dynamic viscosity  gas and pressure gradient ∇ gas .
In fuel cell operation the gas transport along the GDL thickness (from GC to CL) is in general more crucial as compared to the in-plane directions, we therefore restrict our investigations on the effective gas transport to the through-plane direction.

GDL microstructure degradation due to PTFE loss
As already pointed out in Section 1, there exist several GDL degradation routes, resulting in structural and/or wettability alteration.However, the simulative studies of this work focus solely on the effects of PTFE loss.According to Pan et al. [22], the underlying degradation routes can be divided into material dissolution during normal operation, freeze-thaw cycling at start-up/shutdown or radical attack in idling periods.While above mechanisms for additive degradation are already identified, literature on detailed quantitative and spatial analysis of PTFE loss is to our knowledge still not available.Hence, an algorithm is developed in this work to artificially realize additive degradation in reconstructed microstructures by voxel-wise removal of PTFE from the solid surface.Basic assumption of this artificial aging routine is that additive loss exhibits a stochastic nature and proceeds in dependence of the PTFE abundance on the solid surface.Along the GDL thickness, additive voxels are therefore selected for depletion with a probability based on the PTFE surface coverage profile (Fig. 2(a)).This modeling approach results in two structural aging effects: • As hydrophobic additive is removed, the solid surface becomes more hydrophilic.• Loss of PTFE leads to a reduction of the solid volume fraction, i.e., porosity increase.
Starting from the microstructure reconstruction of the pristine GDL (Fig. 2(c)), PTFE is then incrementally removed corresponding to proceeding stages of aging.In this way, 11 aged GDL microstructures are generated for the subsequent simulative studies, with additive surface coverages ranging from 47% (pristine) to 0% (completely degraded).Aging effects of above described scheme are illustrated for three selected degradation stages of the compressed GDL microstructure in Figs.2(c)-2(e).The repercussions on the porosity and the PTFE surface coverage are elucidated in the through-plane profiles in Fig. 2(a) and Fig. 2(b).Prior to LB simulations, the partially aged GDL microstructures are then cut by removal of the GDL/MPL transition zone as described in Section 4.2.This slightly affects the overall porosity and PTFE surface coverage.A summary of the structural properties for the GDL samples at varying stages of aging is provided in Table 2.

MPL microstructure degradation due to increase in macroporosity
In this work, MPL aging is assumed to result in an increase in the macroporosity.While the approach for microstructure reconstruction of the GDL (Section 3.1) and MPL (Section 3.2) is in principle the same, modeling of structural degradation effects is conceptually different for the respective two cases.The GDL microstructure is reconstructed from μ-CT exactly once (Section 3.1), resulting in a binarized microstructure representative for the pristine GDL (Fig. 2(c)).Based on this reconstruction, different stages of artificial aging are then realized by employing the PTFE degradation routine as described in Section 5.1.

Table 2
Structural properties for GDL microstructures at proceeding stages of aging due to PTFE loss.The terms 'compr.' and 'compr.+ cut' refer to the GDL microstructures after compression (Section 3.1, Fig. 2) and after compression and removal of the GDL/MPL transition zone (Section 5.1, Fig. 1 In order to model the structural changes of MPL aging, on the other hand, the MPL microstructure is separately reconstructed from μ-CT for each stage of degradation exhibiting a different macroporosity.Starting with  macro MPL = 20%, a binarized microstructure representing the pristine MPL (Fig. 2(f)) is reconstructed from μ-CT as described in Section 3.2.This reconstruction from μ-CT is then repeated for varying stages of aging, however, individually adjusting the manual thresholding to obtain a higher  macro MPL in each of the binarized MPL microstructures (Fig. 2(g)).In spite of representing aged MPL microstructures, all of these reconstructions are therefore based on the μ-CT imaging data for a pristine impregnated H14+MPL (Section 2).Analogous to the degradation of GDL microstructures (Section 5.1), MPL aging is thus purely virtual in this work.The macroporosities for the reconstructions at different stages of MPL degradation are summarized in Table 3.

Simulative studies on GDL and MPL aging
The impact of GDL and MPL degradation on the liquid water breakthrough characteristics is investigated in three series of LB simulations.In order to isolate the aging effects of GDL and MPL within each of these studies, only one of the two respective microstructures is considered for proceeding stages of aging, while the other is held at a certain degradation state.At first, GDL degradation is investigated for varying degrees of PTFE loss (Table 2) and assuming a pristine MPL ( macro MPL = 20%, Fig. 2(f)).In a second study, the focus resides on degradation effects of the MPL and liquid water transport is simulated for a pristine GDL ( = 47.0%,Fig. 2(c)) microstructure and increasing MPL macroporosities (Table 3).Lastly, GDL aging is studied once again, but assuming a partially degraded MPL with a macroporosity of  macro MPL = 60% (Fig. 2(g)).MPL as primary degradation effect in MPL microstructures: Pristine and degraded MPL microstructures are obtained with binarization by manual thresholding from high-resolution μ-CT (0.96 μm per pixel) of an impregnated H14+MPL.The degradation procedure is illustrated in (f) and (g).

Influence of liquid water flow rate on saturation profiles
During fuel cell operation, liquid water is transported through the porous GDL at very low flow rates and gravitational and inertial effects are usually negligible.Under these circumstances, liquid phase flow patterns can then be characterized by the Capillary number and the dynamic viscosity ratio Ca is a dimensionless number denoting the ratio of viscous to capillary forces and is defined with the velocity  liq and dynamic viscosity  liq of the liquid phase and the surface tension .In Eq. ( 10), μ gas denotes the dynamic viscosity of the gas phase.According to Lenormand et al. [77], transport regimes in porous media can be described as either viscous fingering, stable displacement or capillary fingering, depending on Ca and .With capillary numbers around Ca ∼ 10e −5 − 10e −8 and a viscosity ratio of  = 17.5, liquid water transport within the GDLs of operating fuel cells can be thus classified as capillarydominated flow [78].For GDL thicknesses commonly ranging from 100 to 300 μm [8], breakthrough of the liquid phase therefore can take up to several minutes [79,80].Due to the small time step sizes in microstructure-resolved simulations (such as LBM), this would result in large numbers of time steps beyond computational accessibility.
A first investigation therefore aims at finding the highest flow rate, below which the resulting saturation profile becomes stationary, i.e., independent of the liquid phase velocity.In that case, in operando liquid water transport characteristics become accessible to LB simulations, despite significantly higher flow rates than in real fuel cell operation.As liquid water transport within the GDL microstructure is a capillarydominated process, dynamic effects should be negligible and capillary behavior would therefore be invariant to minor variations in the flow velocity.In a sensitivity study, different flow rates are imposed on the pristine GDL and MPL via the velocity boundary condition on the GDL inlet side (cf.Section 4.2).Liquid water invasion is then simulated and halted upon breakthrough.Subsequently, liquid water distributions are characterized by means of time dependence for breakthrough occurrence and through-plane saturation profiles.
In Fig. 3 liquid water invasion behavior is illustrated for 5 different liquid phase flow rates  liq as given in lattice units (corresponding physical counterparts are provided in Table 4).In order to investigate the similarity of liquid phase invasion behaviors, the time dependence of the breakthrough event is analyzed first.Under the assumption of a capillary-dominated process, liquid water invasion behavior should be insensitive to flow rate variation.In that case, the time to breakthrough would follow a reciprocal dependence on the imposed liquid phase flow velocity, i.e.,  BT ∼ 1∕ liq .Fig. 3(a) shows that this relation holds well for  liq ≤ 5e −5 , suggesting similar two-phase transport characteristics.For higher flow rates, breakthrough occurrence deviates increasingly towards longer time spans, indicating different liquid phase invasion patterns.This first conjecture is confirmed by a comparison of the breakthrough saturation profiles in Fig. 3(b).For the four lowest flow rates, the liquid water distribution along the GDL thickness is found to be similar with only insignificant fluctuations.At  liq = 1e −4 , on the other hand, the saturation profile is observed to be already slightly elevated.For even higher flow rates, significantly larger amounts of liquid water have accumulated within the GDL upon breakthrough.This change of intrusion pattern is also illustrated for three different liquid phase flow rates in Figs.
. The observed tendency of the saturation concomitantly increasing with higher flow velocities can be related to the build-up of dynamic pressure, facilitating the pore-filling process.In this case, the invasion behavior is deteriorated by dynamic effects and hence not representative for the in operando conditions.The implied existence of different flow regimes can be also verified with the Capillary number Ca (cf.Table 4).According to Lenormand et al. ([77], Section 4.2), flow regimes beyond  > 0 and Ca ≤ 1e −4 are classified as capillary fingering.For larger capillary numbers, the influence of inertial effects increases progressively.The flow regime then undergoes a transition, until for Ca > 0 it exhibits the characteristics of stable displacement.
Based on above study, a flow rate of  liq = 5e −5 is found to be sufficiently low to ensure a capillary-dominated regime and is therefore set as default for all further LB simulations in this work.In addition, it is noted that the current density (see Table 4) corresponding to the lowest simulated flow rate  liq = 5e −6 is already near realistic values.Since also higher liquid phase velocities of  liq ≤ 5e −5 lead to similar flow patterns, the findings of the subsequent studies on breakthrough characteristics can be related to real fuel cell operation.

Impact of GDL aging on breakthrough characteristics
In a second study, the effect of PTFE degradation on the liquid water transport characteristics within the GDL is investigated.LB simulations on liquid water intrusion are conducted for the GDL microstructure reconstruction at 11 stages of aging as summarized in Table 2.In order to isolate the GDL aging effects, the MPL is considered to stay pristine throughout this investigation.The simulations are stopped upon breakthrough and the respective partially wetted microstructure samples are characterized by analysis of the liquid water distribution.Fig. 4 shows that the loss of PTFE does not have any significant effect on the breakthrough characteristics.Even though the surface wettability and the porosity of the GDL change for proceeding stages of degradation, this apparently does not result in different liquid water distributions and the overall saturation at breakthrough.This observation could be attributed to the pristine MPL having a dominant influence on the breakthrough characteristics by strongly restricting the water intrusion paths into the GDL.In representing a bottleneck for capillary transport, the MPL would therefore superpose the degradation effects stemming from the aged GDL microstructure.Such superior impact of the MPL properties on the transport characteristics has been also reported by Jinuntuya et al. [81].However, the dominance of the MPL can diminish significantly once it is degraded, as will be illustrated hereinafter.Another explanation for an absent degradation effect might be a characteristic of Freudenberg GDLs to transport liquid water mainly along the GDL thickness and hardly in the in-plane direction [82].With only minor inplane redistribution, the invasion paths for liquid water breakthrough are short and weakly tortuous.The effect of the hydrophobic additive (PTFE) on the GDL fiber surfaces is hence limited and so might be the degradation impact of PTFE loss.Due to the absence of a significant aging effect, a subsequent investigation on the effective gas transport within the partially wetted microstructures is omitted for this data set.

Impact of MPL aging on breakthrough characteristics
In another degradation study of this work, the influence of MPL aging on the liquid water transport is investigated.For this, LB simulations on liquid water transport are conducted for the pristine GDL ( = 42.8%,Fig. 2(c)) and the MPL at 9 different stages (cf.Table 3) of aging from pristine ( macro MPL = 20%, Fig. 2(f)) to fully degraded ( macro MPL = 100%).As apparent from Fig. 5, MPL aging has a major effect on the liquid water breakthrough characteristics.For increasing macroporosities  macro MPL , liquid water transport is found to be progressively promoted, resulting in larger amounts of accumulated water within the GDL (Fig. 5(a)).This trend can be related to the fact that liquid water intrusion into the GDL is strongly restricted by the macropores of the MPL.With higher  macro MPL , size and number of these macropores increase and thus facilitate liquid water transport through larger and also new invasion paths.This circumstance is also illustrated with the liquid water distributions for three stages of proceeding MPL aging in Figs.5(c)-5(e).In addition, a transitional region is located around 40% ≤  macro MPL ≤ 50%, with a steep gradient in the local and overall saturations at breakthrough.Below and above this range of macroporosities, however, the amount of liquid water accumulating in the GDL appears to be less sensitive to changes in  macro MPL .This observation can be related to the circumstance that liquid water transport into the GDL is not just delimited by the MPL macropores, but rather by the effective cross-sectional area formed by the MPL and the adjacent GDL surface layers.Furthermore, Fig. 5(b) shows how the overall saturation  total  increases for proceeding MPL degradation.For the dependence on  macro MPL , an inverse tangent relation is determined via nonlinear curve fitting.
With the liquid water accumulating in the GDL and thus blocking pore space, we observe that the effective gas transport at breakthrough is progressively deteriorated upon proceeding MPL aging (Fig. 6).Starting off with the pristine MPL, the tortuosity  wet gas (Fig. 6(a)) of the partially wetted GDL is determined to be still almost the same as for the dry state with  dry gas .In contrast to this, the effective permeability  wet gas (Fig. 6(b)) is observed to be already significantly reduced as compared to  dry gas .An almost unaffected tortuosity combined with a diminished permeability might again indicate the characteristics of Freudenberg GDLs to transport liquid water along the thickness of the GDL and barely in the in-plane direction [82].For proceeding stages of MPL aging, the tortuosity  wet gas then increases nonlinearly with  macro MPL , whereas the permeability  wet gas declines linearly with higher MPL macroporosities.In the case of the fully degraded MPL ( macro MPL = 100%), however, gas transport appears to be completely inhibited by the presence of liquid water ( wet gas = ∞,  wet gas = 0).

Impact of GDL aging on breakthrough characteristics for a partially degraded MPL
In a final investigation GDL degradation is studied again, but with a partially aged MPL ( macro MPL = 60%, Fig. 2(g)).In contrast to the absence of any degradation impact with the pristine MPL (Fig. 4), we now observe for the partially aged MPL distinct GDL aging effects on the liquid water transport characteristics at breakthrough (Fig. 7).Along with progressive PTFE loss, larger amounts of liquid water are now accumulating, especially in the inlet half of the GDL (Fig. 7(a)).As evident from Fig. 7(b), this results in an overall saturation increase with reciprocal dependence on the surface coverage  of the hydrophobic additive.This trend is similar to the findings of Seidenberger et al. [43], albeit for lower PTFE surface coverages.Furthermore, a comparable nonlinear saturation change is also observed by Wang et al. [83], however, in their work it is in dependence of the additive loading.The observed degradation effect of PTFE loss on the breakthrough characteristics is also illustrated in the liquid phase distributions for three stages of proceeding aging in Figs.7(c)-7(e).Above observations lead to the conclusion that partial degradation of the MPL apparently diminishes its influence on the liquid water transport characteristics.Since the aged MPL is not anymore fully dominating capillary transport, GDL degradation affects the breakthrough characteristics.Eventually, the effective gas transport is deteriorated upon proceeding stages of GDL aging as well (Fig. 8).This is indicated by rising tortuosities  wet gas and decreasing permeabilities  wet gas , both with nonlinear dependencies  2(c)).Liquid water distribution within the GDL is illustrated for three stages of proceeding MPL degradation in (c)-(e).The color code for the respective phases is as follows: Carbon fiber (gray), PTFE (green), 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.)Fig. 6.Impact of MPL degradation on effective gas transport characteristics at breakthrough.Higher liquid phase saturations (Fig. 5) lead to increasing pore blockage for proceeding stages of aging.This deterioration is indicated in a progressively (a) higher tortuosity  wet gas and (b) lower permeability  wet gas of the gas phase in the partially wetted pristine GDL.
on the PTFE surface coverage .Furthermore, we derive reverse degradation trends for  dry gas and  dry gas in the dry GDL.These stem from the increase in GDL porosity  GDL upon PTFE loss as illustrated in Fig. 2(b).

Conclusion
In this work, the impacts of GDL and MPL aging on the breakthrough characteristics of the GDL were investigated.For this, liquid water transport through artificially degraded GDL and MPL microstructures was simulated with a 3D Color-Gradient Lattice Boltzmann model.GDL and MPL microstructures were reconstructed from high-resolution μ-CT data for a PTFE-impregnated Freudenberg H14 and an impregnated H14+MPL, respectively.Degraded GDL microstructures were generated, assuming PTFE loss as a primary degradation mechanism.Aging of a MPL layer was modeled by an increase in macroporosity.
Combining different stages of GDL and MPL aging, liquid water intrusion was then simulated and stopped upon breakthrough.In order to investigate the degradation effects on multiphase transport, the partially wetted GDL microstructures were then characterized by means of local and overall saturation.In addition, repercussions on the effective gas transport were analyzed by calculation of gas phase tortuosities and permeabilities.In the frame of the underlying investigations, the following key observations were made: • A sensitivity study revealed the presence of a capillary-dominated regime for the liquid phase invasion up to velocities of  liq = 1.14e −4 m/s.Below this threshold the flow characteristics become independent of the velocity.This criterion enables a significant reduction of the computational costs for the LBM, as it allows using much higher flow rates than present during fuel cell operation.7) lead to increasing pore blockage for proceeding stages of GDL aging.This deterioration is found in a progressively higher tortuosity (a) and lower permeability (b) for the gas phase in the partially wetted GDL.The secondary degradation effect of porosity increase is embodied in decreasing tortuosities and increasing permeabilities for the dry GDL.
• When combined with a pristine MPL, degradation of the GDL microstructure did not result in any distinct aging effect on the liquid water breakthrough characteristics.These observations were explained with the pristine MPL dominating capillary transport and thus superposing any GDL degradation effects.• MPL aging was observed to promote liquid water transport within the pristine GDL significantly as indicated by higher breakthrough saturations.A steep gradient in the accumulated liquid water was furthermore located at macroporosities around 40% ≤  macro MPL ≤ 50%.With the liquid water accumulating in the GDL pores, the effective gas transport is deteriorated by MPL degradation.This circumstance is reflected in a nonlinearly increasing tortuosity and a linearly decreasing permeability.
• After a certain degree of aging, the MPL loses its full dominance on capillary transport and GDL degradation effects become noticeable as well.For proceeding PTFE loss, breakthrough saturations are found to increase significantly due to the reduction of surface hydrophobicity.As a result of progressive pore blockage, the effective gas transport is hampered upon GDL degradation.This is indicated by a tortuosity increase and a permeability decrease, both following a nonlinear dependence on the PTFE surface coverage.
Above simulative studies provide new insights into the repercussions of both, GDL and MPL aging on the liquid water transport within the GDL.This work thus contributes to the understanding of degradation effects and their impact on PEMFC water management.Combining GDL and MPL microstructures at different stages of degradation, respective aging effects on the GDL breakthrough characteristics were studied.We have hence isolated the impacts of GDL and MPL degradation, which is usually hard to accomplish in experimental work.In addition, we are also not aware of any comparable preexisting simulative study.
The simulation results furthermore suggest that in FC operation degradation of GDL and MPL can remain unnoticeable for a certain period of time, until at some point an abrupt increase in liquid water saturation and therefore cell flooding may occur.Such a highly nonlinear degradation behavior poses a particular challenge for the PEMFC technology, as it can potentially lead to a sudden, unexpected failure during operation.This finding might be also helpful to experimentalists for developing new degradation strategies, as well as for interpreting the results of such degradation tests with respect to the observed effects of GDL aging on cell performance.Simulation experts, on the other hand, could be inspired to study FC degradation on the cell-level by utilizing effective transport parameters as derived in this work.

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.

Fig. 1 .
Fig. 1.(a) The GDL can be differentiated into a bulk region of lower porosity and two surface regions with steep porosity gradients.Being coated onto one GDL side, the MPL penetrates into the large GDL surface pores, resulting in a transition zone between GDL and MPL.For this comparison of porosity profiles the uncompressed microstructures of the impregnated H14 and the impregnated H14+MPL are binarized using the same grayscale threshold.With this simplified binarization the MPL porosity serves only for qualitative analysis at this point.(b) Numerical LBM setup for liquid water intrusion into aged GDL/MPL microstructures: Prior to simulation, the transition zone (see (a)) is cut off the compressed GDL microstructures, reducing the domain size to 500 × 200 × 101 voxels.Subsequently, a single reconstructed MPL layer is appended to the bottom of the cut GDL.The combined GDL/MPL microstructure is then sandwiched between a liquid and a gas phase buffer zone.Liquid water invasion is enforced via a flow rate  liq at the bottom domain boundary, at the top boundary gas phase pressure  gas is kept constant.

Fig. 2 .
Fig. 2. (a)-(e): Modeling of PTFE loss as primary degradation path in GDL microstructures: Impregnated fiber reconstructions are artificially aged by voxel-wise removal of additive from the solid surface.This modeling approach results in (a) reduction of surface hydrophobicity and (b) porosity increase.The degradation procedure is illustrated in (c)-(e).(f)-(g): Modeling of increase in macroporosity  macro MPL as primary degradation effect in MPL microstructures: Pristine and degraded MPL microstructures are obtained with binarization by manual thresholding from high-resolution μ-CT (0.96 μm per pixel) of an impregnated H14+MPL.The degradation procedure is illustrated in (f) and (g).

Fig. 3 .Fig. 4 .
Fig. 3. Sensitivity study on liquid water transport characteristics as a function of the imposed liquid phase flow rate  liq .For the pristine GDL ( = 42.8%,Fig. 2(c)) and pristine MPL ( macro MPL = 20%, Fig. 2(f)) breakthrough is found to occur on (a) different time scales  BT and (b) with varying saturation profiles  local  .Liquid water invasion paths are illustrated for three different liquid phase flow rates in (c)-(e).The color code for the respective phases is as follows: Carbon fiber (gray), PTFE (green), 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.)

Fig. 5 .
Fig. 5. Impact of MPL degradation on liquid water transport characteristics at breakthrough.Increase of MPL macroporosity  macro MPL leads to higher (a) local and (b) overall saturations within the pristine GDL ( = 42.8%,Fig. 2(c)).Liquid water distribution within the GDL is illustrated for three stages of proceeding MPL degradation in (c)-(e).The color code for the respective phases is as follows: Carbon fiber (gray), PTFE (green), 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.)

Fig. 7 .Fig. 8 .
Fig. 7. Impact of GDL degradation due to PTFE loss on liquid water transport characteristics at breakthrough: For a partially aged MPL microstructure ( macro MPL = 60%, Fig. 2(g)) reduction of surface hydrophobicity and porosity increase are now observed to affect the breakthrough characteristics.This results in a significant increase of (a) local (b) and overall saturations.Liquid water distributions are illustrated for three stages of proceeding GDL aging in (c)-(e).The color code for the respective phases is as follows: Carbon fiber (gray), PTFE (green), 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.)

Table 1
Experimental conditions of the μ-CT scanning procedure and grayscale reconstruction process. (a)).

Table 4
Liquid phase velocity and capillary number as employed in the LB simulations of the flow rate sensitivity study.Corresponding physical velocities are obtained with conversion factors of   = 0.96e −6 m/l.u. and   = 4.22e −7 s/l.t (fuel cell operation temperature of  = 80 • C) for the lattice spacing and time, respectively.The capillary number is determined via Ca =  liq  liq  liq ∕ and with the fluid properties as given in Section 4.2.The current density is calculated according to Faraday's law.