LBM simulation of water removal and heat accumulation in PEMFCs: influence of GDL multilayering and determinants of dynamic reaction rates

Proton exchange membrane fuel cells (PEMFCs) electrochemically convert hydrogen and oxygen into water, generating electricity and heat while offering high efficiency, rapid response, and zero-emission operation^[1−4]. PEMFC performance critically depends on effective water and heat management. On the one hand, insufficient water reduces proton conductivity in the membrane, whereas excess water condenses and accumulates in the catalytic layer (CL), gas diffusion layer (GDL), and gas channel (GC), thereby hindering oxygen transport^[3−5]. On the other hand, low temperatures slow electrochemical kinetics, whereas high temperatures accelerate membrane drying and material degradation^[6,7]. Therefore, stable water and heat conditions are essential for sustaining efficient dynamic electrochemical performance in PEMFCs.

Because the oxygen reduction reaction (ORR) occurs at the cathode and its kinetics are slower than those at the anode, the cathode is more prone to non-uniform water transport, which can induce oxygen transport limitations and heat accumulation. As a result, coupled water, oxygen, and heat transport in the cathode GDL often imposes a dominant constraint on performance improvement^[8,9]. Prior studies indicate that optimizing the cathode GDL structure enhances water and heat transport performance^[9−14]. For example, it has been reported by Chu et al.^[9] that variations in GDL porosity strongly affect polarization behavior at high current densities and that higher GDL porosity enhances mass transport. An increase in effective thermal conductivity with decreasing GDL porosity has also been demonstrated by Zamel et al.^[10] through computational fluid dynamics (CFD) simulations. Meanwhile, variations in liquid water content in GDLs also affect overall heat transfer performance^[11]. GDL multilayering with gradient structures has gained increasing attention as a strategy to further optimize local transport relative to single-layer GDLs^[12−14]. Water distribution in multilayer GDLs with a through-plane (TP) pore gradient was analyzed by Ko et al.^[12] using X-ray visualization, and it was found that water content decreases as the pore gradient increases. A multilayer GDL with tapered porosity was proposed by Hao et al.^[13] to enhance the overall heat transfer capacity of PEMFCs. Multilayer GDLs with gradient pores were fabricated by Balakrishnan et al.^[14] using electrospinning, which improved water removal and heat dissipation efficiency, thereby enhancing high current density performance. Although numerous experiments and CFD simulations have shown that multilayer GDL structures enhance overall PEMFC performance, the microscopic mechanisms governing water and heat transport and their effects on the reaction rate remain poorly understood. Accordingly, the effects of non-uniform GDL structural characteristics on multicomponent reactive transport require further investigation.

As a carbon-based porous medium, the GDL exhibits microstructural heterogeneity that significantly influences reactive transport in PEMFCs. The lattice Boltzmann method (LBM) is a mesoscale approach well suited to handling the non-uniform structural characteristics of the GDL and to elucidating the microscopic mechanisms that govern multicomponent reactive transport in PEMFCs^[15−20]. For example, LBM was coupled with deep learning by Wang et al.^[15] to enable large-scale simulations and revealed multiscale water clustering and transport mechanisms across extensive dry and flooded regions in GDLs. A mesoscale model was employed by Lai & Qu^[16] to analyze the effects of GDL porosity, fiber diameter, thickness, wettability, and pressure drop on water dynamics, showing that GDL porosity exerts the greatest influence on water transport. LBM-based models were employed by Habiballahi et al.^[17] and Lai & Qu^[18] to quantify the local TP water distribution and content in multilayer GDLs with gradient porosity, thereby further elucidating the microscopic mechanisms that enhance water removal. However, most mesoscale models simulate water generation assuming either constant-pressure or constant-velocity boundary conditions and often neglect the effects of temperature and oxygen concentration variations on dynamic reaction rates. A LBM model was recently developed by Chen et al.^[19] to relate the dynamic reaction rate to oxygen concentration in PEMFCs. Similarly, a LBM with reactive boundary conditions was used by Guo et al.^[20] to model water and oxygen transport and investigate how gradient porosity affects the reaction rate. However, these studies focused primarily on reactant mass transport and imposed a simplified thermal condition in the computational domain by assigning temperature as either a constant value or a constant temperature difference, without considering reaction heat generation and the associated heat transport during the dynamic reaction process. Non-uniform GDL structural variations modulate water flow, oxygen transport, and heat transfer, which are interrelated and collectively govern dynamic reaction rates. However, prior studies have not clarified the microscopic mechanisms underlying these variations when comparing single-layer and multilayer GDLs, and the combined effects on dynamic reaction rates remain poorly understood. Therefore, it is essential to develop a mesoscale model that captures the microscopic coupling of water and heat transport in single-layer and multilayer GDLs and enables a decoupling analysis to identify the determinants of dynamic reaction rates in PEMFCs. Despite its importance, such analyses remain relatively limited, motivating the present study.

To address these challenges, this study develops an LBM-based mesoscale model in which the electrochemical reaction rate is jointly governed by the oxygen concentration and a temperature-dependent rate constant in PEMFCs and investigates the microscopic mechanisms of water removal and heat accumulation in single-layer and multilayer GDLs across different porosities. Particular attention is paid to elucidating how the coupled processes of water removal and heat accumulation within the GDL affect the dynamic reaction rate, and a decoupled analysis is conducted to quantitatively identify the corresponding determinants. By means of theoretical analyses and numerical testing, this work rigorously confirms the accuracy and reliability of the proposed model. Through detailed examination, this work comprehensively provides a decoupling analysis of the microscopic mechanisms governing water removal and heat accumulation in single-layer and multilayer GDLs and quantifies the impact on the dynamic reaction rate, thereby offering targeted guidance for GDL structural design in PEMFCs.

In this model, the coupling between flow and heat/mass transfer is assumed to be unidirectional, with the numerically computed flow field directly incorporated into the heat and mass transfer models^[21]. In addition, given that the temperature range considered in this work is relatively low and the maximum ratio of the latent heat of water vaporization to the overall reaction enthalpy is about 15% according to the NIST Chemistry WebBook^[22], evaporation and latent heat are considered to have limited effects and are neglected, while water is assumed to remain in the liquid state, consistent with the assumptions adopted in previous studies^[15−20,23]. The multicomponent flow in the PEMFC is described by the Navier–Stokes equations^[8,21]:

where, $ \boldsymbol{u} $ is the velocity (m/s); $ \rho $ is the density (kg/m³); $ p $ is the pressure (Pa); $ \mu $ is the viscosity (Pa$ \cdot $s); and $ {F}_{S} $ is the force term (Pa/m).

The corresponding equation for heat transfer is expressed as follows^[21]:

where, $ \lambda $ is the thermal conductivity; $ {c}_{p} $ is the specific heat at constant pressure; $ T $ is the temperature. The reaction heat source term $ {S}_{T} $ is considered and defined as follows^[21,23]:

where, $ \Delta {H}_{T} $ denotes the reaction heat enthalpy, representing the portion of the total reaction enthalpy $ \Delta H $ that is released as heat. Previous studies have reported that approximately 50% of the energy released during PEMFC operation is dissipated as heat^[24], and the value of $ \Delta {H}_{T} $ is set accordingly in this work. $ {S}^{\alpha } $ is the volumetric source term of the electrochemical reaction, and $ \alpha $ denotes the fluid component.

Owing to the significantly slower diffusion of water relative to oxygen in PEMFCs, the mass transfer model accounts only for oxygen transport^[19,20]. The corresponding governing equation is expressed as follows^[20,21]:

where, $ C $ is the oxygen concentration (mol/m³); $ D $ is the diffusion coefficient (m²/s).

Electrochemical reaction model

In the PEMFC, the electrochemical reaction considered in this work involves the reaction of hydrogen and oxygen to produce water, electricity, and heat, with Pt/C serving as the catalyst. This reaction is expressed as follows^[2]:

The cathodic ORR exhibits slower kinetics than the anodic reaction and is more strongly affected by non-uniform water distribution, oxygen transport, and heat transfer limitations. Therefore, this work focuses on the cathode process and specifically resolves mesoscale multicomponent mass and heat transport in the GDL; therefore, an overall electrochemical reaction expression is adopted as follows^[2,19−21]:

Accordingly, the dynamic electrochemical reaction is modeled as a function of both oxygen concentration and temperature in this work. This treatment provides a more comprehensive description than previous approaches that either employed constant-pressure or constant-velocity boundary conditions^[15−18] or considered reaction rate variations driven solely by oxygen concentration while neglecting temperature effects^[19,20]. The dynamic oxygen reaction rate $ {J}^{oxygen} $ is expressed as the product of the reaction rate constant $ {k}_{elec} $ and the oxygen concentration $ {C}^{oxygen} $, as follows:

where, $ {k}_{elec}\left(T\right) $ is calculated using the Butler–Volmer equation^[25]; $ {i}_{0} $ denotes the exchange current density, which varies with temperature^[26]; $ F $ is the Faraday constant; $ {C}^{ref} $ is the reference oxygen concentration; $ {a}_{e} $ is the transfer coefficient; $ R $ is the universal gas constant; and $ \eta $ is the overpotential. The water dynamic reaction rate, $ {J}^{water} $, is determined from the $ {J}^{oxygen} $, and is expressed as follows:

Accordingly, the expression for $ {S}^{\alpha }\left(x,t\right) $ used in Eqs. (4) and (5) is calculated from the surface flux as follows^[19,20]:

where, $ V $ denotes the total crack volume in the microporous layer (MPL); $ A $ represents the effective reaction area; $ M $ is the molar mass. Accordingly, dynamic reactive transport in the PEMFC is modeled using a LBM-based framework by incorporating these source terms into the evolution equation.

LBM solving model

LBM solving for multicomponent flow

Due to the superior capability of the LBM in handling complex fluid-solid interfaces and multicomponent flow interactions^[22,27−30], this study integrates the LBM with the Shan–Chen model to simulate multicomponent flow in PEMFCs. The evolution equation applies the multiple-relaxation-time (MRT) LBM model to enhance computational stability^[29,30], and is expressed as follows:

where, $ f_{k}^{a}\left(x,t\right) $ denotes the discrete density distribution function at position $ x $ and time $ t $; the superscripts $ a $ and $ eq $ refer to different fluids and the equilibrium state, respectively.

The D2Q9 model is employed in this work, where $ {\overrightarrow{e}}_{k} $ denotes the discrete velocity vectors in nine directions, as illustrated in Fig. 1. $ {M}^{a} $ denotes the orthogonal transformation matrix^[30], defined as:

Figure 1. 

Discrete velocity vectors of the D2Q9 model.

where, $ c $ is the lattice velocity; and $ {\mathit{\Lambda }}^{a} $ denotes the diagonal relaxation matrix, defined as follows^[27]:

where, $ \lambda _{0}^{a}=1.0 $, $ \lambda _{1}^{a}=1.4 $, $ \lambda _{2}^{a}=1.4 $, $ \lambda _{3}^{a}=1.0 $, $ \lambda _{4}^{a}=1.2 $, $ \lambda _{5}^{a}=1.0 $, $ \lambda _{6}^{a}=1.2 $, $ \lambda _{7}^{a}=\lambda _{8}^{a}=\lambda _{\rho }^{a} $, and this set of relaxation parameters has been reported to provide satisfactory numerical stability in prior studies^[20], and $ \lambda _{\rho }^{a} $ is calculated from the kinematic viscosity $ {\upsilon }^{a} $ as follows^[27]:

where, $ {c}_{s} $ denotes the lattice sound speed and is defined as $ {c}_{s}=c/\sqrt{3} $, with $ c={\Delta }x/{\Delta }t $ determined by the prescribed lattice spacing $ {\Delta }x $ and time step $ {\Delta }t $^[21]. Accordingly, $ {c}_{s} $ is set to $ 1/\sqrt{3} $ in this work.

In Eq. (12), $ f_{k}^{a,eq}\left(x,t\right) $ represents the equilibrium distribution function, defined as^[29]:

where, $ {w}_{k} $ denotes the weighting factor; $ {\rho }^{a} $ is the fluid density; and $ {\overrightarrow{u}}^{eq}\left(x,t\right) $ represents the equilibrium velocity, which is defined as follows^[30,31]:

The force term $ F_{k}^{a} $ in Eq. (14) can be expressed as^[31]:

where, $ {\overrightarrow{F}}^{a}\left(x,t\right) $ is calculated using the Shan–Chen model^[29], as expressed below:

where, $ \overrightarrow{F}_{ff}^{a}\left(x,t\right) $ denotes the interfacial tension between different fluids, and $ \overrightarrow{F}_{fw}^{a}\left(x,t\right) $ represents the fluid-solid interaction force that characterizes the wettability of pore surfaces, defined as^[29]:

where, $ {G}^{a\beta } $ denotes the interaction strength between different fluids; $ G_{a}^{as} $ represents the interaction strength between the fluid and solid phases; $ s\left(x\right) $ is a function that distinguishes solid boundaries from pore regions; $ {\psi }_{a}\left(x,t\right) $ is the pseudopotential function dependent on the local density^[29], defined as:

In Eq. (14), $ S_{k}^{a}\left(x,t\right) $ denotes the source term associated with the electrochemical reaction^[19,20], and is defined as:

LBM solving for heat transport

In the heat transport simulation, the D2Q9 model is simplified to the D2Q5 model because previous studies have shown that reducing the number of discrete velocity directions does not compromise numerical accuracy^[21]. The temperature evolution equation is then given by^[21]:

where, $ {h}_{k}\left(x,t\right) $ denotes the discrete temperature distribution function with $ k $ ranging from 0 to 4; and $ h_{k}^{eq}\left(x,t\right) $ represents the corresponding equilibrium distribution function^[21]:

where, $ {K}^{T} $ can be selected between 0 and 1 based on the thermal diffusivity $ {\alpha }^{a} $; $ T $ is the temperature; and $ J_{k}^{T} $ is the rest fraction^[32], defined as:

where, $ J_{0}^{T} $ is set to 0.2; and the temperature $ T $ can be expressed as^[21]:

In Eq. (30), $ {\tau }_{h} $ denotes the relaxation time, and its relationship with the thermal diffusivity $ {\alpha }^{a} $ for different components, including water, oxygen, and solid, is expressed as follows^[21]:

LBM solving for mass transfer

In the mass transfer simulations, the D2Q5 model is likewise employed as described in Eq. (30), and the concentration evolution equation is given by^[33]:

where, $ g_{k}^{a}\left(x,t\right) $ denotes the discrete concentration distribution function with $ k $ ranging from 0 to 4; and $ g_{k}^{a,eq}\left(x,t\right) $ represents the corresponding equilibrium concentration distribution function^[32,33]:

where, $ {K}^{c} $ is a constant, set to 0.5; $ {C}^{a} $ is the concentration; $ J_{k}^{c} $ is the rest fraction, calculated similarly to Eq. (30); and the concentration $ {C}^{a} $ is given by^[33]:

In Eq. (35), $ \tau _{g}^{a} $ denotes the relaxation time, which is related to the diffusion coefficient $ {D}^{a} $ by the following equation^[33]:

where, $ {C}_{Q} $ denotes the lattice correlation coefficient and is set to 0.5. The constant oxygen concentration boundary at the GC top ($ {C}^{oxygen,top} $) is calculated from the partial pressure of oxygen in air and determined to be 7.25 mol m^−3[20]:

Model validation

A computer simulation code was implemented in C++ on a GPU-based parallel computing framework, following Eqs. (1)−(39). Accordingly, numerical simulations were performed to validate the model, including diffusion and reactive transport, heat transport, spurious velocity, a Laplace law bubble test, and fluid-solid contact angle variation, with the detailed results shown in Supplementary File 1 (Figs. S1−S5).

Model parameters, initial and boundary conditions

Figure 2a presents a schematic of PEMFC operation and highlights the mechanisms underlying reaction-generated water and heat. As water is generated and oxygen is consumed at the cathode, the model exclusively focuses on this side. The dynamic reaction rate in the CL is determined by both oxygen concentration and temperature, as defined in Eq. (6), while oxygen diffusing from the GC to the CL is progressively consumed there. Simultaneously, water and heat generated in the CL migrate through the GDL toward the GC. Therefore, the two-dimensional TP GDL effectively captures the dominant relative variation trends in the coupled transport of water, oxygen, and heat, which have also been validated in many experimental studies^[34,35]. Accordingly, the computational domain (Fig. 2b) in this work is defined as a two-dimensional TP slice aligned with the gas flow direction beneath the GC centerline and mainly focuses on the relative variations in the coupled transport of water, oxygen, and heat within the GDL. The domain measures 2,002 μm × 411 μm and is discretized into a uniform grid of 2,002 × 411 lattice cells, comprising a GC (2,002 μm × 252 μm), a GDL (2,002 μm × 150 μm), an MPL (2,002 μm × 8 μm), and a simplified CL (2,002 μm × 1 μm). The electrode contact angle is set to 120°, indicating hydrophobicity. Given that MPL micropores (tens of nanometers) are significantly smaller than crack widths (micrometer scale)^[36], the microporous structure is simplified, whereas crack structures are explicitly represented. The MPL crack depth and width are uniformly set to 4 and 10 μm, respectively^[20], and the CL is modeled as a simplified layer beneath the MPL. Since the conversion between the physical time and the lattice time step is given by $ {\Delta }{t}_{phys}={v}_{LBM}\dfrac{{\Delta }x_{phys}^{2}}{{v}_{phys}} $, the physical time $ {\Delta }{t}_{phys} $ in this work is 4.70 × 10⁻⁸ s. At the bottom of the computational domain, the reaction rate at each lattice node within the CL region is calculated based on the local oxygen concentration and temperature. The reaction rates from all CL lattice nodes are then summed and uniformly distributed throughout the crack volume of the MPL, as shown in Eqs. (10) and (11). During this process, the water and heat generated within the MPL cracks are gradually transported into the GDL. This work also accounts for the distinct oxygen diffusivities in the pore space and in water, thereby representing water blockage effects, with the diffusivity values listed in Table 1. This treatment follows an approach similar to that proposed by Chen et al.^[19,20]. Periodic boundary conditions are applied at the left and right boundaries for water, oxygen, and heat transport. At the top boundary, the computational domain includes only a near-GDL portion of the GC (252 μm), whereas practical GCs are typically on the order of 1 mm in height^[37]. Accordingly, the far-field oxygen concentration $ {C}^{oxygen,top} $ and temperature $ {T}^{top} $ are treated as constants at the top boundary, consistent with prior mesoscale studies^[19,20,38]. Solid surfaces are treated with a non-slip bounce-back condition^[39]. The key simulation parameters are summarized in Table 1.

Figure 2. 

(a) Schematic of PEMFC operation; (b) structure of the computational domain.

A stochastic reconstruction algorithm is employed to generate the GDL pore structure, in which an initial random seed is specified, and carbon fibers are randomly and progressively deposited until the prescribed global porosity is reached. The reconstructed geometry is then converted into a binary structure consisting of solid and pore phases for mesoscale simulations^[42]. Table 2 summarizes the structural parameters used across all simulation cases. These parameters include a single-layer GDL with varying porosity (0.600−0.900, Cases 1−7). Using the porosity of Case 5 as the baseline, multilayer GDLs are then examined by dividing the GDL along the TP direction into three layers of uniform thickness, each with a thickness of 50 μm, denoted as the bottom, middle, and top layers from the GDL-MPL interface to the GDL-CL interface. The TP porosity gradients across the bottom, middle, and top layers are then varied in Cases 8–15, with interlayer porosity differences of −0.100, −0.075, −0.050, −0.025, 0.025, 0.050, 0.075, and 0.100. In all cases, the fiber diameter is sampled from the same normal distribution and constrained to 7–10 μm, consistent with values reported in the literature^[5]. This work further quantifies the pore size distribution (Fig. 3a) and tortuosity (Fig. 3b) for all considered cases. Consistent with extensive experimental evidence^[43], increasing GDL porosity leads to larger overall pore sizes and reduced tortuosity. Figure 3c and d further illustrate the TP porosity distributions corresponding to the different GDL configurations. Meanwhile, a representative two-dimensional GDL simulation yields an effective thermal conductivity of 1.57 W m⁻¹ K⁻¹, compared with 1.8 W m⁻¹ K⁻¹ obtained from three-dimensional experiments with the same porosity and thickness^[35]. The discrepancy is mainly attributed to fewer fiber contacts in two-dimensional representations than in real three-dimensional microstructures, and the deviation remains acceptable. For cross-case comparison, the transient results are evaluated using the same steady-state criteria as those used in prior studies^[20,38], with steady-state reached when water reaches the GC top boundary or when the temporal variation in GDL water saturation becomes sufficiently small.

Figure 3. 

(a) Pore size distribution, and (b) tortuosity for different cases; TP porosity distributions for (c) single-layer GDLs, and (d) multilayer GDLs with varying porosities.

Figure 4 illustrates the water flow characteristics in single-layer GDLs with varying porosity. During dynamic electrochemical reactions, water progressively coalesces into clusters within MPL cracks. Driven by capillary forces, water infiltrates the GDL pores and exhibits three characteristic modes: droplet formation, cluster development, and flow elongation^[44]. Subsequently, water preferentially breaks through pores with minimal capillary resistance, thereby establishing stable transport paths^[45], which are herein termed the dominant water removal paths within the GDL. For instance, in Region A₁ of Fig. 4a, although the relatively low local porosity hinders water transport, the coalescence of initial water clusters enhances the driving force and generates finger-like flow paths toward the GDL-GC interface, which is consistent with flow visualization experiments and thereby supports the reliability of the simulation results^[46]. Conversely, in Region A₂, extremely low local porosity poses a high capillary barrier, thereby preventing water breakthrough. In particular, when GDL porosity falls below 0.700 (Fig. 4b), excessively high tortuosity promotes the development of severe flooding within the GDL (Regions B₁–B₃), markedly hindering oxygen transport. As GDL porosity increases, overall capillary resistance decreases substantially, promoting the development of dominant water removal paths and thus strengthening oxygen supply (Regions B₄−B₆). Meanwhile, the coalescence of breakthrough droplets enhances liquid connectivity within the GDL and forms more continuous water removal paths (Regions C₁−C₃)^[47], whereas the resulting larger droplets in the GC are more difficult to transport by the gas flow. When the porosity reaches 0.9 (Case 7), weakened capillary effects suppress the development of dominant water removal paths (Region B₇) and intensify localized water blockage. Therefore, for single-layer GDLs with varying porosity, capillary resistance at excessively low porosity, pore connectivity and tortuosity, and weakened capillary effects at excessively high porosity jointly influence the water removal characteristics.

Figure 4. 

Evolution of water flow characteristics in (a) a single-layer GDL with porosity 0.800 and (b) Cases 1–7 with varying single-layer GDL porosities.

Figure 5a depicts the temporal evolution of the local water saturation (S_local,y) along the TP direction in the GDL for Case 5. From the MPL-GDL interface to the GDL-GC interface, S_local,y decreases gradually with increasing y, consistent with previous reports^[18,44]. Specifically, near the GDL-GC interface, structural discontinuities and capillary-pressure jumps induce an increase in S_local,y. As porosity increases (Fig. 5b), S_local,y decreases in the MPL-GDL region (Region A₁) but increases in the GDL-GC region (Region A₂). This occurs because higher GDL porosity weakens flow blockage and enhances permeability, thereby progressively mitigating water accumulation and lowering the average GDL pressure (Pre_GDL) (Fig. 5c). The GDL water saturation (S_GDL) correspondingly decreases with increasing porosity (Fig. 5d), but the trend reverses in Case 7. The water removal efficiency (E_water) is defined as the ratio of water saturation in the GC to that in the GDL (S_GC/S_GDL) and is used as a dimensionless drainage indicator that characterizes the relative tendency of water to migrate across the GDL-GC interface into the GC relative to the amount retained in the GDL at the same time, thereby assessing water transport and breakthrough ability. For Cases 1−7, the corresponding E_water values are 0.25, 0.32, 0.44, 0.55, 0.77, 1.01, and 0.86, respectively (Fig. 5e). In this work, each case was regenerated at least five times using a random algorithm to obtain comparable GDL configurations with identical porosity but distinct microstructures, and the maximum deviation among the replicated E_water values was 6.8%, confirming the numerical reproducibility of the simulations. The trend in E_water closely parallels that of the average normalized oxygen concentration in the reaction region (C^oxygen/C^oxygen,top), wherein increasing the porosity of single-layer GDLs enhances water removal efficiency by up to a factor of four and improves oxygen transport capacity by approximately 34%, whereas excessive porosity leads to a decline in efficiency.

Figure 5. 

(a) Temporal evolution of S_local,y along the TP direction for a GDL porosity of 0.800 (Case 5); (b) steady-state S_local,y distributions; (c) Pre_GDL, (d) S_GDL; and (e) E_water and normalized C^oxygen for single-layer GDLs with varying porosities.

Water removal in multilayer GDLs

To evaluate the performance of multilayer GDLs relative to the single-layer case, the GDL is further modeled as a three-layer porosity gradient structure comprising bottom, middle, and top layers along the TP direction. Accordingly, Case 5 (P_TP: 0.800) serves as the baseline. As illustrated in Fig. 6, when the porosity gradient decreases along the TP direction, flow resistance is minimal near the MPL and increases progressively toward the GC. This configuration creates a capillary barrier at the GDL-GC interface, hindering the development of dominant water removal paths (Regions A₁−A₄). Conversely, as the porosity gradient moderately increases along the TP direction (Cases 12, 13), the higher driving force near the GDL-MPL interface offsets the larger capillary resistance, whereas the lower capillary resistance near the GDL-GC interface promotes water breakthrough, thereby markedly reducing water accumulation (Regions A₆−A₇). Correspondingly, the enhanced water mobility increases the number of initial water clusters from six (Case 5) to eight (Case 13). However, when the interlayer porosity difference becomes excessive (Cases 14, 15), the resulting mismatch in capillary entry pressure across adjacent layers disrupts dominant water removal paths and promotes water retention (Regions A₈−A₉). Consequently, water removal is mainly influenced by water invasion near the GDL-MPL interface, breakthrough resistance near the GDL-GC interface, capillary-driven transport assisted by the porosity gradient, and retention induced by excessive interlayer contrast for multilayer GDLs with different porosity gradients.

Figure 6. 

Illustration of water flow evolution under different multilayer GDL porosities.

Therefore, multilayer GDLs can effectively modulate interlayer water flow and distribution characteristics (Fig. 7a−d). When the porosity gradient decreases along the TP direction (Fig. 7a), the S_local,y gradient steepens, and most water is concentrated in the GDL bottom layer. As the porosity gradient increases along the TP direction (Fig. 7b), the S_local,y gradient from bottom to top becomes less steep, indicating the formation of a stable driving force pathway for water flow. However, an excessive interlayer porosity difference (Fig. 7c) generates high interlayer capillary resistance, thereby increasing bottom-layer water content and steepening the S_local,y profile again. Therefore, appropriately increasing the porosity gradient along the TP direction helps establish a capillary-pressure gradient (Fig. 7e), thereby reducing water accumulation in the GDL, yielding S_GDL values of 0.57 and 0.53 (Cases 12, 13). In particular, performance is optimal in Case 13, with E_water increasing by 40% and normalized C^oxygen by 11% relative to the single-layer case (Fig. 7g). However, excessive interlayer porosity differences reduce both E_water and normalized C^oxygen (Cases 14, 15), and the effect is especially pronounced in Case 15, where performance falls below that of the single-layer case.

Figure 7. 

Temporal evolution of S_local,y along the TP direction for multilayer GDLs in (a) Case 8, (b) Case 13, and (c) Case 15; (d) steady-state S_local,y, (e) S_GDL, (f) Pre_GDL, and (g) E_water and normalized C^oxygen for multilayer GDLs with varying porosities.

Heat accumulation with GDL multilayering

Heat accumulation in single-layer GDLs

Dynamic electrochemical reactions in PEMFCs generate substantial water and heat^[6−8,48]. Because water has a high specific heat capacity and carbon fiber exhibits high thermal conductivity^[10,11], spatial heterogeneity in water content and GDL porosity significantly influences heat transfer in PEMFCs. As illustrated in Fig. 8a, a temperature gradient develops within the GDL, with the highest temperature located near the bottom of the GDL, forming a pronounced heat accumulation region that decreases along the TP direction. This trend is also consistent with the findings reported by Wei et al.^[48]. As the reaction proceeds, heat accumulation within the GDL progressively intensifies, with the local heat source primarily concentrated in regions where water transport is restricted (Regions A₁), whereas regions with higher water content exhibit noticeably weaker heat accumulation (Regions A₂), consistent with Fig. 4. Meanwhile, higher porosity in single-layer GDLs reduces local water content and fiber fraction. The concurrent reduction in heat capacity and thermal conduction ability weakens thermal buffering and heat dissipation, thereby promoting a higher temperature rise (Fig. 8b). Notably, in Cases 5, 6 (Regions B₁−B₂), local heat sources are widespread in the bottom of the GDL. Although water content rises at excessively high GDL porosity (Case 7), the carbon fiber fraction decreases further, thereby sustaining pronounced heat accumulation (Region B₃). Figure 8c shows the steady-state average temperatures of the GDL (T_GDL) and its bottom layer (T_Bottom), with the maximum across-case difference in T_GDL reaching 12 K. In this work, each case was regenerated at least five times using a random algorithm to obtain comparable GDL configurations with identical porosity but distinct microstructures, and the resulting T_GDL values varied by no more than 0.36 K, confirming the numerical reproducibility of the simulations.

Figure 8. 

(a) Evolution of temperature distributions in Case 5 with porosity 0.800; steady-state (b) temperature distributions, and (c) T_Bottom and T_GDL for single-layer GDLs with varying porosities.

Heat accumulation in multilayer GDLs

In the temperature fields of multilayer GDLs, as the porosity gradient decreases along the TP direction (Fig. 9a; Cases 8−11), flooding is pronounced, whereas heat accumulation remains relatively low (Regions A₁−A₄). By contrast, when the porosity gradient increases along the TP direction and the interlayer porosity difference is small (Fig. 9b; Cases 12, 13), heat accumulation increases markedly over time, and the localized heat source region expands (Regions B₁−B₂), reaching a maximum at an interlayer porosity difference of 0.050 (Case 13). At this point, T_Bottom is 3.4 K higher than the single-layer value (Fig. 9c). However, when the interlayer porosity difference becomes excessive (Fig. 9b; Cases 14, 15), the temporal increase in heat accumulation weakens, and the extent of localized heat sources gradually shrinks (Regions B₃−B₄). This behavior arises because the GDL bottom layer is the primary site of reaction heat accumulation. The increased bottom layer capillary resistance promotes local water retention (Fig. 9d), while the fiber fraction remains relatively high (Fig. 9e). Together, these factors enhance thermal buffering and heat transfer, thereby significantly reducing the temperature rise (Fig. 9c). Therefore, these results indicate that multilayer GDLs can mitigate heat accumulation relative to single-layer GDLs. Notably, although neglecting evaporation and latent heat may overestimate the magnitude of the temperature rise, it does not affect the interpretation of the relative trends in heat accumulation under different GDL structures. Specifically, the simulation results show that regions with higher water content tend to exhibit lower temperatures, and considering latent heat would further strengthen this trend, which remains consistent with the heat accumulation pattern revealed in this work.

Figure 9. 

Temperature distributions for multilayer GDLs with (a) porosity gradients decreasing and (b) increasing along the TP direction. (c) Steady-state values of T_Bottom and T_GDL, (d) relative water content ratios in the bottom, middle, and top layers and (e) interlayer porosity different among multilayer GDLs with varying porosities.

Correlation between water removal and heat accumulation

Although several studies have investigated the effects of GDL structural variations on water flow or heat transport^[8−14], the coupled mechanisms linking water removal and heat accumulation within the GDL remain insufficiently understood. This study finds that higher porosity increases both E_water and T_GDL in single-layer GDLs (Fig. 10). When porosity becomes excessive (Case 7), E_water decreases, whereas T_GDL remains elevated. This occurs because weakened capillary effects reduce water removal capacity within the GDL, whereas an insufficient fiber fraction markedly diminishes heat dissipation capacity. In multilayer GDLs, Case 5 serves as the reference to partition the domain into Regions A−C. Derived empirical relations for E_water and T_GDL in Region A can readily be used to infer those in the remaining regions. When the porosity gradient decreases along the TP direction (Cases 8−11), E_water and T_GDL are lower and cluster in Region A. When the porosity gradient increases along the TP direction (Cases 12, 13), both E_water and T_GDL are higher and fall primarily in Region B until interlayer porosity differences become excessive (Cases 14, 15). Different from the single-layer case, T_GDL is markedly reduced beyond the turning point (Region C). This occurs because both bottom-layer water content and fiber content increase, thereby enhancing heat dissipation performance. Therefore, the multilayer design achieves a lower degree of heat accumulation and provides relatively improved water removal performance within the GDL compared with the single-layer case. Notably, this work does not aim to directly predict the absolute optimal parameters for practical PEMFCs, but rather to support the comparative interpretation of the trends and the mechanistic conclusions identified in this work.

Figure 10. 

Correlation between E_water and T_GDL in the single-layer and multilayer GDLs.

Insights into determinants of dynamic reaction rates

In PEMFCs, the optimization effect of GDL structural regulation on reaction performance can be effectively evaluated by quantifying variations in $ {J}^{oxygen} $. In single-layer GDLs (Fig. 11a), $ {J}^{oxygen} $ increases with porosity and reaches a maximum in Case 6, whereas further increasing porosity to 0.9 (Case 7) does not lead to a further increase in $ {J}^{oxygen} $. Across the single-layer cases, $ {J}^{oxygen} $ increases by up to threefold. In multilayer GDLs (Fig. 11b), $ {J}^{oxygen} $ is generally lower when the porosity gradient decreases along the TP direction (Cases 8−11). By contrast, increasing the TP porosity gradient markedly improves reaction performance (Cases 12, 13). At an interlayer porosity difference of 0.050 (Case 13), $ {J}^{oxygen} $ reaches its maximum, which is approximately 26% higher than that of the single-layer GDL. However, excessively large interlayer porosity differences (Cases 14, 15) markedly deteriorate reaction performance. For example, at an interlayer porosity difference of 0.100 (Case 15), $ {J}^{oxygen} $ is 16% below the single-layer value. Therefore, variations in porosity in both single-layer and multilayer GDLs have a significant impact on the reaction rate, and the governing factors require further analysis.

Figure 11. 

Variation of $ {J}^{oxygen} $ in the (a) single-layer and (b) multilayer GDLs.

As shown in Eq. (8), $ {J}^{oxygen} $ is governed by the temperature-dependent rate constant $ {k}_{elec}\left(T\right) $ and the oxygen concentration $ {C}^{oxygen} $. Variations in the GDL structure simultaneously modify water and heat transport, which in turn alter $ {C}^{oxygen} $ and $ {k}_{elec}\left(T\right) $ and thereby regulate $ {J}^{oxygen} $. This study takes Case 5 as a unified reference and analyzes the steady-state relative variations of $ {k}_{elec}\left(T\right) $, $ {C}^{oxygen} $, and $ {J}^{oxygen} $ under different GDL structures, thereby quantitatively assessing their relative contributions and identifying which factor acts as the determinant of the $ {J}^{oxygen} $ variation, namely temperature-induced changes in $ {k}_{elec}\left(T\right) $ or variations in $ {C}^{oxygen} $. As shown in Fig. 12, temperature in single-layer GDLs varies markedly with porosity in Cases 1−6, causing $ {k}_{elec}\left(T\right) $ to increase by 32%, whereas $ {C}^{oxygen} $ increases by only 18%, indicating that the improvement in $ {J}^{oxygen} $ is mainly determined by $ {k}_{elec}\left(T\right) $ (Fig. 12a). However, an excessively high porosity (Case 7) weakens capillary effects, impeding water and oxygen transport, leading to a 5% decrease in $ {C}^{oxygen} $. Meanwhile, the reduced fiber content sustains a high degree of heat accumulation (Fig. 8), leaving $ {k}_{elec}\left(T\right) $ nearly unchanged. At this stage, the reduction in $ {J}^{oxygen} $ is primarily determined by the relative decrease in $ {C}^{oxygen} $. In multilayer GDLs, when the porosity gradient decreases along the TP direction (Cases 8−11), a capillary barrier forms that impedes water and oxygen transport, causing $ {C}^{oxygen} $ to vary by only up to 4%, whereas $ {k}_{elec}\left(T\right) $ varies by up to 15%, thereby maintaining the determinant of $ {J}^{oxygen} $ as the temperature-dependent $ {k}_{elec}\left(T\right) $ (Fig. 12b). By contrast, when the porosity gradient increases along the TP direction (Cases 12, 13), a capillary-pressure gradient is established significantly enhancing water removal and oxygen supply, which increases $ {C}^{oxygen} $by up to 11% and correspondingly increases $ {k}_{elec}\left(T\right) $ by 8%. Consequently, the determinant of the improvement in $ {J}^{oxygen} $ shifts from $ {k}_{elec}\left(T\right) $ to $ {C}^{oxygen} $. However, when the interlayer porosity difference becomes excessively large (Cases 14, 15), the high fiber fraction and elevated water content in the bottom layer markedly mitigate heat accumulation (Fig. 9). Under these conditions, $ {C}^{oxygen} $ decreases by 3%, whereas the temperature effect causes $ {k}_{elec}\left(T\right) $ to decrease by 7%, thereby causing the determinant of $ {J}^{oxygen} $ variation reverts to $ {k}_{elec}\left(T\right) $. Therefore, the trade-off between water removal and heat accumulation within the GDL markedly influences the dynamic reaction rate.

Figure 12. 

Effects of $ {C}^{oxygen} $ and $ {k}_{elec}\left(T\right) $ on $ {J}^{oxygen} $ in (a) single-layer and (b) multilayer GDLs. The red regions indicate that $ {k}_{elec}\left(T\right) $ is the determinant of $ {J}^{oxygen} $, whereas the blue regions indicate that $ {C}^{oxygen} $ is the determinant of $ {J}^{oxygen} $.

By quantitatively analyzing the determinants governing reaction rate variations under different GDL structures in Fig. 12, this work further provides a mechanistic understanding of how GDL structural regulation influences the PEMFC reaction rate (Fig. 13), which comprises two principal pathways. One is the oxygen mass transfer pathway, in which GDL structural variations affect water removal efficiency, thereby altering the oxygen concentration, as quantified by $ {C}^{oxygen} $, which in turn strongly influences the reaction rate. The other is the rate constant pathway, in which GDL fiber fraction and local water content jointly govern the degree of reaction heat accumulation, thereby affecting the rate constant, as quantified by $ {k}_{elec}\left(T\right) $, which also exerts a strong influence on the reaction rate. Changes in the reaction rate feed back through water and heat generation fluxes to the water removal and heat accumulation pathways within the GDL, thereby forming a nonlinear closed-loop process. By quantifying the relative variations of $ {C}^{oxygen} $, $ {k}_{elec}\left(T\right) $, and $ {J}^{oxygen} $ for different single-layer and multilayer GDL structures, the effects and relative contributions of the two pathways to the reaction rate can be quantitatively separated, thereby providing targeted guidance for GDL structural design.

Figure 13. 

Schematic of the decoupled analysis linking the GDL structure to the dynamic reaction rate in PEMFCs.

This study develops an LBM-based mesoscale model to investigate the microscopic mechanisms governing water removal and heat accumulation in non-uniform single-layer and multilayer GDLs. The model is further employed to decouple the determinants of the dynamic reaction rate across GDL structural variations and to identify two principal pathways. The main findings are summarized as follows:

(1) Increasing porosity in single-layer GDLs lowers capillary resistance and widens the dominant water removal paths, thereby facilitating water and oxygen transport. However, as the porosity approaches 0.9, weakened capillary effects exacerbate localized water blockage, thereby reducing the water removal performance of the GDL. In multilayer GDLs, an increasing through-plane porosity gradient allows the higher driving force near the GDL-MPL interface to offset the larger capillary resistance, while the lower resistance near the GDL-GC interface promotes water breakthrough. Accordingly, the capillary pressure gradient in multilayer GDLs yields a maximum improvement of 40% in water removal efficiency compared with the single-layer cases. However, when the interlayer porosity difference becomes excessively large, the resulting mismatch in capillary entry pressure across adjacent layers impedes both water removal and oxygen supply.

(2) In single-layer GDLs, higher porosity improves water removal but often intensifies localized heat accumulation because reduced water content and fiber fraction weaken thermal buffering and heat dissipation. In multilayer GDLs, increasing the porosity gradient along the through-plane direction progressively intensifies heat accumulation until the interlayer porosity difference exceeds 0.05. Unlike the single-layer case, heat accumulation is markedly mitigated beyond the turning point. This mitigation arises because increased capillary resistance in the bottom layer promotes local water retention while the fiber fraction remains relatively high. Together, these factors strengthen heat transfer and thermal buffering, thereby mitigating heat accumulation.

(3) The dynamic reaction rate depends jointly on the rate constant affected by heat accumulation and the oxygen concentration linked to water removal efficiency within the GDL. In multilayer GDLs, the maximum reaction rate increases by up to 26% compared with the single-layer case. To elucidate how water removal and heat accumulation affect reaction rate, a decoupled analysis is conducted to identify its determinants. In single-layer GDLs, the determinant of reaction rate is the temperature-dependent rate constant, but shifts to oxygen concentration at excessive porosity. In multilayer GDLs, the determinant of the reaction rate is the temperature-dependent rate constant when the porosity gradient decreases, whereas it shifts to oxygen concentration as the porosity gradient increases, but reverts to the temperature-dependent rate constant when interlayer porosity differences become excessive.