key: cord-1026321-uepkkmnv
authors: Farzaneh, Meisam; Ström, Henrik; Zanini, Filippo; Carmignato, Simone; Sasic, Srdjan; Maggiolo, Dario
title: Pore-Scale Transport and Two-Phase Fluid Structures in Fibrous Porous Layers: Application to Fuel Cells and Beyond
date: 2020-11-23
journal: Transp Porous Media
DOI: 10.1007/s11242-020-01509-7
sha: a56edd804e133d8cf6486ea687299d81636f4c9d
doc_id: 1026321
cord_uid: uepkkmnv

We present pore-scale simulations of two-phase flows in a reconstructed fibrous porous layer. The three-dimensional microstructure of the material, a fuel cell gas diffusion layer, is acquired via X-ray computed tomography and used as input for lattice Boltzmann simulations. We perform a quantitative analysis of the multiphase pore-scale dynamics, and we identify the dominant fluid structures governing mass transport. The results show the existence of three different regimes of transport: a fast inertial dynamics at short times, characterised by a compact uniform front, a viscous-capillary regime at intermediate times, where liquid is transported along a gradually increasing number of preferential flow paths of the size of one–two pores, and a third regime at longer times, where liquid, after having reached the outlet, is exclusively flowing along such flow paths and the two-phase fluid structures are stabilised. We observe that the fibrous layer presents significant variations in its microscopic morphology, which have an important effect on the pore invasion dynamics, and counteract the stabilising viscous force. Liquid transport is indeed affected by the presence of microstructure-induced capillary pressures acting adversely to the flow, leading to capillary fingering transport mechanism and unstable front displacement, even in the absence of hydrophobic treatments of the porous material. We propose a macroscopic model based on an effective contact angle that mimics the effects of the such a dynamic capillary pressure. Finally, we underline the significance of the results for the optimal design of face masks in an effort to mitigate the current COVID-19 pandemic.

Hydrogen is a clean and flexible energy carrier that plays an important role in the decarbonisation of the industry and transport sectors (European Commission 2020). As an example, hydrogen can be used in fuel cells to generate electricity, power or heat. Owing to their many advantages such as high power density, near-zero emissions and low noise, hydrogen-based Polymer Electrolyte Membrane fuel cells (PEMFC) are a reliable alternative to fossil-based fuel engines for transport, industrial and household applications (Cano et al. 2018) .

The core of a PEMFC consists of the polymer electrolyte membrane sandwiched between two thin catalyst layers and two porous electrodes, which are called gas diffusion layers (GDLs). The catalyst layer has a great impact on the cost of fuel cells, and research efforts have been devoted to its improvement (Brown et al. 2019 ). On the cathode catalyst layer surfaces, protons and oxygen participate in an exothermic chemical reaction the products of which are water and heat. Within the pores and interstices of the GDL porous microstructure, the reactant gas and the produced water flow to and from the catalyst layer. The GDL is the key component for regulating the transport phenomena inside the cell and, in turn, its design is an important factor to determine cell performance (Majlan et al. 2018; Balakrishnan et al. 2020) . In particular, the GDL porous microstructure is fundamental in determining the two-phase (liquid-gas) dynamics and distribution in the cell, a complex mechanism of mass transport that affects the management of the produced water (Berg et al. 2004) .

Managing the water balance inside fuel cells is a great challenge. An optimal operation of a fuel cell requires an adequate hydration of the membrane for allowing proton conductivity. On the other side, excessive water fills up the pores and interstices of the GDL, blocking the access of reactant gas to the surfaces of the catalyst layer and decreasing the available active sites for reaction (Culubret et al. 2020) . Consequently, the accumulation of produced water, which is mostly observed at high current densities, limits cell performance (Li et al. 2008) . The imbalance of water inside the cell is particularly important in new types of fuel cells, such as anion exchange membrane fuel cells (Eriksson et al. 2019; Chen 2020 ). An inefficient water management in fuel cells is also recognized as one of the leading mechanisms that induce fuel starvation and cell degradation-two well-known ongoing issues that reduce the lifetime of fuel-cell systems Nandjou et al. 2016) . It is therefore of great interest to investigate the interplay between two-phase flow phenomena and GDL microstructure, in order to establish efficient water management strategies, and increase the performance, robustness and competitiveness of fuel cells in the energy market.

GDLs are composed of a highly porous (porosity > 70% ) woven carbon cloth or carbon paper with 100-500 μm thickness. The diameter of fibres composing the GDLs is typically 10-15 μm , and the mean pore size is 50-120 μm (Gostick et al. 2007; Lu et al. 2010) . The GDL porous microstructure provides a complex network of narrow and dilated pore throats giving rise to local constrictions of the flow paths. This complex topology induces high curvatures of the liquid-gas interface and, given the presence of surface tension, it produces significant capillary forces that act as a resistance to the flow (see e.g. Lu et al. (2019) ). The presence of these capillary forces is a source of pressure build-up and intermittent water dynamics in the GDL. This behaviour often results in impulsive forces that flush out the water to restore the cell operational performance (Li et al. 2008 ).

The strongly intermittent transport phenomena in GDL samples have been experimentally observed. As an example, fluctuations of the pressure drop, induced by the breakdown and re-build of water flow paths, have been measured (Mortazavi and Tajiri 2014; Niknam et al. 2020 ). Water can emerge from different preferential locations on the porous surface when transported through a GDL, an observation that suggests that the spatial distribution of flow paths is also intrinsically dynamic (Lu et al. 2010) . The observed behaviours of liquids penetrating into GDLs are consistent with the physical description of two-phase flows in porous media, as characterised by intermittent fluid connectivities that are strongly related to the porous microstructure (Reynolds et al. 2017) . Recent advances in imaging analysis techniques have allowed the access to pore-scale flow visualisation within fuel cells, for instance, by means of nuclear magnetic resonance (NMR) imaging, neutron imaging and X-ray computed tomography (CT). Even though these experimental methods provide valuable information on microscale water transport behaviour within fuel cells, their application is often challenging due to opacity of the system, limited spatial and temporal resolutions, cost and availability (Bazylak 2009) .

Numerical simulations provide a valuable tool for overcoming the limitations related to experimental techniques, in terms of both investment and temporal costs (Andersson et al. 2016) . As an example, the two-phase flow in fuel cells has been mathematically modelled and simulated using continuum approaches based on the volume averaging of physical quantities (Nam and Kaviany 2003) . However, GDLs suffer from the lack of length scale separation since the thickness of a GDL measures only few pore sizes, making the application of continuum models controversial (Rebai and Prat 2009) . To incorporate the GDL pore-scale structure in numerical simulations, pore network (PN) models (Lenormand et al. 1988; Rebai and Prat 2009) and artificially generated fibrous material have been used Deng et al. 2019; Hao and Cheng 2010) . PN models represent the microstructure of porous media as a two-or three-dimensional network of idealized pores and throats, whereas artificial materials are usually created by means of stochastic techniques based on macroscopic design specifications, such as porosity and fibre size. Such studies can catch the small-scale transport mechanisms usually filtered out in macroscopic continuum models, but they rely on simplified microstructures. Therefore, pore-scale simulation of liquid water transport and distribution inside realistic GDLs geometries can unveil a much richer physical picture. An interesting approach for the investigation of pore-scale two-phase dynamics makes use of X-ray computed tomography for reconstructing porous media as input for microscale simulations (see e.g. Raeini et al. 2014) . Some studies have exploited this strategy focusing on reconstructed GDL geometries (Rama et al. 2012; Gao et al. 2013) . These studies investigated the three-dimensional evolution of liquid invasion in hydrophilic and hydrophobic GDLs under different pressure gradients. The prominent effect of microstructure on two-phase fluid distribution and dynamics has been observed in samples with surface areas of the size of only few pores, i.e. of the order of a hundred micrometres squared. These insights depict the significance of such analyses and motivate further investigations of two-phase flows in GDLs, applying the analysis to samples of larger size. Indeed, a comprehensive analysis of the pore-scale two-phase transport in GDL samples is still lacking, thus limiting the possibility of understanding the microscopic behaviours of the fluids inside the fuel cell.

Recent studies have shown the potential of combining lattice Boltzmann simulations and X-ray computed tomography for the investigation of single-phase flows in large samples of porous battery electrodes, notably providing the computation of electrode conversion and reactive properties (Kok et al. 2019; Maggiolo et al. 2020 ). In the present study, we adopt a similar approach to perform pore-scale two-phase lattice Boltzmann simulations in the porous structure of a commercially available carbon paper GDL for fuel cells (AvCarb MGL370), reconstructed via X-ray computed tomography. Such fibrous layer materials are characterised by much larger lateral dimensions compared to the thickness. The porescale heterogeneities are expected to be statistically significant along the minor dimension, which is of the size of a few pores, and to play an important role in the global performance of the material. We are particularly interested in the effect that the microstructural variations have in determining the dynamics along the medium thickness.

Three samples with a surface area and thickness of the order of 1 square millimetre and 200 μm, respectively, are extracted from the reconstructed materials and used as input for three-dimensional two-phase flow simulations to mimic a liquid-gas (water-vapour) system. The results of the simulations are presented in terms of liquid water dynamics, spatial distribution and flow preferential paths. We capture the two-phase flow dynamics within the GDL at the pore-scale, and we illustrate that the flow experiences a complex dynamics in the form of different spatial patterns, where inertial, capillary and viscous forces all contribute to the intermittent flow behaviour. These results provide important information for expanding the knowledge about the physical mechanisms underlying the intermittent mass transport phenomena in fuel cells, which often cause fluctuation in the performance and unstable operation of a fuel cell (Ji and Wei 2009; St-Pierre et al. 2000) . Finally, we also propose a model to describe the water transport in the GDL based on extracted effective permeability values.

The results on the two-phase flow dynamics and fluid structures, and their relation to the porous microstructure, provide new insight on the pore-scale mechanisms of transport in fibrous materials, also beyond the application to fuel cells. In the final part of the paper, we will share important aspects regarding the optimal design of face masks with an improved filtering performance, and we discuss how the use of certain type of materials in the production of personal protective equipment can possibly contribute to limit the transmission of viruses and mitigate the current COVID-19 and future pandemics.

The manuscript is organized as follows: in Sect. 2, the material reconstruction technique and the numerical methodology for simulating two-phase flows are presented; in Sect. 3, we present and discuss results of the numerical analysis; in Sect. 4, we discuss the relevance of the results for the design of face masks devoted to the mitigation of virus transmission, and in Sect. 5, we summarise the main findings of the research.

The three-dimensional geometry of the GDL material (AvCarb MGL370) is acquired and reconstructed via X-ray CT. For that purpose, a metrological CT system (Nikon Metrology MCT225) is used, which is characterised by micro-focus X-ray source (minimum focal spot size equal to 3 μm ), 16 bit detector with 2000 × 2000 pixels and cabinet ensuring controlled temperature of 20 • C. The used CT scanning parameters are reported in Table 1 . In order to effectively acquire the geometry of microscale fibres, the CT metrological structures resolution has to be maximised (Zanini and Carmignato 2017) . To this end, the voxel size is minimised by reducing the relative distance between the sample and the X-ray source. This is possible because small portions are cut out from the initial material, with dimensions of about 3.2 mm × 3.2 mm × 0.37 mm . Since the thickness (0.37 mm) is lower than the other dimensions, eight samples are packed and scanned with a single acquisition procedure. As a consequence, the scanning time reported in Table 1 is the time needed to scan up to 8 material portions. In addition, to further improve the scanning resolution, the focal spot size is also minimised by keeping the X-ray power below 7 Watts. A large number of bi-dimensional projections (i.e. grey-scale images representing the local attenuation of X-rays occurring when traversing the samples material) are acquired at a different angular position of the samples during their revolution around the rotary table axis. Starting from these projections, a three-dimensional reconstruction is conducted using a filtered back-projection algorithm (Feldkamp et al. 1984 ) and cross-sectional images (as the one seen in Fig. 1 ) are extracted for each sub-sample (with the distance between consecutive sections equal to the voxel size), exported into image stacks and binarised to discriminate between material (black) and pores (white).

Even if the resolution is maximised, the smallest fibres may present discontinuous shapes (as seen in Fig. 1 ), because their microstructure cannot be fully captured by the CT scan. For this reason, a filtering procedure is operated, based on the cubic interpolation of three-dimensional geometry in a grid with increased resolution. Following this image analysis technique, we obtain an enhanced three-dimensional representation of the material, with a final resolution of 1.6 μm . Compared to the original material's specifications, we obtain a slightly higher porosity, from 78 to 80%. We also have to reduce the thickness of the sample, from the original 370 μm to 192 μm in order to avoid unphysical effects at the borders due to the presence of adjacent materials at the two sides of the GDL during the CT acquisition.

Subsequently, we select three portions (representative element volumes, REVs, Fig. 2 ) of the reconstructed GDL to serve as input for the numerical simulations. The three REVs are chosen in order to increase the statistical confidence of the simulation data. The computational domain of each sample is then characterised by a cross-sectional surface area Fig. 1 A cross-sectional view extracted from the CT-reconstructed volume and binarised to discriminate between material (black) and pores (white). The size of the total surface area is ∼ 3.2 mm × 3.2 mm of 819.2 μm × 819.2 μm and a thickness of 192 μm , to which correspond 512 × 512 and 120 computational cells, respectively. The pore size distributions of each materials samples are presented in Fig. 3 where we also depict the pore size distribution averaged over the 3 samples. We observe that the material has a microscopical structure characterized by pore sizes that range from a few microns to tens of microns. The average pore radius results r p = 40 μm , on the basis of which we argue that the size of the domain, with a surface area approximately 10 2 times larger than the average pore size, well captures the spatial variation that characterizes the GDL microstructure.

The lattice Boltzmann method is a numerical approach that allows the study of multiphase flows in porous media at a low computational cost . It is well suited for numerical solutions of complex flows with pore-scale accuracy, and it is therefore able to capture the two-phase fluid-dynamic behaviours determined by the material microstructural characteristics of a fuel cell GDL. The GDL geometry, reconstructed via CT and filtered as explained in Sect. 2.1, is here used as input for pore-scale simulations. We make use of a single-component two-phase flow model to mimic a liquid-gas (water-vapour) system by means of the numerical methodology proposed by Shan and Chen (Shan 1815) for the simulation of surface tension effects. Such a methodology does not rely on the dynamic solution of the two-phase interface. Instead, the surface tension force acting between the two phases is simulated by van der Waals-like intermolecular interactions. The two-phase flow is solved on a three-dimensional regular Cartesian grid (D3Q19) by means of the local solution of the probability distribution function f r , a statistical quantity that indicates the molecular probability of fluid mass and momentum (Succi 2001) . The lattice Boltzmann equation reads as:

where = (x, y, z) is the position vector, t is the time, r is the so-called discrete speed, f eq r is the equilibrium distribution function, and the subscript r = 1 … 19 indicates the lattice directions. The relaxation time is directly linked to the fluid dynamic viscosity as

where c s is the speed of sound and = (( , t) the fluid density. The fluid is forced by an applied body force F r that mimics the effect of a pressure gradient −∇ x P acting on the x direction. The body force is defined as (Guo et al. 2002): where w r represents the weighting coefficient and = ( , t) the Eulerian velocity vector. It should be noted that in our simulations the pressure gradient ∇ x P acts only in the streamwise direction x, whereas it is null along the transverse directions y and z (see also the sketch in Fig. 4 ). The value of the applied pressure gradient is discussed in the next subsection.

The equilibrium distribution function is defined along the r-th direction as:

(1) (Rabbani et al. 2014) . It is observed that all the three REVs have similar pore size characteristics where eq = + ( − 1∕2)F sc (Chen et al. 2014) . The interaction force between different phases F sc ( , t) is based on a density-dependent pseudopotential function ( ) (Shan 1815):

where G determine the strength of interaction between the two phases ( G = −5.5 in the present model). The form of the pseudopotential function, defined as ( ) = 1 − e − , yields a non-ideal equation of state that allows the coexistence of the two phases:

One of the known limitations of the two-phase lattice Boltzmann models is the emergence of numerical instabilities and spurious currents in cases with high density and dynamic viscosity ratios between the phases. In the present study, the dynamic viscosity ratio is set as M = ∕ g = 35 (with and g the viscosities of liquid and gas, respectively), a value that we found as a good compromise between accuracy of the solution ( M = 100 for water-air systems) and numerical stability. The reader is referred to the work presented in Pettersson et al. (2020) for the two-phase flow algorithm validation.

The fluid-wall interaction at the three-phase contact line (liquid-gas-solid), in terms of chemical properties of the surface (wetting), can be tuned by means of an averaging of the pseudopotential values in the neighbouring nodes N:

The resulting contact angle is determined by the value of the parameter w (De Maio et al. 2011; Vaitukaitis et al. 2020 ). The contact angles in fuel cells GDLs have been measured to Fig. 4 a Schematic of the computational geometry in the present work; The flow cross section A is 819.2μm × 819.2μm and the GDL thickness L pm is 192 μm b Representation of the liquid front penetrating through the space between two fibres and the induced liquid-gas interface curvature of radius r c . The value of r c can be approximated with the pore throat radius, i.e. r c ≈ r p

A x widely vary, over a range ≈ 80 • -150 • on the basis of temperature and hydrophobic agent content (Lim and Wang 2004) . The value of the contact angle is known to affect fuel cells performance (see e.g. Park and Popov (2009) ) but less is known about the effect of the porous microstructure. Both wettability and medium geometrical disorder are expected to greatly impact fluid displacement in porous media (see e.g. Wang et al. (2019) ), thus motivating the need of a deeper understanding of the effects of pore-scale microstructural traits on multiphase flows in fuel cells. In the present study, we address such a scientific question by simulating the two-phase flow dynamics within a GDL without hydrophobic treatments, and the equilibrium contact angle is set to eq = 90 • .

Finally, the macroscopic fluid density and momentum are determined through the statistical averaging of the molecular probability distribution functions as (Succi 2001 ):

We investigate the dynamic infiltration into the GDL of a volume of water V ∼ 0.13 mm 3 . Such a volume is equal to the total volume of the investigated porous samples AL pm , so that we ensure that all the media can possibly be filled with water, i.e. a situation where the cell is experiencing high flooding. A sketch of the computational domain is presented in Fig. 4a . Along the streamwise x direction, we impose a pressure gradient that drives the water through the GDL microstructure and periodic boundary conditions at the domain boundaries. At the boundaries along the transverse directions y and z, the symmetry is guaranteed by imposing free-slip boundary conditions.

When the water produced on the catalyst layer approaches a pore throat, such as the constriction induced by intersecting fibres (see e.g. Fig. 4b ), it is temporary held up and its transport impeded until the local water pressure exceeds the capillary pressure induced by the liquid-gas interface curvature. The competition between the values of the pressure experienced by the water and the pressure threshold induced by the pore geometry is the mechanism that determines if the pore throat is infiltrated or not. Based on the mean pore radius of our sample, the capillary pressure threshold can be in the first instance estimated as P t ∼ 2 ∕r p = 3.6 kPa (with = 0.072 N/m the water-air surface tension). Such a value is indeed consistent with experimental observations of water penetrating a GDL and with the corresponding measured pressure build-up (Mortazavi and Tajiri 2014) .

We then assume that the driving force for water motion is the pressure gradient −∇ x P determined by the difference between the pressure build-up at the GDL entrance (at the catalyst layer-GDL interface) and the pressure of the gas phase at the end of the porous domain (at the GDL-gas channels interface). The gas channel is assumed to be at atmospheric pressure. Following this reasoning, the pressure gradient results −∇ x P = P t ∕L pm = −187 kPa∕cm . This is the value of the pressure gradient that forces the flow in our simulations, implemented through equations 1 and 2. We neglect gravitational effects, since in fuel cell applications, the associated volumetric force g ∼ 0.1 kPa∕cm results three orders of magnitude lower. It should be noted that the value of the water pressure at the catalyst layer-GDL interface can vary greatly during cell operation. For instance,

such a variation can be due to the intermittent production of water that a cell experiences or to the complex effects that microstructural and chemical properties of the catalyst layer induce on the thermodynamical properties of water (Eikerling and Kornyshev 1998) . The presence of a microporous layer in between the catalyst layer and the GDL can complicate the picture even more. However, the water production and behaviour in the catalyst layer are out of the scope of this study. Since we are interested in following the time-dependent dynamics of water induced by the pressure build-up in the cell electrodes, we consistently apply the pressure conditions relevant to such an investigation. Under such conditions, as simulations will show, liquid penetrates in some pores of the medium and will accumulate in others, producing a rich spatial and temporal behaviour that unveils important insights about water transport in fuel cell GDLs.

We first qualitatively look at the dynamics of water infiltrating the reconstructed GDL. Figure 5 shows an example of water infiltration for one of the three investigated samples, at three different times. We observe the predicted heterogeneous behaviour, with some pores filled much faster than others, and the presence of pronounced curvatures of the liquid-gas interface induced by the geometrical constraints of the porous medium. We also observe that at time t * ∼ 0.8 , the liquid phase has reached the GDL outlet, emerging in the form of spatially distributed fingers.

Three snapshots of the displacing liquid water from lattice Boltzmann simulations in sample 3: a t * = 0.2 , b t * = 0.5 and c t * = 0.8 . The visualization illustrates that the liquid phase is squeezed into the microscopic void space of the fibrous material and transported through the porous medium microstructure. The grey-coloured structure represents the solid fibres. At time instant t * ∼ 0.5 , the breakthrough of the liquid-water front occurs and water emerges from the surface of the GDL From a macroscopic perspective, the liquid infiltration can be quantified by the percentage of volume of pores invaded by water at each instant, i.e. by the porous medium saturation S(t). In Fig. 6 , we report the variation of saturation with time. For this figure, as for the other ones throughout the paper, the computed quantities refer to the average value among the three different REVs (see also Sect. 2.1). For comparative purposes, we non-dimensionalise the time as t * = t∕t c , where t c = L pm ∕U c is the time that the considered liquid would take to fill a tube of radius r p and length L pm , when it experiences an average characteristic velocity U c = −∇ x P r 2 p ∕(4 w ) (i.e. the Poiseuille velocity averaged in time during infiltration of a liquid in a pipe at a constant pressure drop). From this first comparison, we observe that at the latter characteristic time, i.e. when t * = 1 , the porous medium is not completely filled and S(t * = 1) < 1 . At a time t * ∼ 0.5 , we observe the appearance of the first liquid cluster on the GDL outlet surface, a time that we herein refer to as the breakthrough time. At a time t * ∼ 1 , roughly 1% of the total water has emerged along certain preferential flow paths. In the next subsection, we will further discuss the mechanism that determines the spatial distribution of these preferential flow paths. Figure 6 also shows that the liquid dynamics is characterised by three distinctive regimes: (i) a ballistic one where the saturation increases as S(t) ∝ t * 2 , (ii) a less pronounced infiltration dynamics, where saturation increases with the square root of time, S ∝ t * 1∕2 , and (iii) a final stage where the saturation approaches a constant value. This behaviour points out that the initial infiltration is inertia-dominated, with the average liquid penetration height h(t) = S(t)L pm advancing with an increasing velocity dh∕dt . We note that this inertial regime occurs when the liquid invades the initial part of the GDL (approximately the 20%), and it can also be affected by the water dynamics in the catalyst layer. We then argue that in the presence of high pressure build-ups at the catalyst layer-GDL interface, this mechanism can significantly impact the water transport in fuel cell GDLs. Fig. 6 Water saturation as a function of time (left ordinate axis). The liquid front initially propagates with a rate ∝ t * 2 , due to inertial forces, whereas at the later stages, when viscous effects becomes important, the liquid penetrates ∝ t * 1∕2 . The liquid saturation at long times is fitted by Eq. 12. The figure also reports on the right ordinate axis the cumulative fraction of the liquid emerging from the GDL outlet surface, exit = V exit ∕V (dashed line), with V the total liquid volume in the numerical domain. The latter plot is showed till 2.5t * , when the water that emerges from the outlet is reconnected with the inlet liquid front due to the periodic boundary conditions When the liquid penetrates the porous medium, the amount of fibre surface wetted by the liquid increases and, in turn, the viscous resistance becomes increasingly important. At longer times, the macroscopic liquid transport is then determined by the competition between the pressure drop −∇ x P L pm and the viscous forces acting on the liquid w h . The macroscopic liquid infiltration rate does not increase, but it rather decreases proportionally to h:

Such a competition between the pressure and viscous forces leads to the scaling S ∝ t * 1∕2 , which is the classical result of Darcy's law for a liquid infiltrating a porous medium (that is, the equivalent of the Lucas-Wahsburn solution for a flow driven by a pressure drop). We observe a direct transition from the initial inertial regime S ∝ t * 2 to the viscous regime S(t) ∝ t * 1∕2 , and thus a lack of the often observed convective regime S ∝ t * . This kind of a transition is due to the lack of separation of the time scales characterising the inertial time t i = √ r 3 p ∕ and the viscous time t v = r 2 p ∕ , that is, to the relatively high value of the Ohnesorge number Oh = t i ∕t v = √ 2 ∕( r p ) ≈ 10 −1 . This observation points out that the viscous forces start to play a dominant role at relatively short times and that the convective regime, which is usually observed at t i < t < t v , occurs in a short time interval (Stange et al. 2003) . Since in such a short period preceding the viscous regime, an inertial-viscous transition regime is also expected (Fries and Dreyer 2008) , we deduce that the convective, purely inertial regime cannot be observed here. After a time t * ≥ 2 , the saturation growth rate is reduced and tends asymptotically to zero, dS∕dt → 0 . As we will see in the next section, the transport behaviour at long times is in fact characterised by the occurrence of preferential flow paths that connect the inlet and outlet of the medium, and the global mechanism of medium filling is inhibited.

The inertial and viscous forces are not the only forces contributing to water transport. As we previously argued, the microstructural traits of the GDL induce strong curvatures of the liquid-gas interface and give rise to capillary forces. We measure the time-dependent average capillary pressure P c (t) = P l (t) − P g (t) − P v (t)∕2 , as the difference between the spatially averaged liquid and gas pressures P l,g (t) = 1∕V l,g ∫ V l,g p( , t) dV (with V l and V g the liquid and gas volumes, respectively), minus half of the estimated viscous pressure drop in the liquid from Darcy's law P v (t) = u l h(t) 4∕r 2 p (with u l the liquid-phase velocity). For this estimation, we assume a linear decrease in the pressure along the flow direction because of the viscous forces in the liquid, and we also assume negligible the gas viscous forces; following these assumptions, we can write P c (t) + P v (t) = P in (t) − P out (t) , with P in (t) and P out (t) indicating the pressures at the inlet and outlet of the medium, and derive P c (t) = P in − P v − P g = P l − P v ∕2 − P g . A positive value of P c indicates that the liquid-gas interface is convex with respect to the streamwise direction x, and it represents a resistance to the flow. On the other side, a negative or a null value point to a concave or flat interface and do not provide resistance to the flow. Figure 7 presents the variation of capillary pressure and the values of the capillary number Ca(t) = dh∕dt w ∕ during the liquid infiltration. The capillary pressure is nondimensionalised with the estimated capillary pressure based on the average pore size P t .

The computed values are plotted as a function of the saturation of the medium. In the early stage, the capillary pressure increases as the liquid enters the GDL, i.e. dP c ∕dS > 0 . We can imagine this increment as a result of the initial inertial stage, where the strong

dynamics induces incrementally larger curvatures of the interface. This argument is further corroborated by the increasing value of the capillary number. After the initial growth and after having reached a maximum, the capillary pressure slightly decreases. The observed maximum at P c ∼ 1.2P t suggests the presence of pore throats smaller than the average pore size, which induce high interface curvatures and capillary pressures, when the liquid front has advanced, on average, roughly the 20% of the medium length. To the observed decrease in the capillary pressure at S(t) > 0.2 corresponds a decrease in the liquid infiltration rate, quantified by the capillary number Ca(t) , and, in turn, the beginning of the regime where viscous forces incrementally contribute to the flow resistance. The regime transition occurs at a time comparable with the nondimensional viscous time t * v ∼ 0.64 . However, the viscous time characterising such a transition can significantly vary depending on various conditions defining the two-phase system, even in simple geometries, and thus the formulation of t v chosen here offers us only an approximated indication of the expected flow behaviour (see e.g. Stange et al. (2003) for an interesting discussion on the topic). We indeed have to consider the additional effect that the microstructure could induce on the partitioning between inertial, viscous and capillary forces. We then argue that the transition between flow regimes is likely to be also induced by a significant geometrical constriction, i.e. by the presence of a significant amount of small pore throats. Given the observed variation of pore radius values in the samples (see Fig. 3 ), we expect both small and larger pore throats delineating the morphology of such a constriction. The liquid phase breaks through the largest pore (with the smallest capillary resistance), the interface curvature is reduced and in turn the capillary pressure. The infiltration velocity also decreases as a consequence of the increased viscous resistance. By looking at the observed macroscopic behaviour of the liquid infiltration shown in Fig. 6 , we immediately recognize that to S(t) ∼ 0.2 approximately corresponds the transition from (a) (b) Fig. 7 Right panel a capillary pressure and capillary number plotted as functions of water saturation. In the inertia-dominated region, S < 0.2 , inertial effects dominate, the liquid accelerates and the capillary pressure sharply increases, i.e. dCa∕dS > 0 and dP c ∕dS > 0 . For S > 0.2 , the flow dynamics is characterized by both capillary and viscous resistance forces. In this regime, we observe, dCa∕dS < 0 and dP c ∕dS ∼ 0 . Left panels: b front width h ′ as function of time, and c the non-dimensional pressure drop along the front width P f ∕P t as function of the front width h ′ the ballistic dynamics to the later regime S(t) ∝ t * 1∕2 . In the next section, we will further evaluate the impact of the GDL morphology on the dynamic transition between two-phase flow regimes. In the late stage, the capillary pressure approaches a constant value of P c ∼ P t while the capillary number tends to zero. The latter observation does not imply that the liquid flow is null, but rather that the saturation does not change in time and the liquid follows preferential paths, so that Ca(t) ∝ dh∕dt → 0 . The initial trend suggests that the capillary pressure experiences a dynamic behaviour, which results proportional to the average velocity of infiltration and dP c ∕dS ∝ dCa∕dS , as already observed in the literature (see e.g. (Weitz et al. 1987) ). It is important to stress that for S(t) > 0.2 , both the viscous and capillary forces play an important role in determining the liquid motion, so that this regime can be referred to as a capillary-viscous one.

In order to further characterise the two-phase interface displacement, we estimate the front width (or roughness) as the root mean square of the front fluctuations as h � (t) = ( ∑ A (h * (y, z, t) − h(t)) 2 ∕(A )) 1∕2 (with h * the local front position along x). We observe that in the viscous-capillary regime, the front width increases in time, as can be observed in Fig. 7b . We also observe that in such a flow region, the pressure drop along the front width P f (h � ) (between the base h − h � and tip h + h � ) decreases almost linearly with h ′ , in Fig. 7c . These two observations suggest that after the initial invasion event, the infiltration occurs in an unstable manner. A decrease in the pressure difference between base and tip of the front indeed suggests that as the front becomes rougher, the probability of invasion at the tip increases with respect to that at the base, promoting unstable front displacement. In another equivalent perspective, the viscous resistance (pressure drop) along the front width decreases during kinetic roughening. A similar destabilising mechanism has been discussed for instance by Xu et al. (1998) . In the next section, we will observe that it is the microstructure variation along the through-thickness direction that is responsible for such a peculiar front dynamics.

Pore-scale numerical simulations provide much deeper insights than do the models based on macroscopic quantities. As it can be observed from Fig. 5 , the spatial distribution of liquid in GDLs results heterogeneous and, as discussed in the previous section, the front exhibits kinetic roughening. The heterogeneous redistribution of water inside the porous medium can be observed along the primary direction of the flow, as shown in Fig. 8 . The liquid phase does not penetrate uniformly, but rather the liquid front presents a smooth shape, as indicated by the decrease in the saturation along the thickness of the medium. This decrement indicates that after the breakthrough time t * ≥ 0.5 , the liquid phase exhibits a compact uniform front at the GDL entrance, whereas farther along the GDL, it is redistributed along different preferential flow paths (see also the right panels of Fig. 8) . These preferential finger-shaped paths are fluid structures typically observed during twophase fluids invasion in porous media. For instance, viscous fingering refers to the situation when a less dense fluid displaces a denser one and, for sufficiently high invasion velocities, viscous forces destabilise the invading front and sustain the formation of fingers at the interface (Méheust et al. 2002) . On the contrary, when a denser fluid invades a medium initially filled with a less dense one, the viscous forces tend to stabilise the front (Aker et al. 2000) . In fuel cells, the GDL is usually initially filled with a gaseous fluid (hydrogen or oxygen). As the cell starts operating, liquid water (more dense) is produced and invades the medium. In such a situation, with M = ∕ g ≫ 1 , viscous forces contribute to stabilise the front while the distribution of pore-scale capillary thresholds determines the onset of fingers, a phenomenon known as capillary fingering. The occurrence of capillary fingering is known to depend also on other physical parameters, such as the capillary number and the hydrophobicity of the material, the latter determining the shape of the liquid-gas curvature at the pore throats (Lenormand et al. 1988; Cottin et al. 2010) .

For instance, the transition from uniform infiltration to capillary fingering has been observed in media with uniform randomly distributed pore sizes, when changing the liquid-gas-solid contact angle, from hydrophilic to hydrophobic (Jung et al. 2016) , or when controlling the competition between pore size gradients and medium disorder (Lu et al. 2019 ). Here, we observe a similar transition, but strongly governed by the porous microstructure. The high value of the capillary pressure recorded (see Fig. 7 ) can be interpreted as the result of an "effective" contact angle, which, rather than by the chemical properties of the solid surface (have in mind that in our simulations, the equilibrium contact angle is eq = 90 • ), is induced by the combination of strong dynamics of the capillary pressure and the presence of significant local constrictions of the flow path, which leads to a high curvature of the two-phase interface. Subsequently, the interface experiences a mechanical instability and a successive sudden invasion event occurs across certain pore throats. Such an effective contact angle appears then to induce noncooperative burst instabilities that dominate the interface advance to form ramified fluid structures, as observed for porous media with contact angles ≳ 120 • (Jung et al. 2016) . In other words, the invading fluid can be considered as "non-wetting" as we observe the occurrence of capillary fingering even with non-hydrophobic surfaces. The mechanism of microstructure-induced capillary thresholds counteracts the viscous forces to destabilise the front as also confirmed by the scaling of the pressure drop at the front width observed in Fig. 7 . In Sect. 3.4, we will quantify such an effect through the formulation of an effective model. A similar physical approach, which makes use of an effective contact angle, has been used for example for the understanding and characterisation of two-phase flow dynamics in soils, see e.g. Wessel (1988) . In Fig. 9 , we report the variation of planar porosity along the domain. The variation of the porosity suggests the presence of pore constrictions followed by pore enlargements. Such a microstructural traits is likely to induce a positive permeability gradient along the flow direction, thus to destabilise the front after it invades the first pore constriction at S ∼ 2 , and to induce the capillary fingering observed in Fig. 8 and the scaling depicted in Fig. 7 (Xu et al. 1998) . This mechanism can also be intuitively understood by considering that the viscous pressure drop is inversely proportional to permeability (Darcy's law) and that thus a sufficiently high increase in the permeability along the flow direction can cause a decrease in the viscous drop along the front. It also confirms the important role that pore-scale statistical variations along the medium thickness can have in determining a fluid displacement, given the lack of separation of length scales characterising the pore and minor medium dimensions. Notably, we observe a minimum value of the porosity at x∕L pm ∼ 0.25 , an observation that sustains our hypothesis that a significant change in the two-phase flow behaviour coincides with the occurrence of a change in the pores morphology (see also Fig. 7) . In particular, the pore radius PDFs computed at different cross sections depicted in the bottom of Fig. 9 show the more frequent occurrence of small pore throats at x∕L pm ∼ 0.2 . Such small pores of radius ∼ 15 μm are significantly smaller than the average pore radius r p = 40 μm . We suggest that the presence of pores with a small radius in a certain region of the GDL can be due to an inhomogeneous distribution of the binder, which is used to bound together the polymer matrix during GDL fabrication. A bimodal distribution of pore sizes has been for instance observed in reconstructed GDLs and associated with substantial differences in the binder amount in different (a) (b) (c) Fig. 9 Porosity, , as a function of GDL normalized thickness x/L (top) and the probability distribution functions (PDFs) of pore radii at three cross sections (bottom) for the three samples. In the bottom figure, the vertical and horizontal axes are the relative frequency and pore radius, respectively. The PDFs are presented at x∕L pm ∼ 0.1 (black lines), x∕L pm ∼ 0.25 (blue lines) and x∕L pm ∼ 0.8 (red lines) regions (Zenyuk et al. 2016) , a difference that here we observe to be significant in the through-plane direction. The transition occurs thus both in space and time: liquid invasion occurs rapidly and uniformly for short times, until the front reaches x∕L pm ∼ 0.2 (just before the observed constriction at x∕L pm ∼ 0.25 ). Later in time and farther in space, the liquid is divided into clusters; infiltration occurs firstly on the large pores and then along various flow paths or fingers in an unstable manner.

The breakthrough time defines the moment when the first portion of liquid reaches the outlet. It is then interesting to observe the spatial configuration of the liquid phase at this instant, since it provides information about how the porous structure has affected the mass transport from the inlet to the outlet of the medium. We should also emphasise that the breakthrough time ( t * ∼ 0.5 ) occurs shortly after the transition between the ballistic and viscous-capillary regimes ( t * ∼ 0.3 ), an insight that suggests a fast dynamics of pore invasion after the liquid front overcomes the capillary barrier induced by the pore constriction at x∕L pm ∼ 0.2.

In Fig. 10 , the spatial configuration of the liquid fingers at various cross sections after the constriction is shown for one of the REVs at the breakthrough time t * ∼ 0.5 . After x∕L pm ∼ 0.2 , the average dimensionless diameter d f ∕(2r p ) and number n f of the fingers gradually decrease along the cross sections of the medium. This trend is quantified in Fig. 11 , where the latter quantities are plotted along the medium thickness. We firstly note that the initial part of the GDL ( x∕L pm < 0.2 ) is almost uniformly filled with water, since d 2 f ∕4 n f ∼ A . In this region, rather than fingers, we observe clusters of liquid separated by fibres, with the medium almost fully filled, as also shown in Fig. 8 , where S(x∕L pm < 0.2) ∼ 1.

Closer to the GDL outlet, the diameter of the fingers decreases from an average size of d f ∼ 160 μm = 4r p at x∕L pm ∼ 0.5 , to d f ∼ 80 μm = 2r p at x∕L pm ∼ 1 . In Fig. 11 , a decrease in the number of fingers along the GDL thickness is also observed. The quantification of the liquid structures thus offers us another important insight: at the breakthrough time, the distribution of liquid invasion gradually changes in space along Fig. 10 Liquid distribution at the breakthrough time at three slices along the streamwise direction for sample 3; a x∕L pm = 0.4 , b x∕L pm = 0.6 and c x∕L pm = 0.8 . The white domains represent liquid clusters (fingers), while the black ones are solids or gas. The simulations illustrate that the invading fingers of the liquid water become thinner as the liquid advances in the porous medium the medium, from a collaborative pore filling behaviour, with clusters of the size of few pores d f (x∕L pm ∼ 0.5) ∼ 4r p , to a preferential flow along a single pore, when d f (x∕L pm = 1) ∼ 2r p , as more and more pore throats contribute to impede liquid transport along the direction of the flow. Interestingly, this finding motivates a possible new design of GDLs, as characterized by an increasing size of the pore diameters, as a recent study also suggested (Balakrishnan et al. 2020 ). Such a design should promote the water transport in the GDL by decreasing the capillary pressure thresholds, induced by the pore throats, along the flow direction. However, such an innovative design should also be carefully tailored to take into account the effects of the spatial distribution of pore sizes (i.e. the amount disorder of the medium), a parameter which is found to be important for the control of capillary fingering phenomena (Lu et al. 2019) .

The observed spatial redistribution of the liquid phase from clusters to single invaded pores at the breakthrough time suggests that water primarily emerges from individual pores distributed on the GDL outlet surface. To better understand the characteristics of the structures characterising the liquid transport along such preferential paths, from the inlet to the outlet of the GDL, we quantify the shape and dynamics of the liquid fingers at the outlet cross section, after the breakthrough time, for t * ≥ 0.5 . The computed values of d f are reported in Fig. 11 , where the average distance between emerging fingers f is also computed, while in Fig. 12 the emerging fingers at the outlet cross section are displayed for a visual comparison. As previously suggested, the liquid phase emerges in the form of fingers with the size comparable with the average pore sample diameter ∼ 2r p . In the times following the breakthrough, the average size of the fingers slightly increases, as well as their number (not shown here), as more liquid reaches the outlet of the medium along a gradually increasing number of pores. The time-dependent increase in finger dimension on the outlet surface suggests a powerlaw scaling d f ∝ t * 1∕D , being D the fractal dimension of the emerging cluster. This observation is consistent with classic percolation-invasion studies that predict the more probable invasion site occurring at a distance ∝ t * 1∕D from the previous event, with D ∼ 1.82 (Furuberg et al. 1988; Wilkinson and Willemsen 1983) . A least-square fit would give us D ∼ 2.3 , even though we note that such a scaling is limited to just half a decade. On the contrary, we find that the occurence of emerging fingers on the GDL surface appears to be uncorrelated, since their average computed distance fluctuates around a constant value f ∼ 2r p , as reported in the bottom-right panel of Fig. 12. 

We have observed three regimes that characterise the mass transport within the GDL. Among them, the ballistic regime occurs at short times and only in the entrance region of the GDL. In this section, we aim at modelling the water transport from a macroscopic perspective, and we therefore delineate a model based on the long time regimes that dominate the great part of the two-phase flow dynamics in GDLs.

During the viscous-capillary stage, the liquid flows through an increasing number of preferential pores, and it experiences an increasing viscous resistance induced by the augmented liquid-solid interface. As observed in Fig. 6 , the liquid transport globally follows a power law S(t) ∝ t * 1∕2 until a significant amount of water has reached the outlet of the GDL and dS∕dt ∼ 0 . We can express the transport of water inside the GDL until this time ( t * ∼ 2 ) by means of a model based on Darcy's formulation. The latter assumes a linear relation between the pressure drop and the flow rate, as expressed in equation 10, which is quantified by the permeability of the medium. Darcy's equation describes the flow in terms of average values, which are obtained at length scales much larger than the structural heterogeneities of the porous materials. Therefore, a model based on Darcy's equation must inevitably rely on effective quantities that convey information from micro-to macroscales. To build up a Darcy-based model, we then extract the effective quantities that govern the dynamics of liquid infiltration in the capillary-viscous regime, where we observe a Darcian behaviour of the fluid, i.e. S(t) ∝ t * 1∕2 . In this regime, the flow dynamics is induced by the pressure drop −∇ x PL pm acting on the average penetration length h; the capillary pressure induced by the porous microstructure adds an additional resistance to the flow, which cannot be neglected, as discussed in Sect. 3.1. We can then define an effective permeability K eff that embeds the effect of the microstructure and rewrite equation 10 as:

where K eff = 1.2 10 −10 m 2 is the effective permeability extracted by fitting the values reported in Fig. 6 . The same equation can be rewritten in a dimensionless form as:

with K * eff = K eff ∕(r 2 p ∕4) = 0.30 , which points out that, from a macroscopic perspective, the liquid infiltration occurs similarly to the infiltration in a pipe with an effective radius r p,eff ∼ 0.55 r p ∼ 22 μm . The value of the effective radius, lower than half of the observed average finger diameter, indicates the adverse effects of the capillary pressure on the flow dynamics, which is then characterised by an increased resistance to the flow. Recalling the picture of an effective contact angle induced by the microstructure discussed in Sect. 3.2, we can estimate this angle from the value of the calculated effective radius. The liquid-gas curvature induced by the microstructure provides an additional resistance to the flow in the form of an effective capillary pressure P c,eff = −∇ x PL pm (1 − r 2 p,eff ∕r 2 p ) . This value ( P c,eff ∼ 0.7 P t ) is very close to the observed value of capillary pressure at long times (see e.g. Fig. 7) . From the effective capillary pressure, the effective contact angle is then calculated as: a value that we notice places the two-phase flow dynamics in the regime of noncooperative burst instabilities for two-phase flows in porous media (Jung et al. 2016) .

It is interesting to observe that by using the calculated effective contact angle or, equivalently, the effective pore radius, we determine the characteristic invasion velocity as U c,eff = r 2 p,eff ∕r 2 p ∼ 0.3 U c and the characteristic time for liquid to reach the GDL outlet as t c,eff ∼ 3 t c . Indeed, at around this time, for t * > 2 , the mechanism changes by which water is transported, since enough preferential flow paths connecting the inlet and outlet of the medium are formed. In this later stage, the liquid is transported preferentially along the liquid fingers and the total saturation of the medium approaches a constant value S ∼ 0.75 . We can estimate the volumetric flow rate at the outlet in this later stage by means of the product between the average liquid velocity u f = − d 2 f ∇ x P∕(16 w ) and the liquid area at the outlet A f = n f d 2 f ∕4 . The calculation leads to a volumetric flow rate q out = u f A f ∼ 5 cm 3 ∕s , and a capillary number Ca out (t) ∼ 0.5 , at the time t * = 2.5 . These are useful quantities for the prediction of water emerging from the GDL and transported within the fuel cells gas channels.

In a macroscopic effective model that predicts the two-phase mass transport in fuel cells, Eq. 11 can describe the motion of the fluid in the bulk of the GDL, via an effective contact angle given in Eq. 13. Furthermore, the calculated flow rate and the capillary number at the outlet can serve as microscopic input parameters for the boundary conditions of transport at the GDL-gas channels interface. As a final note, we point out that the value of Ca out (t) , which is close to unity, underlines the fast dynamics that water experiences when the liquid reaches the GDL outlet and, in turn, the prominent role that inertia plays in the formation of liquid droplets emerging from the GDL surface.

Because of the urgent need for sharing research results and information in a global effort to deal with the current COVID-19 pandemic, we finally want to highlight the significance of the present results for the understanding and development of effective personal protective equipment (PPE). Among them, face masks play an important role for source control (used by infected persons) or prevention of COVID-19 (used by healthy persons), whether they are surgical masks used in health care environments or nonmedical masks used in the community setting (World Health Organization 2020). It is crucial to optimise the face masks effectiveness in filtering pathogen-bearing respiratory droplets, without significantly compromising their comfort and accessibility to the public (cost). A possible solution lies in the identification of the optimal microstructural design of the masks' filter layers, which in turn facilitates both the decision-making process to issue recommendations with regards to wearing certain types of masks, and the production of new low-cost effective materials. The fibrous porous structure of the gas diffusion layer investigated here strongly resembles many common fabrics materials that can be used as PPE, being composed of layers of interconnecting fibres. For the sake of clarity, we here define "large respiratory droplets" as the large pathogen-bearing droplets characterised by diameters d b significantly larger than the mean pore size of the filter layer, that is, with d b ≫ 2r p . In such a way, they can be well represented by a liquid front impacting a surface area of the size of the here investigated samples, A ≈ (20r p ) 2 . For instance, for face masks with a mean pore size of ≈ 50 m, the size of the large respiratory droplets would be ≈ 500 m, a value that characterises a small yet important portion of cough-generated droplets, since this portion is more likely to contain pathogenic organisms (Duguid 1946) . The present results unveil three very important insights about the transport of large pathogen-bearing through filter layers of face masks.

(i) The fibrous microstructure, even if not hydrophobic, induces an effective resistance to the their transport. Such an effective resistance, which we quantified via an effective contact angle in Eq. 13, is particularly important when the expiratory event induces a high value of the pressure difference between the two sides of the mask, and the resulting pressure drop is comparable with the characteristic capillary pressure of the microstructure P t = 2 ∕r p . Such a scenario is for instance encountered when a mask-wearing subject coughs and the oesophageal pressure exhibits very high values, in the order of tens of kilopascal (Lavietes et al. 1998 ). (ii) Since, as we have observed, pore throats incrementally contribute to impede liquid transport, the design of a layered porous structure with a decreasing size of the pore diameters along the flow direction can augment the filtration mechanism of large respiratory droplets, without necessarily increasing the thickness of the filter layer. Since the filtration mechanism in protective masks can be important both in inward and outward directions (e.g. during inhaling and exhaling), such a gradient in pore sizes could be designed as function of thickness relative to the middle plane of the porous material. A similar design should however also take into account the microstructure-induced pressure drop and thus the provided comfort (breathability). The existence of an effective capillary pressure opposing the liquid transport in neutrally wetted porous media and the expected positive effect of microstructures with variable pore sizes can possibly explain why layered designs of commonly available fabric materials, devoid of protective hydrophobic layers, have been found to offer protection about equivalent to the filtration offered by 5-layer N95 masks (Lustig et al. 2020 ). (iii) Finally, we also observe the emergence of the liquid phase from the outer surface of the porous medium, in the form of fingers with diameter sizes comparable with the mean pore size of the material. These fingers are rather homogeneously spatially distributed and their characteristic velocity is quantified via Ca out ∼ 0.5 . Such an information is important for both the possible targeted design of masks with acceptable limits in the size of transmitted droplets, and for the delineation of health risk assessments related to the expected size and velocity of pathogen-bearing droplets expelled from face masks.

We have presented lattice Boltzmann simulations of two-phase flows performed in a CT reconstructed fuel cell gas diffusion layer (GDL). We have been able to identify the porescale characteristics of the two-phase flow dynamics and spatial distribution within the porous medium, with a micrometre accuracy. We have sampled three different regions of the reconstructed material to allow a statistically robust analysis. We have identified three regimes that characterise the two-phase flow dynamics in GDLs. At short times, for t * ≤ 0.3 , the liquid phase rapidly invades the medium rather uniformly with saturation following a ballistic behaviour S(t) ∝ t * 2 . When the liquid front reaches the first significant flow constriction, represented by an increase in the number of small pore throats, it exhibits a different behaviour. We have observed that to such a change in the porous morphology (presumably related to a local increase in the binder in the gas diffusion layer) corresponds a maximum value of the computed capillary pressure. The capillary pressure is found to act adversely to the flow, being induced by a convex shape of the liquid-gas interface when approaching the pore throats. Such an effect is observed in a three-phase (liquid-gas-solid) system with a neutral contact angle, and it is therefore purely induced by the interplay between flow dynamics and microstructure. From this moment, at intermediate times 0.3 < t * ≤ 2.5 , the liquid phase breaks through certain pores and forms preferential flow paths of the size of one-two pores. We have observed the liquid phase to be transported mainly through the mechanism of capillary fingering during this stage, with an increasing number of preferential flow paths contributing to the global transport. During this stage, the invasion dynamics exhibits an unstable displacement and a kinetic roughening of the front, which is sustained by the sudden enlargement of the pore throats. At the same time, we have observed a gradual decrease in the liquid velocity as a consequence of the gradually increasing liquid-solid surface that provides viscous resistance, so that the saturation follows a power law S(t) ∝ t * 1∕2 . At long times, for t * > 2.5 , when a sufficient amount of liquid has reached the outlet, the two-phase fluid structure stabilises and the liquid flows exclusively following the previously formed preferential flow paths.

On the basis of these observations, we have proposed a model to mimic and predict the liquid transport along the gas diffusion layer. By making use of an effective contact angle that quantifies the effect of the measured capillary pressure, we have been able to build up a simple equation of motion, based on Darcy's equation, that well approximates the transport of water within the fuel cell material. We argue in the paper that such an effective contact angle does not only explain the noncooperative pore filling and unstable front displacement behaviours observed, but also the strong intermittent dynamics of two-phase flows and pressure signals often measured in fuel cells during flooding conditions.

In the last part of the manuscript, we have highlighted the significance of these results for the optimal design of filter layers of protective masks. We have explained the mechanisms governing the transport of large respiratory droplets, produced by intensive expiratory events, through fibrous porous microstructures. We have discussed how such forms of pathogens transmission can be efficaciously mitigated also through the use of nonmedical face masks made of common fabrics materials, provided that the proper characteristics of the different layers that compose them are selected.

Funding Open access funding provided by Chalmers University of Technology. The research project to which this work belongs has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant agreement No 790744 and from the Chalmers' Areas of Advance Transport. The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC).

An excerpt of data is available on the EUDAT B2SHARE platform at https ://b2sha re.eudat .eu/recor ds/32e06 d2305 39414 d8afb 146f5 95ca6 ea under a Creative Common Attribution-NonCommercial-ShareAlike 4.0 International License.

The code used for the numerical simulations is accessible at https ://gitla b.com/dario m/ lbdm under a Creative Common Attribution-NonCommercial-ShareAlike 4.0 International License.

Conflict of interest The authors declare that they have no conflict of interest.

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Viscous stabilization of 2d drainage displacements with trapping

A review of cell-scale multiphase flow modeling, including water management, in polymer electrolyte fuel cells

Designing tailored gas diffusion layers with pore size gradients via electrospinning for polymer electrolyte membrane fuel cells

Liquid water visualization in PEM fuel cells: a review

Water management in PEM fuel cells

Unravelling the surface chemistry and structure in highly active sputtered pt3y catalyst films for the oxygen reduction reaction

Batteries and fuel cells for emerging electric vehicle markets

The reactant starvation of the proton exchange membrane fuel cells for vehicular applications: a review

Pore-scale flow and mass transport in gas diffusion layer of proton exchange membrane fuel cell with interdigitated flow fields

A critical review of the pseudopotential multiphase lattice Boltzmann model: methods and applications

Water balancing

Drainage in two-dimensional porous media: from capillary fingering to viscous flow

Dynamic modeling of the effect of water management on polymer electrolyte fuel cells performance

A new boundary condition for three-dimensional lattice Boltzmann simulations of capillary filling in rough micro-channels

Lattice Boltzmann simulation of liquid water transport inside and at interface of gas diffusion and micro-porous layers of pem fuel cells

The size and the duration of air-carriage of respiratory droplets and droplet-nuclei

Modelling the performance of the cathode catalyst layer of polymer electrolyte fuel cells

Quantifying water transport in anion exchange membrane fuel cells

Europe's moment: Repair and prepare for the next generation COM/2020/456 final

Practical cone-beam algorithm

The transition from inertial to viscous flow in capillary rise

Dynamics of invasion percolation

Lattice Boltzmann simulation of water and gas flow in porous gas diffusion layers in fuel cells reconstructed from micro-tomography

Pore network modeling of fibrous gas diffusion layers for polymer electrolyte membrane fuel cells

Discrete lattice effects on the forcing term in the lattice Boltzmann method

Lattice Boltzmann simulations of water transport in gas diffusion layer of a polymer electrolyte membrane fuel cell

A review of water management in polymer electrolyte membrane fuel cells

Wettability controls slow immiscible displacement through local interfacial instabilities

Mass transfer in fibrous media with varying anisotropy for flow battery electrodes: direct numerical simulations with 3d x-ray computed tomography

Airway dynamics, oesophageal pressure and cough

Numerical models and experiments on immiscible displacements in porous media

A review of water flooding issues in the proton exchange membrane fuel cell

Effects of hydrophobic polymer content in GDL on power performance of a PEM fuel cell

Controlling capillary fingering using pore size gradients in disordered media

Water management studies in PEM fuel cells, part iii: dynamic breakthrough and intermittent drainage characteristics from GDLs with and without MPLs

Effectiveness of common fabrics to block aqueous aerosols of virus-like nanoparticles

Solute transport and reaction in porous electrodes at high Schmidt numbers

Electrode for proton exchange membrane fuel cells: a review

Interface scaling in a two-dimensional porous medium under combined viscous, gravity, and capillary effects

Liquid water breakthrough pressure through gas diffusion layer of proton exchange membrane fuel cell

Effective diffusivity and water-saturation distribution in single-and two-layer PEMFC diffusion medium

Impact of heat and water management on proton exchange membrane fuel cells degradation in automotive application

Signature analysis of two-phase flow pressure drop in proton exchange membrane fuel cell flow channels

Effect of cathode GDL characteristics on mass transport in PEM fuel cells

On the impact of porous media microstructure on rainfall infiltration of thin homogeneous green roof growth substrates

An automated simple algorithm for realistic pore network extraction from micro-tomography images

Direct simulations of two-phase flow on micro-CT images of porous media and upscaling of pore-scale forces

Simulation of liquid water breakthrough in a nanotomography reconstruction of a carbon paper gas-diffusion layer

Scale effect and two-phase flow in a thin hydrophobic porous layer. Application to water transport in gas diffusion layers of proton exchange membrane fuel cells

Dynamic fluid connectivity during steadystate multiphase flow in a sandstone

Lattice Boltzmann model for simulating flows with multiple phases and components

Relationships between water management, contamination and lifetime degradation in PEFC

Capillary driven flow in circular cylindrical tubes

The lattice Boltzmann equation: for fluid dynamics and beyond

Water transport and absorption in pharmaceutical tablets-a numerical study

Disorder characterization of porous media and its effect on fluid displacement

Dynamic capillary pressure in porous media: origin of the viscousfingering length scale

On using the effective contact angle and the water drop penetration time for classification of water repellency in dune soils

Invasion percolation: a new form of percolation theory

World Health Organization: Advice on the use of masks in the context of covid-19: interim guidance

Invasion percolation with viscous forces

Two-spheres method for evaluating the metrological structural resolution in dimensional computed tomography

Gas-diffusion-layer structural properties under compression via X-ray tomography

Comprehensive comparison of pore-scale models for multiphase flow in porous media