Water Transport in the Gas Diffusion Layer of a Polymer Electrolyte Fuel Cell: Dynamic Pore-Network Modeling

The pore-scale modeling is a powerful tool for increasing our understanding of water transport in the fibrous gas diffusion layer (GDL) of a polymer electrolyte fuel cell (PEFC). In this work, a new dynamic pore-network model for air-water flow in the GDL is developed. It incorporates water vapor transport coupled with liquid water by a phase change model. One important feature of our pore-network model is that a recently developed semi-implicit scheme for the update of water saturation is used. It provides good numerical stability in modeling liquid water transport in the GDL even for very small capillary number values. A number of case studies are conducted to illustrate several important mechanisms of water transport in the GDL, such as cyclic processes of local drainage and imbibition, channeling and capillary-fingering evolutions of water flow pattern, and eruptive water transport. Moreover, we show that liquid water separation in the GDL between the ribs and gas channel (GC) is formed under dry GC condition, which is qualitatively in agreement with in situ experimental results.

Water management plays an important role in the durability and performance of polymer electrolyte fuel cells (PEFCs). A delicate water balance must be maintained to keep the membrane hydrated and at the same time avoid liquid water flooding in the porous components. [1][2][3] Over the past two decades, numerous studies have been conducted for improving our understanding of water transport mechanisms within the PEFC. Compared to the anode side, water transport in the cathode has usually attracted more attention. This is because the cathode is prone to being flooded first, and the oxygen reduction reaction (ORR) is very sluggish. 3 The fibrous gas diffusion layer (GDL), as a key component, plays an important role in removing excessive water from the catalyst layer (CL) to the gas channel (GC). Up to now, much effort has been invested in gaining insights into pore-scale processes of liquid water transport in the GDL. [4][5][6][7][8][9][10] Despite the great success in in situ observations of water distribution, 8,9,[11][12][13][14] water dynamics in the GDL is still not well understood. Also, the cost of those experiments is quite high. This has stimulated many pore-scale numerical studies of the GDL over the past few years. [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33] Several distinctive features of liquid water transport in the GDL make its pore-scale modeling extremely difficult. First, a small flow rate of liquid water, with a capillary number value of around 10 −8 , causes the modeling computationally expensive. Second, the complexity 34 of the pore structure and mixed wettability of the GDL give rise to the difficulty in developing an elegant pore-scale model. Nevertheless, many achievements still have been booked.
Mukherjee et al. 4 were the first to conduct the Lattice-Boltzmann modeling of air-water flow in a stochastically reconstructed nonwoven GDL. The authors illustrated the structure-wettability influence on the underlying water transport and flooding dynamics. But, the computational effort was extremely heavy. Compared to direct simulations, [29][30][31][32][33] the pore-network modeling has attracted more attention, [16][17][18][19][20][21][22][23][24][25][26][27][28] mainly because it is computationally cheap. The basic idea of a pore-network model is to first construct a pore network which can adequately represent the porous medium of interest; then, flow and transport are solved in the constructed network, usually with the help of some local rules. In such a way, the computational effort can be reduced considerably. There have been many invasionpercolation-based quasi-static pore network models used in the GDL studies, [15][16][17][18][19][20][21][22][23][24][25][26] which are very computationally cheap. These models are usually used to provide material properties like capillary pressure and relative permeability curves of the GDL, 15,16,22 which are required in the Darcy-scale modeling. 2,35,36 Also, given the fact that liquid water transport in the GDL is strongly capillarity-dominant, the quasi-static z E-mail: chaozhong.qin@gmail.com; c.qin@uu.nl models are able to predict the flow pattern of liquid water. For example, Lee et al. 25,26 developed a pore-network model combining invasion-percolation path finding and subsequent viscous two-phase flow calculation to investigate steady-state water distributions in the GDL. Although the model is computationally effective, it is worth noting that cyclic processes of local drainage and imbibition may prevail over the whole network; thus a steady-state viscous flow at the pore scale does not exist. This phenomenon will be later illustrated by our dynamic pore network modeling.
Compared to quasi-static models, a dynamic model has the following advantages. First, it can facilitate versatile boundary conditions. Second, it can couple with water vapor transport through a phase change model. Third, it can capture important transport mechanisms of liquid water in the GDL, such as eruptive transport of liquid water near the GDL-GC interface. 7,9 But, the computational effort of a dynamic pore-network model is usually heavy. There are two kinds of dynamic pore-network model, called single-pressure algorithm model 5,6,28 and two-pressure algorithm model. 27,[37][38][39] By simplifying the problem, the single-pressure approach reduces computational effort. But, it has some disadvantages. 37 For example, no local capillary pressure for pore bodies can be defined. This means that no information about the interface behavior under dynamic conditions can be gained. As a result, the phase change between water vapor and liquid water will not be well modelled. Moreover, this approach exhibits some inconsistent behavior in fluids occupancy in the network, particularly for small capillary number values. More information about the comparison between the two approaches can be found in Ref. 37. In the literature, there are a few studies of the GDL with the singlepressure algorithm. Sinha and Wang 28 developed a single-pressure pore-network model to study the dynamics of liquid water transport in the GDL under realistic fuel cell operating conditions. Hinebaugh and Bazylak 5 developed a similar dynamic pore-network model of the GDL. They stochastically placed the condensation locations of water vapor in the GDL to investigate its effect on the profile of throughplane liquid water distribution. Later, Medici and Allen 6 included water vapor and heat equations in their dynamic pore-network model. They studied the effect of different GC conditions such as temperature and relative humidity on the water transport in the GDL. It is worth noting that all above-mentioned single-pressure pore-network models assumed no corner flow of air in the network which resulted in the trapping of air phase. But, recent studies showed that immobile air clusters are not usually formed in the GDL. 40 On the other hand, there are very few studies of the GDL with the two-pressure algorithm. This may be because the two-pressure approach suffers from the issue of numerical stability at very small capillary number values. To the best of our knowledge, Lee et al. 27 were the first to present a Figure 1. SEM images of (left) carbon paper (Toray 060), and (right) carbon cloth (E-Tek). 18 The permission of using the graphs is issued by Elsevier.
two-pressure pore-network model to study the water dynamics in the hydrophobic GDL. But, they used an explicit scheme for updating the water saturation.
Recently, Joekar-Niasar et al. 37 developed a new two-pressure dynamic pore-network model. A semi-implicit scheme for the update of water saturation was used by introducing capillary diffusion. It was reported that this model provides more numerical stability, compared to previous models, particularly for unfavorable viscosity ratios and small capillary number values. In this work, we have employed a similar computational algorithm in developing a dynamic pore-network model of the GDL. The model helps us understand the liquid water transport mechanisms in the GDL. Moreover, we include water vapor transport throughout the network, which is directly coupled with liquid water dynamics through a phase change model. This allows us to study the role of phase change in the water removal process under partially humidified GC conditions.

Pore Network Generation
Carbon paper and carbon cloth, as shown in Fig. 1, are two commonly used GDLs. Their pore structures are very different from geological porous media such as soil and sandstone. 34 In addition, pore throats cannot be easily identified in these highly-porous fibrous materials. Up to now, much effort has been invested in representing a GDL with a proper pore network. 5,18,27,38 For example, Thompson et al. 38 used Voronoi diagrams to construct the pore networks of fibrous porous media. Luo et al. 18 extracted the pore networks from their stochastically generated GDLs. But, no verification has been given that those networks adequately represented a realistic GDL. Alternatively, in the past, most researchers 6,16 19-24 employed a regular lattice-based pore network, which was calibrated primarily by medium permeability and pore-size distribution.
In this work, we focus on the flow dynamics of water in the GDL, rather than provide a better representation of the GDL. Therefore, as in most previous studies, we use a regular pore network with cubic pore bodies and pore throats of square cross sections. This allows us to account for the corner flow of air throughout the whole network; and thus no air can be trapped. The size of a cubic pore body is characterized by the radius of the inscribed sphere. The radii of inscribed spheres in pore bodies are assumed to have a truncated lognormal distribution with mean radius of 11 μm 28 . The inscribed radius of pore throat is determined by the following formula: 37 [2] where R i and R j are the inscribed radii of the two connected pore bodies, d is the constant distance (set equal to 25 μm) between lattice grids, and n is a parameter controlling the ratio between the radii of pore bodies and pore throats. Here, its value is set to 1.0, in order to match the permeability of a carbon-paper GDL. 22 Finally, the resulting distributions of the inscribed radius of pore body, the inscribed radius of pore throat, and the length of pore throat are plotted in Fig. 2.
To reduce computational effort, in most case studies, we use a twodimensional (2D) pore network for representing the cross section of a GDL as shown in Fig. 3. The network has 960 pore bodies and 1680 pore throats including inlet and outlet boundaries. This represents a computational domain of 250 μm by 2000 μm in the X and Y directions, respectively. The network parameters, listed in Table I, are assumed mainly based on the information in Ref. 22 and Ref. 28. It is noted that we do not intent to represent a very specific Toray GDL. But, we attempt to include the main material features based on available literature data. Three realizations (R1, R2, and R3) of the 2D network are generated with the through-plane permeability of around 10 −11 m 2 representing an uncompressed carbon-paper GDL. The compression effect on water transport will be further studied. As to the anisotropy feature of the GDL with respect to intrinsic permeability, there may be two methods to obtain larger in-plane permeability in the pore network generation. One is to increase the in-plane coordinate number (i.e. more pore throats in the in-plane directions); the other is to increase the sizes of in-plane pore throats. However, in this work, for  simplicity, we focus on the through-plane water transport as many previous studies. [23][24][25][26][27][28] This anisotropy feature is not considered in our present network generation. Finally, one case study with a 3D pore network is also done in this work. The information about the 3D network will be given later on.

Dynamic Pore-Network Model for Air-Water Flow and Phase Change
The following main assumptions are made in the development of the dynamic pore-network model.
r The volume of pore throats is negligible compared to the volume of pore bodies, because the total volume of pore throats in the generated network is less than 5% of that of pore bodies. Thus, the time required for filling a pore throat is negligible compared to that of a pore body.
r The hydraulic and diffusive resistances in pore bodies are negligible compared to those in pore throats. This assumption is generally used by most quasi-static and dynamic pore-network models.
r Air and water are incompressible and immiscible. The change of air viscosity due to the change in water vapor concentration is negligible.
r The surface wettability of the GDL is hydrophobic and homogenous (i.e. a constant contact angle is assumed).
For pore body i, the volume conservation of each phase is given as: where V i [L 3 ] is the volume of the pore body, s α i [-] is the α-phase saturation in the pore body, j denotes a neighboring pore body, N i is the coordinate number of the pore body, Q α i j [L 3 /T] is the α-phase flux through pore throat i j, g denotes the wetting phase of air, l denotes the nonwetting phase of liquid water, ρ α is the α-phase mass density, and r α i [M/T] denotes the phase change rate. Its expression will be given later on.
The flux of each phase through pore throat i j is calculated by Hagen-Poiseuille equation: where K α i j [TL 4 /M] is the α-phase conductance in the pore throat, p α i is the α-phase pressure in pore body i, and p α j is the α-phase pressure in pore body j. The capillary pressure in the pore body is related to the phase pressures by: In addition, it is computationally efficient to work with a mixture pressure for obtaining the air and water pressure fields over the pore network. 37 The mixture pressure in a pore body can be defined as: Combining Equations 3-6, and with the constraint, s g i + s l i = 1, the mixture pressure equation is obtained as: Once the mixture pressure field is obtained, the air and water pressures are calculated by (A1 and A2). Then, the flow fields are given by Eq. 4. The water saturation in a pore body can be updated by Eq. 3, with the help of a first-order phase change model: where k [L/T] is the coefficient of the phase change rate between water vapor and liquid water, a gl [L 2 ] is the interfacial area between air and water, C wv i [M/L 3 ] is the mass concentration of water vapor, and C sat i is the saturated water vapor concentration. Finally, the mass conservation of water vapor in pore body i is given as: ) unless CC License in place (see abstract  where V g i is the volume of air phase in the pore body, D wv i j [L 2 /T] is the diffusivity of water vapor, A g i j is the flow area for air in pore throat i j, and L i j is the pore-throat length. In Eq. 9, the l.h.s. denotes the temporal accumulation of water vapor in the air phase. On the r.h.s., the first and second terms represent the advective fluxes out of and into pore body i, respectively; the third term represents the diffusive flux exchanged between pore body i and its surrounding pore bodies; and the last term represents the phase change rate. Usually, the diffusive transport of water vapor is dominant in the GDL. But, for the generality of the model development, we keep the advection terms in Eq. 9. In this pore-network model, we select three primary variables, namely, mixture pressure (p i ), liquid water saturation (s l i ), and water vapor concentration (C wv i ). The numerical strategy is briefly described as follows. In each time step, first, Eq. 7 is solved for the pressure field. Then, liquid water saturation is updated by Eq. 3. Finally, Eq. 9 is solved for the water vapor concentration. For the detail of numerical implementation and all constitutive equations used in this model, one can refer to the Appendix and the work by Joekar-Niasar et al. 37

Description of Numerical Simulations
On the cathode side of a PEFC, the GDL is sandwiched between the micro porous layer (MPL) and the bipolar plate (BP). Thus, proper conditions need to be provided for those boundaries. During some in situ observation experiments, 9,11 little liquid water was seen in the MPL. One possible explanation is that liquid water quickly evaporates at the CL-MPL interface and diffuses through the MPL to the GDL. This is because the fine-textured MPL is highly hydrophobic so that no condensation can take place there. But, at the MPL-GDL interface, water vapor starts to condense into liquid water due to the coarse and less hydrophobic GDL. Therefore, at that boundary, a uniform water vapor flux may be assumed, while at the boundary with the BP, a zero-flux boundary condition is imposed under the ribs, and the GC is assumed to be free of liquid water due to a high air flow rate. Furthermore, it is assumed that the inlet vapor totally condenses into liquid water inside the first layer of the GDL pores. In such a way, we can simplify the complicated condensation process at the MPL-GDL interface.
In detail, the used boundary conditions and numerical treatments are described as follows. The inlet and outlet pore bodies are assumed to be devoid of liquid water with zero mixture pressure. At the bottom of the GDL, no-flux boundary conditions are specified for air flow, liquid water flow, and water vapor transport, because water is assumed to be introduced into the first layer of the GDL as a source of liquid water. At the top of the GDL, as shown in Fig. 3, under the ribs, noflux boundary conditions are used for air flow, water flow, and vapor transport; under the channel, capillary pressure is set to zero, and a vapor concentration is specified corresponding to a relative humidity (RH) value in the GC. Finally, periodic boundary conditions are used for all side pore bodies.
A number of case studies, listed in Table II, are conducted in this work. In the case studies under totally wet GC condition (i.e. RH = 100%), the employed operating and physical parameters are listed in Table III. It is noted that in this work we assume a fully hydrophobic GDL. This is because with a contact angle smaller than 135 • the air-water configuration for a specific value of water saturation in a cubic pore body is very complex. 46 Further numerical fittings may be needed to obtain this information. However, compared to a less hydrophobic GDL, it may be expected that the assumption of a fully hydrophobic GDL does not significantly impact our understanding of several important water transport mechanisms. 27 Also, it is worth noting that mixed-wettability of the GDL considerably influences water transport and distribution in the GDL. 21,24,47 But, this feature is not implemented in the present pore-network model. In the case studies of unsaturated GC, we vary the RH value in the GC to investigate its effect on the water transport in the GDL. Finally, a case study with a 3D pore network of the GDL is conducted.

Results and Discussion
Initially, the pore network is saturated with air. Then, air is slowly replaced by liquid water that enters from the bottom of the GDL. This process is called primary drainage in porous media research. Under realistic PEFC operating conditions (e.g. with a current density of 2 A/cm 2 ), it can take a minute to obtain the breakthrough of liquid water at the GDL-GC interface. But, our preliminary simulations show that the time step determined by Eq. A17 is mostly on the order of 10 −7 s. This means that the computational effort will be very heavy for a complete drainage process, even for the 2D pore network. Our preliminary simulations also show that the time step is independent of the water injection rate. Thus, one possible way to reduce the computational effort would be to speed up the drainage process by increasing water injection rate, if this does not significantly influence our understanding of water transport mechanisms in the GDL.
To verify this conjecture, four case studies were conducted for different water injection rates, corresponding to assumed current densities of 20, 200, 500, and 1000 A/cm 2 . Fig. 4 displays the distributions of liquid water in the pore network at the end of the simulations. Note that each simulation has its own flow time to ensure that at the end the same amount of liquid water is injected into the network. It is seen that increasing the water injection rate (up to a hundred times larger than a current density of 2 A/cm 2 ) has negligible effect on the final distribution pattern of liquid water. But, in the case study with a current density of 1000 A/cm 2 , an extra breakthrough location of liquid water is observed marked by a dashed circle in Fig. 4. This indicates that, at such a high flow rate, the water dynamics starts to be strong enough to slightly influence the water pathways. Fig. 5 displays the temporal evolutions of global capillary pressure and water saturation for the above four case studies. Note that the global (i.e. networkscale) capillary pressure is defined as the average of interfacial-areaweighted local capillary pressures over the whole pore network. 37 For    comparison, in each case the physical time is normalized by its own flow time (given in Fig. 4). First, no significant difference in the evolution curve of capillary pressure can be found. This may indicate that for the four different injection rates, the water flow patterns and even local air-water configurations are more or less the same with respect to the same water saturation of the network. This is further confirmed by the temporal evolutions of water distribution in the whole network (cf. supplementary materials S1, S2, S3, and S4). Second, the frequent fluctuation of the value of capillary pressure indicates that at the pore scale no real steady-state liquid water transport exists in the GDL. Third, it is seen that after breakthrough time, the global water saturation stays at a constant value of around 0.55 in each case. This means that the formed water pathways have a high capacity of removing liquid water out of the GDL, even for an assumed current density of 500 A/cm 2 . It seems to be against the experimental observation that increasing current density usually deteriorates the water flooding in the GDL. But, it should be noted that, in this work, we assume that the GC is free of liquid water. In reality, the GC flooding always limits the transport of liquid water in the GDL. In other words, increasing current density gives rise to more water coverage of the GDL surface, which in turn increases the flooding level inside the GDL. 28 Therefore, to investigate the effect of current density (i.e. water intake rate) on water flooding in the GDL, a GC flooding model must be coupled with the GDL modeling. To sum up, the above case studies confirm that we can reduce the computational effort considerably by moderately increasing the water injection rate. This does not significantly influence our understanding of liquid water transport in the GDL under the assumption that the GC is free of liquid water. For instance, with our dynamic pore-network model, it takes around 24 hours to run the case study with a assumed current density of 200 A/cm 2 on a single core of a workstation (Intel Xeon CPU E5-2670 2.6GHz), while it takes around 120 hours for the case study with a current density of 20 A/cm 2 . In the following case studies, the water injection rate corresponding to a current density of 200 A/cm 2 will be used by default, unless otherwise stated.
Usually, in the modeling of a PEFC, the two-phase Darcy's equation is used for the GDL. 2,36 But, a few studies 19,35,44 have pointed out that the Darcy-based continuum model is not applicable to the GDL, because the GDL thickness is comparable to the size of a representative volume element (REV). 45 Moreover, our dynamic pore-network modeling provides another evidence of the failure of the Darcy's equation in the GDL modeling. The supplementary material S2 shows that all invaded pore bodies experience a cyclic process of local drainage and imbibition, along with the movement of water invading front in the network. This implies that, under a small water injection rate (i.e. small capillary number value), Haines jumps occurring in the water invading front can influence liquid water transport over the whole network. In addition, it is found that (cf. supplementary materials S1, S2, and S3) increasing the water injection rate would increase the frequency of the above-mentioned cyclic processes, rather than change the water saturation and transport pathways. Obviously, this phenomenon cannot be modeled with the Darcy's equation. Fig. 6 shows the temporal evolutions of global capillary pressure, water pressure, and air pressure in case 2. It is seen that the air pressure and the difference between water pressure and capillary pressure are negligible in the drainage process. The global dynamics of liquid water transport in the GDL is very weak. It indicates that, with proper boundary conditions, a quasi-static pore-network model may be used to fast provide the information of liquid water transport pathways in the GDL. 26 This is further confirmed by Fig. 7, which shows the water pathways obtained by the quasi-static porenetwork model developed by Lee et al. 25,26 These water pathways are Figure 7. Pathways of liquid water transport in the GDL, simulated by the quasi-static pore-network model developed by Lee et al. 26 The invaded pore bodies and pore throats are marked by the black color; and the information on water saturation is not given.  identical to those obtained by the present dynamic pore-network model (see Fig. 4). Fig. 8 shows the distributions of water saturation in two other realizations (R2 and R3) at 0.6 s. Together with the water distribution in case 2, it is seen that breakthrough locations of liquid water are mainly determined by the pore structure of the GDL. No preferential locations are observed under totally wet GC condition. Fig. 9 shows the cross-section-averaged water saturation distributions along the through-plane direction in the three realizations. The highest water saturation is found at the MPL-GDL interface. This is due to the assumption that the inlet water vapor fully condenses into liquid water at the inlet. Turhan et al. 13 experimentally showed a highly nonlinear through-plane distribution of liquid water with a peak near the centre of the GDL. There would be two possible explanations to this kind of water distribution. First, it may be due to the fact that only a small number of water invading sites are formed at the MPL-GDL interface. 18,19 Second, it may indicate that non-equilibrium phase change between water vapor and liquid water exists in the GDL as an important source term for liquid water flow. 13 But, to capture this phenomenon, an advanced condensation model needs to be included in our dynamic pore-network model which will be further studied. Up to now, there has been much discussion on the flow pattern of liquid water in the GDL. 7,15,24 Based on the values of capillary number (10 −8 ) and viscosity ratio (17.5), the water flow falls in the regime of capillary fingering. 48 On the other hand, channeling flow (i.e. individual water clusters) was also observed in some experiments. 7 According to our dynamic pore-network modeling (cf. supplementary materials S2, S5, and S6), the flow pattern of liquid water in the pore network is found to be developed in the following way. At the beginning, several individual channeling clusters are formed at the bottom of the GDL, due to the uniform water input there. Then, these clusters expand and merge under capillary action while moving toward the top of the GDL. Therefore, we conclude that channeling flow and capillary fingering together govern the flow pattern of liquid water in the GDL. 7 In the operation of a PEFC, a temperature gradient in the GDL and an unsaturated GC condition are usually encountered. 6,36 It is widely recognized that water vapor diffusion assists with the water removal in the GDL. It is coupled with liquid water transport by the phase change. In this part, we mainly focus on the study of the effect of an unsaturated GC on the liquid water distribution in the GDL. To this aim, several assumptions and numerical treatments are used in the model setup. First, to reduce the computational effort for liquid water transport, the water injection rate corresponding to a current density of 200 A/cm 2 is employed. On the other hand, the diffusivity of water vapor is increased by 100 times to in turn mimic the operating condition of a current density of 2 A/cm 2 . Second, the last term in Eq. 7 is dropped in order to eliminate the unphysical advection of air phase due to the increased phase change rate (100 times). Third, we assume that the GDL is under isothermal conditions. To make a qualitative comparison with the experimental work of Boillat et al., 8 we use the following operating conditions and physical parameters. Two case studies with the GC RH of 90% and 75% are conducted. The cell operating temperature is set to 333 K. The phase change rate is set to 100 m/s to ensure an equilibrium condition between water vapor and liquid water in the GDL. Fig. 10 shows the water saturation distributions under the two different GC conditions. The flow time is 2 s to ensure that the water flow pattern is fully developed. When the GC RH is 90%, the left breakthrough location of liquid water vanishes compared to the water distribution in case 2 (see the second graph in Fig. 4). Also, a few more pore bodies under the channel are absent with liquid water. When the GC RH is reduced to 75%, liquid water mostly resides under the ribs while no liquid water is seen under the channel. A strong water separation in the GDL is formed between the ribs and the channel. Our numerical findings are qualitatively in agreement with the in situ observations of liquid water distribution using high resolution neutron radiography (see Fig. 4 in Ref. 8). Fig. 11 shows the corresponding distributions of cross-section-averaged water saturation along the in-plane direction. It is seen that decreasing the GC RH reduces the water flooding under the channel drastically, because most water can be removed out of the GDL in the vapor form. But, there is much less impact of the GC RH on the water flooding under the ribs because of the transport limitation of water vapor. This can be further confirmed by Fig. 12 which shows the water vapor distributions in the network. Large concentration gradients exist under the channel, while the area under the ribs is nearly saturated with water vapor.
Finally, the results of the case study with a 3D pore-network structure are presented. The 3D pore network has 9600 pore bodies and 24800 pore throats, representing a GDL with the dimensions of 250 μm, 1000 μm, and 500 μm in the X, Y, and Z directions, respectively. The same pore-size distribution as for the 2D pore network is used. The through-plane permeability of the generated 3D pore network is 1.02 × 10 −11 m 2 . The water injection rate corresponding to a current density of 1000 A/cm 2 is used for reducing the computational effort. The totally wet GC condition is considered. This 3D case study is aimed to gain a more realistic picture of liquid water transport in the GDL as well as confirm the findings obtained from the 2D porenetwork modeling. It takes around 120 hours to run this case study on a single core of the workstation. From Fig. 13, it is seen that the global water saturation reaches a steady-state value of around 0.3 12 , which is much smaller than that in the 2D pore-network modeling. Obviously, this is because in a 3D porous structure it is much easier for liquid water to invade toward the exit of the GDL. Similarly, as shown in Fig. 14, the in-plane cross-section-averaged water saturation is much lower than that in the 2D modeling, except at the inlet region. From  Fig. 15, we can observe a number of breakthrough locations in the middle and sides of the GC. Highly nonlinear local water saturation distribution indicates local drainage and imbibition everywhere in the domain. This is further confirmed by the temporal evolution of water saturation in the network (cf. supplementary material S7). In addition, from supplementary material S7, we observe periodic burst flow of several water clusters at the breakthrough locations, after the water flow pattern is formed in the GDL. This is very similar to the experimental observation in Ref. 9, which was called the eruptive transport of liquid water near the GDL-GC interface. But, it is worth noting that, in their experimental study, the growth of water droplets in the GC was present, which is expected to enhance the dynamics of water eruptive transport in the GDL and increase its periodic time. 9

Conclusions and Future Perspectives
A novel dynamic pore-network model is developed to study water transport in GDL. It also can model the phase change between liquid water and water vapor. Several important transport mechanisms of water in the GDL are numerically illustrated. First, it is found that, in invaded pore bodies, cyclic processes of local drainage and imbibition prevail in the network over the whole flow process. Therefore, no steady-state water transport at the pore scale exists in the GDL. This also indicates that the traditional Darcy-based continuum model is not applicable to the GDL. Second, our case studies show that channeling flow and capillary fingering together govern the formation of water flow pathways in the GDL, which are also observed in in situ experiments. Third, periodic eruptive water transport is captured near the GDL-GC interface. Under a dry GC condition, water vapor diffusion plays a dominant role in removing water out of the GDL. A strong water separation in the GDL is formed between the ribs and the channel.
To gain more insights into water transport in the GDL, several improvements for the present dynamic pore-network model should be pursued. First, the feature of mixed wettability of the GDL should be considered, which is usually due to an imperfect surface coating and the aging of the GDL. Second, the contact angle hysteresis in a pore body may play an important role in local water dynamics, due to the fact that cyclic processes of local drainage and imbibition prevail during liquid water transport in the GDL. Last but not least, water dynamics in the GC should be considered, in order to study its effect on water transport in the GDL. which is just the case for the GDL pore network (see Fig. 2). Combing Eq. A11 and Eq. A12, the final form of the discretized equation of water saturation is given as: ⎛ ⎝ V i t + where the superscript t + t denotes the value at current time step, the superscript t denotes the value at last time step. All other quantities are calculated from last time step.
Here, this scheme is called the semi-implicit scheme. 37 The fully implicit scheme is used for the discretization of water vapor transport Equation 9. Its final form is given as: where B i j = D wv i j A g i j /L i j The determination of the time step in a two-pressure dynamic pore-network model is very important to the numerical stability. Also it needs to comply with the physics of liquid water flow in the GDL. First, we introduce the minimum wetting phase saturation s g i,min in a pore body as follows. It is impossible to displace the wetting phase from the corners of a cubic pore completely. Thus, we assume that each pore body has a minimum wetting saturation that depends on the imposed global pressure difference p global as well as the blockage of the invading fluid (i.e. liquid water). The capillary blockage of the invading fluid in a pore body, p c i,block , is defined to be the minimum entry pressure of all pore throats that are connected to this pore body and not yet invaded. So, with the help of Eq. A3, the local minimum air saturation in a pore body is defined as: Then, the time step is determined by the time of filling of pore bodies by air or water phases. In the model, we allow for the local drainage and imbibition. So, we calculated t i for all pore bodies by: 37