Pore Network Modelling of Capillary Transport and Relative Diffusivity in Gas Diffusion Layers with Patterned Wettability

Engineering the wettability and microstructure of gas diffusion layers offers a powerful strategy to optimize water management in polymer electrolyte fuel cells. To this goal, we recently developed a radiation grafting technique to synthesize GDLs with patterned wettability. Despite the promise of this approach, current designs are empirically-driven which hampers the development of advanced wettability patterns. To design materials with improved transport characteristics over a range of water saturations, physically representative models can be employed for the bottom-up design of gas diffusion layers with local variations in hydrophilicity. In this paper, pore network models using topology and size information extracted from high resolution tomographic images of three common gas diffusion layer materials are presented with patterned wettability. We study the in ﬂ uence of the substrate microstructure, the hydrophobic coating load, and the hydrophilic pattern width. It is shown that tuning the wettability leads to enhanced phase separation and increased diffusive transport which is attributed to decreased gas phase tortuosity. The network model elaborates on previous experimental studies, shedding light on the effectiveness of the radiation pattern transference and the importance of matching the masking pattern with the substrate microstructure. The work opens the door for exploration of advanced patterns, coupled with ﬂ ow from gas ﬂ ow ﬁ eld designs.

Despite the potential of polymer electrolyte fuel cells (PEFCs) to provide pollution-free, fast-refueling, and long-range transportation, their current elevated cost hampers widespread commercialization. 1,2 Among the several strategies explored, [3][4][5][6][7] improving the complex water management within the electrochemical cell remains a promising strategy to realizing further improvements in power density and thereby reducing stack volume and weight. In a PEFC, the water balance needs to be carefully maintained. Humidification of the proton conductive membrane and ionomer in the catalyst layer is achieved by feeding humidified gases. Simultaneously, on the cathode side, water is electrochemically produced and must be removed through the flow fields, as water accumulation within the porous layers negatively impacts the performance by reducing gas diffusivity (i.e. flooding). [8][9][10] Avoiding flooding during PEFC operation is partially achieved through control of relevant operating conditions (e.g. temperature, pressure, relative humidity) that minimize water condensation within porous transport layers. [11][12][13][14][15] In addition, highly engineered porous materials are leveraged to optimize multiphase transport. In this context, the gas diffusion layers (GDLs) are central components that need to fulfill a number of important requirements, namely transport of reactants to catalytic layers and removal of electrochemically produced water, conduction of electrons and heat, and provision of mechanical support to the membrane-electrode assembly. State-ofthe-art GDLs are bilayered materials comprising a macroporous substrate and a microporous layer (MPL). The former is composed of an arrangement of carbon fibers that are held together by applying a binder or other methods (i.e. hydroentanglement). The density of fibers, their diameter and shape, and the manufacturing procedure determines the final material microstructure. On the other hand, the application of hydrophobic coatings, including their type, load, and application process impacts the final material surface chemistry and wettability. Optimizing fuel cell performance thus requires a detailed understanding of the interaction of GDL microstructure and surface chemistry on each of these roles. Specifically, it is critical to understand the air-water capillary pressure behavior and the impact of liquid water saturation on gas-phase diffusivity. Diffusion is the main mode of transport for the reactant gasses through the GDLs, which typically operate with 20%-30% saturation, 16 though this number can increase at higher power output and inlet gas relative humidity and also varies locally for sections of the GDL underneath the current collector ribs and lands. 17,18 Minimizing diffusive transport resistance is therefore key to overcome the mass transport limitations of PEFCs-as displayed by their polarization curves at high current densities-and will enable increased specific power output.
Recently there have been numerous advancements in the design of diffusion media for PEFCs, including the development of novel microporous layers, 7,19-21 new coating formulations, 22,23 and materials with local hydrophilicity. [24][25][26][27] The latter approach is particularly appealing as it enables engineering water and gas pathways at the micrometric scale within the electrochemical cell. Forner-Cuenca et al. [28][29][30] recently developed a synthetic methodology, based on electron radiation grafting, to synthesize GDLs with patterned wettability. During the fabrication procedure, the GDL is first coated with a hydrophobic polymer and then irradiated with electrons using masks with slits, thereby activating the polymer coating in specified locations and patterns dictated by the mask. Next, the materials are immersed in a solution containing hydrophilic monomers at moderate temperature to initiate a radical polymerization reaction at the activated polymer surface. The resulting material features an arrangement of hydrophilic and hydrophobic regions corresponding to exposed and masked material. Key advantages of the method include (1) the flexibility of the pattern design (including throughplane and in-plane), 29 (2) customization of the resulting functional properties by controlling the polymer chemistry and reaction conditions, 28 (3) the anticipated durability due to the formation of covalent bonds, 28 and (4) compatibility with large-scale continuous manufacturing process (i.e. roll-to-roll). 31 The synthesized materials have been successfully employed to tackle two distinct challenges in PEFCs. First, GDLs with patterned wettability have been assembled on the cathode side and resulted in decreased mass transfer overpotential, 32 increased power density, 30 and stable pressure drop in combination with interdigitated flow fields. 33 Second, GDLs modified with super-hydrophilic patterns on the anode compartment enabled PEFCs with evaporative cooling. 34,35 In this concept, liquid water is brought through a channel in the anodic flow field where it is wicked by the hydrophilic pathways on the GDL to simultaneously humidify the membrane and cool down the cell by leveraging evaporation heat.
Although promising, the design of patterned GDLs has largely been empirically driven, which requires large experimental resources and slows down progress. Furthermore, the wealth of substrate choice and possible treatments opens multiple avenues for explorative research. Simulations can aid the design process by enabling the exploration of the effects of optimizing multiple important design parameters (e.g. microstructure, coating technique, surface chemistry). The present study aims to explain the observed experimental results and provide a framework for aiding future optimization of the procedure.
Pore network modelling (PNM) has been shown to be an effective tool for investigating multiphase transport processes in porous materials. 36 The concept is to discretize the void space into individual pores connected by throats which form the sites and bonds (or nodes and edges) of a graph representation of the pore space. Percolation theory is then applied to calculate the distribution of phases according to rules of invasion percolation, 37 which correspond to quasi-static flow in the capillary regime. 38,39 During invasion percolation, interfaces between invading and defending phases are moved between neighboring pores in discrete steps following a path of least capillary resistance that is determined by the geometry of the pore space and the local wettability. Once phase configurations are determined, further transport simulations to determine the effective diffusivity of the network at various saturations can be performed. The main advantage of the pore network approach is the simplicity of the multiphase algorithms compared with direct numerical simulation (DNS) techniques such as volume of fluid algorithms that solve the full Navier-Stokes equations 40 and Lattice Boltzmann models. 41 DNS explicitly resolves the multiphase fluid interfaces, as opposed to rule-based percolation algorithms and require computational resources many orders of magnitude greater than PNMs. Pore network models also offer a large reduction in the number of degrees of freedom when solving the transport equations without severe loss of accuracy. 42 Together, these features enable the analysis of much larger and more representative domains which is useful for the present study, as a representative domain containing multiple masked regions spans several mm-a size that would be too large for DNS to tackle all at once.
In this work, pore network simulations were performed using three-dimensional GDL microstructures obtained with X-ray computed tomography, using the OpenPNM framework. 43 A detailed model of capillarity was developed to describe the menisci movements in fibrous GDL, considering toroidal throats with sinusoidal profiles as well as non-piston-like pore filling events. 44 Using this comprehensive model of capillary transport, the effect of the material microstructure, the coating load, and the radiation pattern geometry on the water distribution were investigated. To validate the model, we compare the simulations with previous experimental results. 29 Finally, diffusivity calculations were performed to predict improvements in diffusive transport compared to the un-patterned substrates in the presence of liquid water.

Experimental Methods
Material preparation.-Three commercial GDL substrate materials, without microporous layer (MPL), were investigated. Toray TGP-H-060, Freudenberg H23 and Sigracet SGL 24AA. The use of MPLs has proven to enhance water management but are difficult to image and their study in conjunction with GDL patterning is beyond the scope of this work. For convenience, the substrates will henceforth be referred to by the manufacturer only as Toray, Freudenberg and SGL throughout the paper and in figure legends. Fluoroethylene propylene (FEP) coating was applied in-house using the dip-coating procedure, as previously reported. 25 In this work, 30 and 70% wt. of FEP was applied by submerging the GDL substrates in the coating solution of FEP (FEPD121 DuPont 55% wt. solids) and ultrapure water to different concentrations. The same thermal treatment was applied in all cases, consisting of three steps: solvent evaporation, surfactant evaporation, and polymer sintering. Further details can be found in previous work. [28][29][30] Next, the materials were exposed to a laboratory electron beam (COMET, Switzerland) with a 200 keV acceleration voltage and 50 kGy dose. 2 mm thick stainless-steel masks were employed for the 500-930 and 250-460 μm patterns and 0.5 mm thick stainless-steel mask for the 100-500 μm was used to systematically prevent radiation from penetrating the samples and create patterns of irradiated and pristine FEP. Irradiated samples were subsequently immersed into a tubular reactor with pure N-vinylformamide (NVF) at 70°C during 60 min under oxygen-free atmosphere. 28 An in-house built goniometer was previously used to determine the surface contact angles for the GDL sections with grafted hydrophilic monomers. The lowest achievable contact angle was identified to be around 20°for FEP-grafted-pNVF on ideally flat surfaces. 28 However, to achieve this inside a GDL pore requires conformal and smooth coverage of all the surrounding fibers which is unlikely in practice and a more realistic minimum is judged to be around 40°in the present study.
X-ray computed tomography.-For each GDL substrate material, nine disks of 6 mm diameter were cut out and analyzed for weight and thickness. The sample closest to the respective mean weight and thickness was selected and imaged in an uncompressed state using X-ray computed tomography (CT) with a General Electric Nanotom m. The scans were carried out using a tube voltage of 60 kV and a current of 240 μA to acquire 1900 projections with an average exposure time of 2.5 s (∼1.6 h total scan time). Radiographic projections were reconstructed using a back-projection algorithm, resulting in reconstructed images with a cubic voxel side dimension of 1.8 μm. Boundary effects from sample cutting were excluded during processing by cropping a 4.5 mm diameter cylindrical section from the center of the scanned images. A binarized representation of the solid structure and the pore space was obtained from the raw gray-scale images by a combination of segmentation and morphological operations as described previously. 45 No differentiation between fiber and binder was possible in the current setup and both were assigned to be a single solid phase. Moreover, the sub-voxel porosity of the binder could not be resolved with the present resolution.
All the substrates have fiber diameters of approximately 10 μm, however in the case of SGL and Toray a significant amount of binder material has also been added by the manufacturer as evidenced by the through-plane slices in Fig. 1. Pore size distributions (PSD) for each network are provided for comparison in Fig. 1, which were calculated by summing the total fraction of the pore space volume attributed to pores with diameters within the binned ranges. The distributions are similar for Toray and Freudenberg, but SGL has a very different PSD with noticeably larger diameters. This can be attributed to the fact that SGL has significantly higher porosity near the surfaces compared to Freudenberg and Toray yet similar fiber sizes, meaning the void spaces are generally larger. Toray has the largest percentage of pores and throats in the intermediate range and the Freudenberg pores are the smallest in average in this series. The images were obtained for the substrates with 30% wt. coating load only. In order to simulate the percolation and diffusion through materials with higher coating loads (i.e. 70% wt.) the fibers were morphologically expanded uniformly by one voxel. This resulted in a decrease in porosity of the image that agreed closely with the actual experimental values. Importantly, this also reduced the size of the pores and throats, thereby reducing the effective diffusivity of the resulting network, in agreement with the observed performance.

Modelling
Pore network creation.-Algorithm to create a 3D pore network from the X-ray CT data.-All the pore network modeling results reported in this work were obtained with the open source package OpenPNM, 43 using networks extracted from tomograms with tools in the PoreSpy package. 46 Network extraction was performed on the binarized images with a watershed segmentation algorithm known as SNOW (Sub-Network of an Over-segmented Watershed), which has been shown to work well for highly porous, anisotropic media such as GDLs. 47,48 Anisotropy is an important characteristic of GDLs, resulting in approximately double the transport capacity in-plane compared with through-plane due to the alignment of the carbon fibers. 49,50 The SNOW extraction algorithm proceeds in two main steps. Firstly, pore regions are identified using a watershed filter, providing an image with all the voxels assigned a pore index. The region image is then processed to create a network representation by identifying the key properties of the network such as pore volume, centroid coordinates, inscribed sphere diameter and throat area, perimeter, centroid and length, among others. A full list can be found elsewhere. 47 The region image from the first step is useful for visualization by back-populating it with information from the percolation simulations, such as water configurations colorized by the invasion sequence.
In order to apply boundary conditions to the network when running simulations it is useful to create fictitious boundary pores at the edges of the domain with zero volume. 51 The procedure is simple enough when dealing with traditional cubic networks with flat and well-defined faces, but with a random topology or extracted network the task is more challenging. The automatic placement of boundary pores was achieved in this study by padding the region image prior to extraction to create fictitious pore regions around the image as described by Khan et al. 42 Determination of the contact angle.-Modelling percolation in the extracted networks requires specifying an invading phase contact angle distribution throughout the network. The contact angle is applied on a pore-by-pore basis and can be changed in three ways (see Table I): (1) coating with FEP, (2) irradiating the FEP and grafting hydrophilic polymer, and (3) random variability due to surface roughness. There is an inherent variability in the wettability of the material due to the effectiveness of covering the fiber surfaces that surround each pore and this was accounted for by introducing three potential scenarios for the fibers: uncoated, partially coated, and maximally coated. The probability of picking from the scenarios was determined by the coating load (X) as displayed in Table I. From experimental studies, 70% wt. is assumed to result in the maximal coating and other loadings are normalized and scaled linearly by this factor. 28 The present study investigates a 30% wt. material load making X = 3/7 where probabilities for individual pores being coated are given in the table below and 0% wt. where the probability of having any coating is zero.
The uncoated GDL contact angle (U) is assumed to be 80°for all substrates, 52 the coated contact angle (C x ) can be changed by radiation and grafting so is a function of in-plane position. Fully exposed sections become hydrophilic after grafting with minimum achievable contact angle assumed to be 40°; while non-irradiated regions remaining hydrophobic feature a contact angle of 100°. 28 Beam diffraction effects (e.g. forward scattering and backscattering) are accounted for by linearly scaling the coated contact angle with the distance from a mask edge (at the center of each transition region) as illustrated in Fig. 2 given the following parameters: the hydrophilic and hydrophobic mask widths and a transition width over which to linearly change the contact angle due to beam broadening. The transition width over which the contact angle is linearly scaled between hydrophilic and hydrophobic extremes is simplified to be 250 μm, based on the worst-case scenario. 53 Accurate simulation of the wettability transition regions would require a more detailed model that incorporates the dependence of beam broadening on electron beam energy, substrate thickness, density, and morphology but this exercise is outside the scope of the present work. Remarkably, the radiation grafting method has been proven to render patterns with negligible scattering using higher electron energies 53,54 . As explained in the Results and Discussion section, we hypothesize that the PNM partially fails to accurately reproduce the actual porous scaffold connectivity, which results in overpredictions of the pattern segregation. The introduced transition region on the contact angle, thus, accounts for the connectivity effect and for potential scattering events.
Once the contact angle map is determined for a given pattern width and coating load, the in-plane coordinates of the pores are used to set the contact angle for each pore. Finally contact angles are randomly adjusted up to ±10°, to account for the surface roughness of the fiber, binder and coating. It should be noted that the use of the word "maximal" to describe the coating was deliberate and different from the word "fully." As mentioned previously, it is likely that all fiber surfaces are not fully coated and consequently available for grating and sections of the carbonaceous bare substrate are likely to persist throughout the bulk with the dip coating technique. The contact angle measurements performed previously were carried out on the top surface of the GDL which is likely to have larger coating coverage than a representative section of the interior. Nevertheless, we hypothesize that improved coating techniques could result in near perfect coverage of fiber surfaces and, for this purpose, we include some comparative simulations using a minimum contact angle of 20°f or the hydrophilic regions.
Governing equations.-Capillary pressure.-Percolation models rely on determining a critical capillary pressure, above which a phase may overcome capillary forces and invade a pore through a connecting throat. The critical pressure is partly determined by the shape and size of the throat. The classic capillary pressure equation defined by Washburn is derived from the Young-Laplace equation assuming straight walled where s is the surface tension, q is the contact angle between phases as measured through the liquid phase and r is the radius of the tube. In this scenario, wettability is purely characterized by the intrinsic wall contact angle, with a threshold occurring at 90°, above which the invading phase is said to be non-wetting and requires positive pressure to invade; and below which the invading phase is said to be wetting and will imbibe into the pore at negative capillary pressures. However, previous studies have shown that the classic model does not accurately capture percolation in fibrous media which do not have straight walls. Spontaneous imbibition, predicted for materials with intrinsic contact angles of less than 90°, is not observed and the magnitude of the predicted pressure is often a poor match to experimental data given realistic contact angles. 44 The difference in behavior is attributed to the curvature of the fiber surfaces, 44,51 which create an effective contact angle that is different to the intrinsic one. A better model, first proposed by Purcell, 55 is to consider the solid surface surrounding each throat to be toroidal with a circular profile, thus accounting for the converging-diverging nature of the fiber walls. 40,44,51,55,56 In the present study, the toroidal model is adapted to have a sinusoidal profile 29  The reason for adapting Purcell's model to have a sinusoidal profile is because it allows for spontaneous imbibition when the intrinsic contact angle is very low. In previous work concerning neutrally-wettable fibrous media, 44 the Purcell model was found to be a good fit but no spontaneous imbibition had been observed. In the case of the experimental results for the altered highly-wettable GDLs studied here, some spontaneous imbibition at negative capillary pressures 29 is, in fact, observed when hydrophilic membranes are used to control entry of water at the bottom face of the GDL. For throats with solid profiles that are circular or elliptical, a ranges between 2 / p  if the meniscus is allowed to transition through the entire length of the throat without interaction with neighboring menisci and other solid objects. The implication is that even when the intrinsic contact angle (q) is close to zero, the capillary pressure required to invade the throat is greater than zero because  is always less than 2.
/ p With a sinusoidal profile, the slope of the wall features a maximum and, therefore, breakthrough pressure can be negative if the intrinsic contact angle is low enough. Setting the solid profile scaling factors a and b to both equal the fiber radius (e.g. 5 μm for most GDL materials), the sinusoidal model predicts that away from the apex, a increases to a maximum 4 max a p = and so when 2, and P 0 c < i.e. spontaneous imbibition may occur for highly wetting fluids as one would expect. In other words, the intrinsic wetting can overcome the geometrical barrier. The exact value of the maximum filling angle in the more general case will depend upon the scaling factors but occurs at the inflection point of the profile.
Effective and relative diffusivity.-Fuel cell polarization has three contributions which are well documented, 57-59 namely activation, ohmic, and mass transfer overpotentials. The mass transport limitations mainly occur at higher operating current densities when the reactant gasses cannot be delivered to the reactions sites as a sufficient rate, so the current generation becomes diffusion-limited. Increasing the effective diffusivity by reducing the tortuosity of the gas phase is the key motivation for altering the wettability locally. In other words, confining the liquid water to specific locations such that is does not increase gas phase tortuosity results in higher mass transport rates.
Diffusion in the pore network is calculated by applying Fick's law for diffusive flux between pores i and j can as follows: Figure 3. Illustration of the meniscus penetrating a toroidal throat with a sinusoidal profile along the axis of the throat. Table I. Details of the coating and base coating contact angles. X is the normalized coating load with maximum being 70% wt. C x is the intrinsic contact angle assumed to apply as a function of in-plane position, depending on the masking pattern as displayed in Fig. 2, and U is the uncoated baseline contact angle of carbon.

FEP coating
Probability of pore being coated Probability of pore being coated at 30% wt. (X = 3/7) Pore contact angle (q) Fully Coated where C is the molar concentration in pores i and j, and g D ij , is the diffusion conductance between the pores, to be elaborated on below. Performing a mass balance around a single pore i with N neighbors, assuming steady-state and no reaction, gives: Diffusive conductance between pores is defined as follows: where f s is a scale factor depending on the saturation or phase occupancy inside each pore and throat (1 × 10 −6 for air when liquidfilled). This factor is a feature of multiphase pore network modelling and is essentially used to limit the conductance of throats only to the occupying phase. The value is chosen to be very small rather than zero for numerical convergence purposes. f s could be made into a variable and used to model corner and film flow if the residual saturation is well connected but we do not consider these transport mechanisms in the present study. N r is the molecular density of the gas, i f and l i are the cross-sectional area and conducting length of the pore or throat, respectively and D ab is the diffusion coefficient of gas species a diffusing through stagnant species b (∼2 × 10 −5 [m 2 s −1 ] for oxygen in air at STP). Combining the conductance for the three elements of a conduit consisting of half-pore i, throat t, and half-pore j is done by assuming they act as resistors in series: where e is the porosity and t is the tortuosity or the network. Finally, the relative diffusion coefficient is the ratio of the effective diffusion coefficient in a partially saturated and dry state:

Results and Discussion
Capillary pressure and saturation distributions.-Capillary invasion simulations were performed for each of the GDL materials by injecting water from the bottom face, akin to water being produced at the catalyst layer and invading into the GDL to reach the flow field channels. Results of these simulations are presented in terms of "water thickness," which is comparable to neutron imaging studies, and a saturation vscapillary pressure relation. The images were produced by mapping the pore network invasion patterns back onto the regions of the tomography image as determined by the network extraction, then colorizing according to the total number of water-filled voxels in each column. The water thickness images provide a useful visual indication of how effectively the contact angle treatment separates the phases. However, the saturation vscapillary pressure curves provide the strongest validation and these are also included with comparison to experimental data from. 29 Influence of substrate microstructure.- Figure 4 provides a comparison of the different microstructures for a specific radiation mask (500-930) and coating load (30% wt.). The SGL material has larger pores compared to Toray and Freudenberg and a wider PSD and these larger pores are randomly distributed in the in-plane direction ( Figure 1j). As capillary pressure is both dependent on the contact angle and the pore size, water imbibition in to the SGL substrate results in a more random saturation pattern. Freudenberg ( Figure 1l) and Toray (Figure 1k) have tighter PSDs and smaller pores on average compared with SGL and so the contact angle changes have a much more pronounced effect on the separation of phases.
The phase separation for a certain capillary pressure can be seen by the two distinct saturation curves for the different regions of the modelling domain. The hydrophilic section is labelled by those pores having coordinates under the section of the mask exposed to radiation and hydrophobic is everything else. Where phase separation is better, the two main vertical parts of the saturation curve happen at more distinct capillary pressures. This was clearly evident in both the experiment and simulation for Freudenberg and Toray, but experimental results for SGL provided no such distinction and so a domain average was previously provided. 29 In general, the agreement between experiment and simulation is good, providing validation of the simulation parameters and algorithms. However, for the SGL case the simulation approach partially fails to predict the experimental data, in which no water segmentation was observed. We hypothesize that the microstructure of the SGL paper, specifically its surface morphology, can explain these differences between model and experiment. First, the carbon fibers that compose the SGL mat feature a large fraction of a graphitic-like binder on the surfaces. 60,61-63 The binder size features (i.e. submicrometer) lie bellow the effective resolution of the imaging technique, which causes inaccuracies on the modeling domain. Second, wetting phenomena on rough, heterogeneous surfaces is complex and local Cassie-Baxter effects might break the smooth progression of menisci that percolate through the material. The presented model is likely valid when the solid surfaces are smooth. Further research is required to adapt PNMs to accurately simulate substrates with rough fiber materials and bimodal pore size distributions.
A further discrepancy between the simulation and experimental data can be seen at low saturation for the Toray sample. As GDL materials are quite thin and have relatively large pores compared to their thickness then surface and edge effects can have a significant impact on the overall average saturation. It is possible that during the experiment, there was a preferential accumulation of water at the interface between the GDL and the hydrophilic membrane.
Influence of pattern width.-To compare the effect of radiation mask widths on the saturation, the Toray substrate with 70% wt. coating load is chosen for simulation as the most experimental data has been previously collected for this material. The match for each radiation pattern is generally very good using the base contact angles of 100°and 40°for the different regions as shown in Fig. 5. There is some disagreement for the smallest pattern width (100-500) where the simulation data shows better separation than the experimental, which could feature a limitation of the PNM approach to capture the network connectivity throghout the entire GDL thickness. However, as the 100-500 case has already proven to be ineffective experimentally for these particular GDL microstructures, 30 it does not warrant further study and verification. The variation between the capillary pressure curves clearly shows better separation with wider mask regions, highlights the need to match the radiation mask to the substrate as the average pore size determines the effectiveness of pattern transfer-larger pores are more likely to straddle the regions of the mask when the slits are of comparable width to the pore diameters.
Complete fiber surface coverage.-For all previous simulations, we assumed that the carbon fiber surfaces are only partially covered with FEP, resulting in a minimum contact angle in the exposed regions of 40°. If it was possible to obtain a conformal coating (e.g. through more advanced coating techniques), then the hydrophilic contact angle might be closer to 20°, based on previous measurements. 28 As an example, results using 20°as the base contact angle for the hydrophilic regions whilst keeping the hydrophobic contact angle 100°are presented in Fig. 6. The percolation simulation using the sinusoidal capillary model results in spontaneous imbibition for the 70% wt. coating load when applied to a Toray substrate using a radiation mask with widths of 500-930 μm. The trend is also repeated for other masks and substrates with the higher coating load. The implications for fuel cell operation could be that water is wicked from the membrane and transported to the gas channels at low capillary pressure. Such behavior is desirable for the evaporative cooling application, 34,35 in which liquid water brough in through the anodic flow field channels needs to be wicked into the GDL to evaporate (i.e. cell cooling), and to humidify the ion exchange membrane.
Relative diffusivity.-After validating the capillary invasion patterns, we then simulated gas diffusion rates through the network as a function of water saturation. Figure 7 shows the relative diffusivity data, Eq. 9, for each substrate and each coating load (with the exception of SGL where only 30% wt. coating load is shown). The SGL case with 70% wt. FEP coating load featured poor performance in experiments and was not included in previous results. 29 Both Toray and Freudenberg show marked improvements in relative diffusivity for patterned cases compared with their unpatterned simulations.
The performance improvement of Freudenberg is slightly less than Toray for increased relative diffusivity with patterning. However, the baseline for Freudenberg is better than Toray indicating that phase-tortuosity is already slightly superior and so the patterning has less room to improve. We hypothesize that this could be related to the method of fabrication of the Freudenberg which uses hydroentanglement to create a partially aligned microstructures.
The best-case scenario for relative diffusivity is to decrease linearly with increasing water saturation, as this means a reduction in total flux due to the presence of water blockages (i.e. porosity reduction) but a negligible increase in tortuosity. The former is unavoidable, but the latter can be controlled by the wettability pattern. A power-law envelope is plotted with dashed lines representing the worst and best-case scenarios for each substrate. Remarkably, the relative saturation of the Toray and Freudenberg materials approach the linear power law values for saturations around 0.3 with the higher coating loads. Figure 8 shows a comparison of all radiation patterns for the Toray substrate with both coating loads. Figures 8a and 8b are shown in terms of the absolute diffusivity value to highlight a problem that is hidden when comparing relative diffusivity curves, allowing for a more representative comparison between different substrates. The relative diffusivity for all Toray cases is also shown in Fig. 8c but for clarity we show the increase in relative diffusivity above the base case with no radiation patterning. Simply looking at relative diffusivity, one would conclude that the most effective case is Toray with 70% wt. and 500-930 patterning. At 0.3 saturation, there is an increase of 35% of the absolute diffusivity with no saturation, D .  between Figs. 8a and 8b, it can be seen that the higher coating load leads to worse absolute diffusive transport because more solid material is present in the substrate leading to reduced porosity and average pore sizes.
For all radiation mask widths, some increase in diffusive transport is seen compared with the baseline for both Toray and Freudenberg. This trend continues as saturation increases up to around 0.7, beyond which it performs essentially the same as the base case. However, this high level of saturation is above that expected in an operating fuel cell GDL with the diffusivity increase coinciding with the most relevant saturation level to fuel cell operation. The benefit from patterning the material is greatest for the widest patterns and decreases with decreasing widths as it approaches the average pore size. For all substrates, the coating load and therefore the surface coverage plays a big factor in determining the phase separation, which motivates research into conformal, thin coatings to keep overall porosity high.
In-operando studies 30 of overall fuel cell performance, comparing the 500 μm pattern with the 250 μm pattern and a baseline (or no pattern), have revealed similar trends. The findings of the present study suggest that phase tortuosity is a key target for optimization. However, the effect of the surrounding fuel cell components, such as the microporous layer and gas flow channels are not considered here and are anticipated to influence the saturation distribution in the GDL. The MPL serves to limit the number of liquid injection points from the catalyst layer into the GDL 64-69 and the gas channels serve to remove liquid water from the GDL via convection and evaporation. The ribs of the current collectors also create cooler sites where liquid may condense if the relative humidity in the gas phase is high enough. Very hydrophilic regions will serve as spongesand could be used to direct water from the MPL/catalyst regions and also the condensation sites into the center of the gas channels where it can be removed; therefore the alignment of the pattern with the gas channel design is likely to be an important design consideration. Simple parallel channels are facile platforms to tailor GDL wettability patterns, though more geometrically complex designs 68 may be more challenging to accommodate but potentially more rewarding. Future modelling efforts should consider these advanced engineering parameters.

Conclusions
Optimizing complex multiphase transport through porous diffusion layers in PEFCs can enable further improvements in power density and durability. To this end, diffusion media with patterned  wettability have emerged as a promising solution. And, while current designs are informed by experimental studies, simulations can be used to accelerate the discovery of new materials. To this end, we investigated the influence of several material parameters (i.e. substrate microstructure, coating load, and pattern width) using a microstructure-informed pore network model with distributed wettability to compute capillary pressure characteristic and diffusive transport. The model accurately matches capillary pressure experiments over most experimental parameters, but partially fails to describe imbibition for the narrowest hydrophilic patterns. We hypothesize that connectivity of the real structure might not be fully captured by the pore network topology.
The microstructure, coating load, and radiation pattern width all have a significant effect on the capillary transport and consequently the diffusive transport. We find that the substrate with the broader pore size distribution (i.e. SGL) has large pores (>200 μm) that tend to fill at lower pressures irrespective of the applied contact angle, resulting in the poorest separation of the phases and higher gas phase tortuosity. Substrates with tighter pore size distributions and average smaller pores (i.e. Toray and Freudenberg) enable more effective transfer of the wettability pattern. Of the materials investigated, modified Toray paper resulted in the greatest relative improvement in diffusive transport when partially saturated. Interestingly, we found that fully-coated Toray with 500 μm width hydrophilic domains features nearly optimal diffusivity-saturation characteristic at saturations of ∼30%, which is representative of PEFCs operation at higher current densities. This finding motivates research into conformal, thin hydrophobic coatings for gas diffusion electrodes. Now that a robust modelling strategy has been developed, it will be possible to carry out a topological optimization of the pattern geometry, including wettability modifications in the through-plane direction. Future theoretical work should include other relevant physical phenomena (e.g. phase changes, electrochemical reactions) and relevant surrounding components (e.g. flow fields, microporous layers) to more accurately capture the device performance.
[ · ] a = ¢ From geometrical consideration: The radius of curvature of meniscus inside the throat: