Steady-State Water Drainage by Oxygen in Anodic Porous Transport Layer of Electrolyzers: A 2D Pore Network Study

: Recently, pore network modelling has been attracting attention in the investigation of electrolysis. This study focuses on a 2D pore network model with the purpose to study the drainage of water by oxygen in anodic porous transport layers (PTL). The oxygen gas produced at the anode catalyst layer by the oxidation of water ﬂows counter currently to the educt through the PTL. When it invades the water-ﬁlled pores of the PTL, the liquid is drained from the porous medium. For the pore network model presented here, we assume that this process occurs in distinct steps and applies classical rules of invasion percolation with quasi-static drainage. As the invasion occurs in the capillary-dominated regime, it is dictated by the pore structure and the pore size distribution. Viscous and liquid ﬁlm ﬂows are neglected and gravity forces are disregarded. The curvature of the two-phase interface within the pores, which essentially dictates the invasion process, is computed from the Young Laplace equation. We show and discuss results from Monte Carlo pore network simulations and compare them qualitatively to microﬂuidic experiments from literature. The invasion patterns of di ﬀ erent types of PTLs, i.e., felt, foam, sintered, are compared with pore network simulations. In addition to this, we study the impact of pore size distribution on the phase patterns of oxygen and water inside the pore network. Based on these results, it can be recommended that pore network modeling is a valuable tool to study the correlation between kinetic losses of water electrolysis processes and current density.


Introduction
With innovations in the energy sector and a need for clean fuel, research is in progress to exploit the potential of hydrogen as an efficient energy source. Exemplarily, fuel cells can utilize hydrogen to produce electricity, and it can also serve as a fuel for internal combustion engines [1]. For the production of hydrogen, electrolyzer technology serves as a very promising and viable option. The purity of the produced hydrogen can be almost 100 vol % [2]. This way, it can be integrated with other renewable resources to offer a broad range of ecologically clean methods for hydrogen production [3][4][5]. Shortly, electrolyzers and fuel cells will be able to help alleviate the effects of fossil and nuclear fuel consumption [2,6,7]. This, however, implies efficient performance of electrolyzers.
A trade-off between the production rates and the efficiency of an electrolyzer is still to be met. For this technology to be commercial, the cost and hence the efficiency needs to be optimized. Among the electrolyzer technologies, polymer electrolyte membrane (PEM) electrolyzers have an edge over the other varieties because of high energy efficiency and compact design [8][9][10]. Power needed for the electrolyzer operation can be obtained from renewable sources too and such systems have also been gaining a lot of attention recently [11]. Polymer electrolyte membrane electrolytic cells (PEMECs) are very common when coupling the technology with other renewable sources like solar or wind energy.
High current densities are necessary in order to obtain high production rates. The problem is that the high current density operation results in a decrease of electrical efficiency. This is for example, reported in [12][13][14][15][16] and it is mainly explained with the kinetic losses associated with the mass transfer resistances through the porous transport layers (PTL) [17]. The counter-flow of the two phases causes hindrance against each other, and this mass transfer resistance causes a decrease in the overall electrolyzer performance [18,19].
It is strongly assumed that the high oxygen production rates achieved at the high current densities allow for the formation of gas bubbles that can invade the PTL and partially drain the water [20,21]. The size of gas bubbles increases with the increase in current density [12,15,22]. This leads to the formation of gas fingers penetrating the PTL from the anode catalyst layer, where the oxygen is produced, to the water supply channel from where the oxygen is removed [8]. The development of these gas fingers obviously results in the reduction of the overall water loading of the PTL. This can severely affect the water permeability, especially if the surface saturation of the PTL with water is significantly reduced by gas-filled pores [23]. On the other hand, the efficiency of the oxygen transport in the opposite direction depends mainly on the tortuosity of the evolving gas branches.
According to the work done by Suerman et al. [24], 25% of the total cell overpotential could be contributed to the mass transport losses and these losses are mostly credited to the oxygen withdrawal from the catalyst surface and the PTLs. Yigit et al. [25] reported that at current densities less than 0.7 A/cm 2 , the hydrogen production rates were very low and at values above 1 A/cm 2 the electrical efficiency decreased. Larger pores or high porosity values could easily mitigate this mass transport problem on the side of the gas phase, but it would also decrease the electron transport and affect the efficiency [8]. In contrast, small porosity values would hinder the gas removal and increase the entrapped water amount within the catalyst layer, and thus, decrease the rate of reaction [8]. For this purpose, current research aims at the structural optimization of PTLs with respect to efficient mass and electron transfer.
Various experimental methods [12,[14][15][16]20,[26][27][28] and modeling approaches [29][30][31][32][33] are already established to analyze the key factors of the PTL, such as flow regimes, structure, porosity, pore size distribution (PSD), permeability, corrosion resistance, and electrical conductivity. From these studies, the mass transfer limitations discussed before are generally either assigned to the flow regimes of the gas phase (e.g., in Dedigama et al. [12] and Zhang et al. [34]) or to the structure of the PTL. The latter has been studied in [14,15,26,29] (Table 1). In Ito et al. [15] an experimental study on a 27 cm 2 PEM electrolyzer cell was presented that investigates the influence of porosity and pore diameter. In this study, Ti-felt and Ti-powder prepared PTLs with different porosities and different pore structures were used. According to this study, the optimum pore diameter would be 10 µm, whereas no significant effect was seen for a porosity value greater than 50%. Grigoriev et al. [14] estimated an optimum value of 12-13 µm and 30%-50% as the optimum value of porosity using polarization curves for Ti-sintered PTLs with different properties. Findings presented in Kang et al. [26] in an experimental study based on thin/well tunable PTLs in a conventional cell suggested high porosity values and small pore sizes. Ojung et al. [29] used a semi-empirical model in their investigations to study PEM cell without flow channels. They varied pore sizes between 5-30 µm. They observed a decrease in performance at 5 µm and greater than 11 µm there was no significant improvement shown by their model. They also concluded that in a system without flow channels, porosity would not influence the performance at fixed pore diameter and an optimum value of 60% was observed. Hwang et al. [35] also showed with the experiments on reversible fuel cells using Ti-felt that for mean pore sizes around 30-50 µm porosity is an insignificant factor. Besides this, explanations for the structure dependence of the electro activity can also be derived from the consideration of flow regimes. Dedigama et al. [12] studied the flow regime within the PTL using electrochemical impedance spectroscopy (EIS) and thermal imaging. They found a reduction in the mass transfer limitation when passing from the bubble to the slug flow regime. In agreement with that, Zhang et al. [34] observed a decrease in the efficiency with an increase of the mass flow rate of water. Aubras et al. [31] showed that the porosity of the anode PTL affects the non-coalesced bubble regime. According to this study, higher porosity can enhance coalescence of oxygen bubbles and increase the performance of electrolyzer. Han et al. [36] also showed an increase in performance linked to an increase in porosity. In addition to that, more recent studies [33,37] revealed the importance of the interaction of the two involved hierarchical porous structures at the anode electrode, namely the PTL and the catalyst layer. Exemplarily, it is demonstrated by Lee et al. [37] using micromodel experiments that the pore sizes control the burst velocity of gas resulting in the application of a thin, low porosity region at the inlet in order to reduce the gas saturations in the PTL.
The majority of the data suggests that a relatively lower value of pore size (around 10 µm) is favorable and no significant conclusion can be drawn about the porosity value. Some authors suggest a high porosity value but others suggest that higher porosity would result in a slug flow, which can lead to inefficient mass transfer and a decrease in efficiency. On the contrary, coalesced bubble regime is also suggested to enhance the performance. In our view, other properties besides mean pore size and porosity might influence the invasion process more significantly. In this paper, we aim at approaching open questions by means of pore network (PN) modeling. In detail, we will consider the role of PSD, which has a higher significance for the invasion in PTL than porosity.
From the above discussions, it can be concluded that advanced manufacturing processes are required to tune the PTL performance. The reader may refer to [26,[38][39][40] for various examples of PTL production techniques that optimize the structure in terms of electrical efficiency and also in terms of material consumption and costs. According to Lettenmeier et al. [39] vacuum plasma spraying can be used to manufacture a PTL with a gradient in pore size along the thickness. It is possible to obtain an average pore size of 10 µm close to the bipolar plate and 5 µm at the electrode side, with the help of this technique (Figure 1). This coating technique can also be used to alter other properties of the PTLs suiting to the need. In a very recent study, Lee et al. [33] showed the effect of porosity gradient on the performance of the electrolyzer. The low porosity region was towards the membrane side and the high porosity region on the water inlet side (Figure 1). They observed high water permeation despite high oxygen saturations. Mo et al. [38] used electron beam melting (EBM) to mitigate the cost and manufacturing issues of the PTL. They showed 8% improvement in the performance compared to the conventionally woven PTLs. They obtained smooth surfaces at both ends of the PTL, thereby reducing the contact resistance between PTL and catalyst layer. Schuler et. al. [41] verified this impact of the interface between PTL and the catalyst layer clearly in their work.  [33]. (b) Pore size gradient within PTL as in the study from [39].

Pore Network Models
The available methods for studying two-phase flow can be divided into continuum and discrete models. The continuum models are usually formulated for macroscale continuous phases employing effective parameters and are thus not suitable for explaining the discrete processes that occur at pore scale. So, in order to gain a deeper understanding of the invasion processes under the action of capillarity, viscous or gravity forces, pore scale models are usually preferred. Pore network models (PNMS) are a type among various pore scale models available, e.g. Lattice Boltzmann models which are mostly available on the scale of one pore. PNMs are discrete models that represent the pore space by a lattice of pores and throats (e.g., [42][43][44][45][46]). In comparison to other methods, which are computationally more expensive in terms of discretization of the physical domain and solving of the governing equations, PNMs can be used to study larger systems computationally more efficiently. Besides fundamental physical studies, PNMs are also used for the parameterization of macroscopic models, e.g. to predict the capillary pressure curve, relative permeability curves, as well as saturation curves [46].
PNMs are generally distinguished between quasi-static and dynamic models [47]. For the simulation of the steady-states in capillary-dominated applications [33,[48][49][50][51] , quasi-static models are used. In the absence of dynamic effects, e.g. driven by viscous forces, the entry capillary pressure of pores and throats controls the displacement of liquid (drainage) or gas (imbibition) in such applications [52]. In this work, the quasi-static model from [46] is applied for the simulation of the drainage of water by oxygen.
The objective of this study is to achieve a fundamental understanding of gas and liquid transport within the PTLs and the pore structure dependence of the invasion patterns. A regular 2D network of pores and throats is used to illustrate the porous PTL. The void space of this network comprises of spherical pores and cylindrical throats. Invasion percolation rules for quasi-static capillary invasion are employed to simulate the displacement of water by oxygen. The displacement mechanism is controlled by the capillary pressure curves of water that are obtained from the Young-Laplace equation:

2
(1) where Pl is the liquid pressure, Patm is the atmospheric pressure, σ in kg s −2 is the surface tension, θ is the wetting angle, and r is the radius of the channel. Invasions or displacements occur when they are energetically more favorable. More clearly, the pressure difference between the wetting fluid (water) and the non-wetting fluid (oxygen) leads to the formation of a curvature with a radius depending on the pressure difference. In the case of drainage, the non-wetting phase is the one with higher pressure to be forced through the porous structure. The incremental increase of the pressure inside the gas phase results in the invasion of the liquid filled

Pore Network Models
The available methods for studying two-phase flow can be divided into continuum and discrete models. The continuum models are usually formulated for macroscale continuous phases employing effective parameters and are thus not suitable for explaining the discrete processes that occur at pore scale. So, in order to gain a deeper understanding of the invasion processes under the action of capillarity, viscous or gravity forces, pore scale models are usually preferred. Pore network models (PNMS) are a type among various pore scale models available, e.g., Lattice Boltzmann models which are mostly available on the scale of one pore. PNMs are discrete models that represent the pore space by a lattice of pores and throats (e.g., [42][43][44][45][46]). In comparison to other methods, which are computationally more expensive in terms of discretization of the physical domain and solving of the governing equations, PNMs can be used to study larger systems computationally more efficiently. Besides fundamental physical studies, PNMs are also used for the parameterization of macroscopic models, e.g., to predict the capillary pressure curve, relative permeability curves, as well as saturation curves [46].
PNMs are generally distinguished between quasi-static and dynamic models [47]. For the simulation of the steady-states in capillary-dominated applications [33,[48][49][50][51], quasi-static models are used. In the absence of dynamic effects, e.g., driven by viscous forces, the entry capillary pressure of pores and throats controls the displacement of liquid (drainage) or gas (imbibition) in such applications [52]. In this work, the quasi-static model from [46] is applied for the simulation of the drainage of water by oxygen.
The objective of this study is to achieve a fundamental understanding of gas and liquid transport within the PTLs and the pore structure dependence of the invasion patterns. A regular 2D network of pores and throats is used to illustrate the porous PTL. The void space of this network comprises of spherical pores and cylindrical throats. Invasion percolation rules for quasi-static capillary invasion are employed to simulate the displacement of water by oxygen. The displacement mechanism is controlled by the capillary pressure curves of water that are obtained from the Young-Laplace equation: where P l is the liquid pressure, P atm is the atmospheric pressure, σ in kg s −2 is the surface tension, θ is the wetting angle, and r is the radius of the channel. Invasions or displacements occur when they are energetically more favorable. More clearly, the pressure difference between the wetting fluid (water) and the non-wetting fluid (oxygen) leads to the formation of a curvature with a radius depending on the pressure difference. In the case of drainage, the non-wetting phase is the one with higher pressure to be forced through the porous structure. The incremental increase of the pressure inside the gas phase results in the invasion of the liquid filled pores and throats, with liquid pressures depending on their radius and wettability (Equation (1)). This means, that the stepwise invasion of the interconnected pore space results in the formation of distinct invasion patterns that depend on the PSD and the connectivity of pores and throats.
In such networks, the porosity can be increased basically in two ways, namely by (i) increasing the number of throats of the same dimensions, and (ii) changing the distribution of the throat sizes ( Figure 2). In the example illustrated in Figure 3, porosity and mean throat diameter are kept constant and only the distribution varies following the invasion pressure curve computed using Equation (1). As can be seen, the differences in the structural organization of the PN are expected to result in different overall liquid saturations and different gas-liquid distributions [53]. The larger pores and throats are invaded by the gas phase while the smaller ones remain liquid saturated. pores and throats, with liquid pressures depending on their radius and wettability (Equation (1)). This means, that the stepwise invasion of the interconnected pore space results in the formation of distinct invasion patterns that depend on the PSD and the connectivity of pores and throats. In such networks, the porosity can be increased basically in two ways, namely by i) increasing the number of throats of the same dimensions, and ii) changing the distribution of the throat sizes ( Figure 2). In the example illustrated in Figure 3, porosity and mean throat diameter are kept constant and only the distribution varies following the invasion pressure curve computed using Equation (1). As can be seen, the differences in the structural organization of the PN are expected to result in different overall liquid saturations and different gas-liquid distributions [53]. The larger pores and throats are invaded by the gas phase while the smaller ones remain liquid saturated.

Model Description
The model is comprised of two parts; part one includes the determination of the data structure for the definition of the geometry of the void and solid space, and the second part contains the equations of the drainage algorithm and the cluster labeling ( Figure 4). The data structure contains the information about the connections between the pores and throats in the network ( Figure 5). This information is used in the drainage algorithm for the stepwise calculation of the successive invasion. The Hoshen-Kopelman algorithm [54,55] is then applied to identify invading and isolated liquid clusters. pores and throats, with liquid pressures depending on their radius and wettability (Equation (1)). This means, that the stepwise invasion of the interconnected pore space results in the formation of distinct invasion patterns that depend on the PSD and the connectivity of pores and throats. In such networks, the porosity can be increased basically in two ways, namely by i) increasing the number of throats of the same dimensions, and ii) changing the distribution of the throat sizes ( Figure 2). In the example illustrated in Figure 3, porosity and mean throat diameter are kept constant and only the distribution varies following the invasion pressure curve computed using Equation (1). As can be seen, the differences in the structural organization of the PN are expected to result in different overall liquid saturations and different gas-liquid distributions [53]. The larger pores and throats are invaded by the gas phase while the smaller ones remain liquid saturated.

Model Description
The model is comprised of two parts; part one includes the determination of the data structure for the definition of the geometry of the void and solid space, and the second part contains the equations of the drainage algorithm and the cluster labeling ( Figure 4). The data structure contains the information about the connections between the pores and throats in the network ( Figure 5). This information is used in the drainage algorithm for the stepwise calculation of the successive invasion. The Hoshen-Kopelman algorithm [54,55] is then applied to identify invading and isolated liquid clusters.

Model Description
The model is comprised of two parts; part one includes the determination of the data structure for the definition of the geometry of the void and solid space, and the second part contains the equations of the drainage algorithm and the cluster labeling ( Figure 4). The data structure contains the information about the connections between the pores and throats in the network ( Figure 5). This information is used in the drainage algorithm for the stepwise calculation of the successive invasion. The Hoshen-Kopelman algorithm [54,55] is then applied to identify invading and isolated liquid clusters.

Network Generation
Pore and throat radii ( , ) are randomly distributed around a given average value with defined standard deviations. The geometrical arrangement of throats and neighbor pores with the relevant geometry parameters is illustrated in Figure 6.

Invasion Algorithm
After the geometric parameters are specified, active, i.e. invading clusters with their menisci are identified and the maximum liquid pressure is computed within the active clusters using Equation (1). The algorithm then selects the largest accessible throat or pore for gas invasion following the rules

Network Generation
Pore and throat radii ( , ) are randomly distributed around a given average value with defined standard deviations. The geometrical arrangement of throats and neighbor pores with the relevant geometry parameters is illustrated in Figure 6.

Invasion Algorithm
After the geometric parameters are specified, active, i.e. invading clusters with their menisci are identified and the maximum liquid pressure is computed within the active clusters using Equation (1). The algorithm then selects the largest accessible throat or pore for gas invasion following the rules

Network Generation
Pore and throat radii (r p , r t ) are randomly distributed around a given average value with defined standard deviations. The geometrical arrangement of throats and neighbor pores with the relevant geometry parameters is illustrated in Figure 6.

Network Generation
Pore and throat radii ( , ) are randomly distributed around a given average value with defined standard deviations. The geometrical arrangement of throats and neighbor pores with the relevant geometry parameters is illustrated in Figure 6.

Invasion Algorithm
After the geometric parameters are specified, active, i.e. invading clusters with their menisci are identified and the maximum liquid pressure is computed within the active clusters using Equation (1). The algorithm then selects the largest accessible throat or pore for gas invasion following the rules

Invasion Algorithm
After the geometric parameters are specified, active, i.e., invading clusters with their menisci are identified and the maximum liquid pressure is computed within the active clusters using Equation (1). The algorithm then selects the largest accessible throat or pore for gas invasion following the rules specified in Vorhauer et al. [46]. As invasion proceeds, entrapped clusters are formed, which are permanently isolated due to the incompressibility of the fluids.

Cluster Labeling
During this process, labeling of liquid clusters is used to identify the pores and throats that are connected to each other and form a pathway for liquid and gas transport. At the start of the invasion process, the network is completely occupied by liquid, and there is only one cluster spanning the whole network and conducting the liquid phase. With initiation of invasion, numerous liquid clusters can form that are distinguished into liquid-conducting clusters (connected to bottom and top side of the PN) and isolated clusters. The Hoshen-Kopelman algorithm [54,55] is used for the labeling of these clusters. Clusters are reviewed and relabeled after each invasion step to update the connections between the pores and throats.
PN initially saturated with water.

4.
Oxygen is injected at the top side and water is removed from the bottom side.

5.
There is no mixing or diffusion between the two phases. 6.
Viscous, gravity and liquid film flow are neglected. 7.
Piston type throat invasion computed based on the Young-Laplace equation.

8.
No further invasion occurs after breakthrough of the gas phase.

Pore Network Simulation of Microfluidic Experiments
A pore and throat network was conceived using the parameters of microfluidic experiments from Arbabi et al. [20]. The PNM was constructed based on the image data extracted from Figure 7. The experimental image in Figure 7 is then used to compare the flow path of gas within the micromodel qualitatively with our own simulation results. In Figure 7, gas pores and throats are highlighted in blue (pores) and yellow (throats) while liquid pores are in red and liquid throats in black. In this investigation, we are interested if the PNM introduced above is able to predict the experimentally estimated invasion path. It is remarked that invasion is dictated by the interface curvature of pore and throat menisci, wherefore the pore and throat sizes are of interest here. They were determined from the experimental image and transferred into the data structure of the PNM. Although the pore sizes are significantly larger than the throat sizes in this example, pore invasion pressures were not matched with the pore sizes. Instead the pore sizes were randomly adjusted. As can be seen below, the simulation leads already to a very good agreement of the results. However, in a future study, it would be preferable to track experimentally the different invasion pressure thresholds of pores and throats based on the interface curvature of liquid menisci.
Processes 2020, 8, x FOR PEER REVIEW 7 of 18 specified in Vorhauer et al. [46]. As invasion proceeds, entrapped clusters are formed, which are permanently isolated due to the incompressibility of the fluids.

Cluster Labeling
During this process, labeling of liquid clusters is used to identify the pores and throats that are connected to each other and form a pathway for liquid and gas transport. At the start of the invasion process, the network is completely occupied by liquid, and there is only one cluster spanning the whole network and conducting the liquid phase. With initiation of invasion, numerous liquid clusters can form that are distinguished into liquid-conducting clusters (connected to bottom and top side of the PN) and isolated clusters. The Hoshen-Kopelman algorithm [54,55] is used for the labeling of these clusters. Clusters are reviewed and relabeled after each invasion step to update the connections between the pores and throats.

Pore Network Simulation of Microfluidic Experiments
A pore and throat network was conceived using the parameters of microfluidic experiments from Arbabi et al. [20]. The PNM was constructed based on the image data extracted from Figure 7. The experimental image in Figure 7 is then used to compare the flow path of gas within the micromodel qualitatively with our own simulation results. In Figure 7, gas pores and throats are highlighted in blue (pores) and yellow (throats) while liquid pores are in red and liquid throats in black. In this investigation, we are interested if the PNM introduced above is able to predict the experimentally estimated invasion path. It is remarked that invasion is dictated by the interface curvature of pore and throat menisci, wherefore the pore and throat sizes are of interest here. They were determined from the experimental image and transferred into the data structure of the PNM. Although the pore sizes are significantly larger than the throat sizes in this example, pore invasion pressures were not matched with the pore sizes. Instead the pore sizes were randomly adjusted. As can be seen below, the simulation leads already to a very good agreement of the results. However, in a future study, it would be preferable to track experimentally the different invasion pressure thresholds of pores and throats based on the interface curvature of liquid menisci.   The data of pore and throat sizes was implemented in the PNM to compute the successive invasion of the PN. The result of the simulation is shown in Figure 8. It is observed that the invasion pattern simulated with the PNM is identical with the experimental image (Figure 7). The perfect agreement reveals the suitability of the model structure and model assumptions and the ability of PNMs to predict the quasi-static drainage patterns. Though, in a future study it would be of interest, if also the stepwise invasion of the PN can be accurately predicted, not only for idealized 2D structures but also in larger and more realistic 3D structures. circles and throats. The image shows the steady-state invasion pattern after breakthrough of the gas phase from inlet (at the bottom) to water channel (at the top).
The data of pore and throat sizes was implemented in the PNM to compute the successive invasion of the PN. The result of the simulation is shown in Figure 8. It is observed that the invasion pattern simulated with the PNM is identical with the experimental image (Figure 7). The perfect agreement reveals the suitability of the model structure and model assumptions and the ability of PNMs to predict the quasi-static drainage patterns. Though, in a future study it would be of interest, if also the stepwise invasion of the PN can be accurately predicted, not only for idealized 2D structures but also in larger and more realistic 3D structures.

Impact of Pore Size Distribution in Monomodal PNs
The pore size properties of sintered PTL were used to study the effect of PSD on the invasion patterns and the steady-state saturation of the PTL at breakthrough of the gas phase. For this purpose, a pore network (PN) with the properties summarized in Table 2 were used: The standard deviation values were varied from 0.5 μm to 3 μm for a mean throat size of 17 μm (Figure 9), and pore sizes were used with a constant standard deviation value of 2 μm. Monte Carlo simulations were done for each data point so that the given gas saturation values in Figure 10 and Table 3 are an average of 20 simulations. In general, it is observed that the final gas saturation increases at breakthrough with widening of the radius distribution.

Impact of Pore Size Distribution in Monomodal PNs
The pore size properties of sintered PTL were used to study the effect of PSD on the invasion patterns and the steady-state saturation of the PTL at breakthrough of the gas phase. For this purpose, a pore network (PN) with the properties summarized in Table 2 were used. The standard deviation values were varied from 0.5 µm to 3 µm for a mean throat size of 17 µm (Figure 9), and pore sizes were used with a constant standard deviation value of 2 µm. Monte Carlo simulations were done for each data point so that the given gas saturation values in Figure 10 and Table 3 are an average of 20 simulations. In general, it is observed that the final gas saturation increases at breakthrough with widening of the radius distribution. The simulation results clearly reveal the impact of PSD on the final distribution of liquid and gas phase. While the average throat size was kept constant, the porosity of the PN increased slightly with increasing standard deviation of throat sizes (Table 3). Thus, following the literature discussions summarized above, it might be anticipated that higher porosities at constant mean throat or pore sizes are a result of an increasing standard deviation of pore and throat sizes in the referenced situations. As shown in Figure 10 and further analyzed below, the variation of PSD affects the invasion and thus the gas saturation. In detail, a higher gas saturation is obtained for broader PSDs. It is to be noted that higher PSDs than presented here can only be studied with a greater mean value of the throat size as will be discussed below.    The simulation results clearly reveal the impact of PSD on the final distribution of liquid and gas phase. While the average throat size was kept constant, the porosity of the PN increased slightly with increasing standard deviation of throat sizes (Table 3). Thus, following the literature discussions summarized above, it might be anticipated that higher porosities at constant mean throat or pore sizes are a result of an increasing standard deviation of pore and throat sizes in the referenced situations. As shown in Figure 10 and further analyzed below, the variation of PSD affects the invasion and thus the gas saturation. In detail, a higher gas saturation is obtained for broader PSDs. It is to be noted that higher PSDs than presented here can only be studied with a greater mean value of the throat size as will be discussed below.    The simulation results clearly reveal the impact of PSD on the final distribution of liquid and gas phase. While the average throat size was kept constant, the porosity of the PN increased slightly with increasing standard deviation of throat sizes (Table 3). Thus, following the literature discussions summarized above, it might be anticipated that higher porosities at constant mean throat or pore sizes are a result of an increasing standard deviation of pore and throat sizes in the referenced situations. As shown in Figure 10 and further analyzed below, the variation of PSD affects the invasion and thus the gas saturation. In detail, a higher gas saturation is obtained for broader PSDs. It is to be noted that higher PSDs than presented here can only be studied with a greater mean value of the throat size as will be discussed below.

Bimodal Pore Size Distributions
The previous results analyze the influence of the PSD on the example of monomodal PNs, i.e., only one peak in the histograms in Figure 9. As can be seen from Table 3, the porosity is only marginally affected in these cases. Due to this, the phase distribution patterns change only slightly in Figure 10. In reality, pore structures often obey bimodal PSDs, i.e., with two peaks in the histogram ( Figure 11) and with much stronger impact on porosity. The influence of macro pores is highlighted in Figure 3.
( Figure 11) and with much stronger impact on porosity. The influence of macro pores is highlighted in Figure 3.
As can be seen, for the monomodal PN in Figure 12, widespread invasion patterns with a high number of invasions is not achieved. This is in contrast to the bimodal PN in Figure 13. Monte Carlo simulations yielded an average gas saturation of 26% for the monomodal network and 38% for the bimodal network. This shows that widening of the PSD (Figure 9) and the introduction of macro pores (Figure 11),both, result in a change of the invasion process. This change is more significant in the second situation. While capillary fingering with narrow single gas branches is rather favored by narrow PSDs, widening of the invasion front with higher gas saturations is obtained by larger PSDs, and larger porosities. However, it can also be shown that the widening of the front can also be achieved when the porosity is kept constant and only the PSD is adjusted (Figure 13). The monomodal and bimodal networks used in simulations have a constant porosity of 71%. Figures 12  and 13 also show the saturation profiles of these simulations. They reveal the importance of the consideration of the PSD for characterization of the invasion process. Figure 11. Bimodal throat size distribution with standard deviation 2.0 μm for smaller (the first peak) and 2.5 μm for larger (the second peak) throats. Figure 11. Bimodal throat size distribution with standard deviation 2.0 µm for smaller (the first peak) and 2.5 µm for larger (the second peak) throats.

(a)
As can be seen, for the monomodal PN in Figure 12, widespread invasion patterns with a high number of invasions is not achieved. This is in contrast to the bimodal PN in Figure 13. Monte Carlo simulations yielded an average gas saturation of 26% for the monomodal network and 38% for the bimodal network. This shows that widening of the PSD (Figure 9) and the introduction of macro pores (Figure 11), both, result in a change of the invasion process. This change is more significant in the second situation. While capillary fingering with narrow single gas branches is rather favored by narrow PSDs, widening of the invasion front with higher gas saturations is obtained by larger PSDs, and larger porosities. However, it can also be shown that the widening of the front can also be achieved when the porosity is kept constant and only the PSD is adjusted (Figure 13). The monomodal and bimodal networks used in simulations have a constant porosity of 71%. Figures 12 and 13 also show the saturation profiles of these simulations. They reveal the importance of the consideration of the PSD for characterization of the invasion process.

Bimodal Pore Size Distributions
The previous results analyze the influence of the PSD on the example of monomodal PNs, i.e. only one peak in the histograms in Figure 9. As can be seen from Table 3, the porosity is only marginally affected in these cases. Due to this, the phase distribution patterns change only slightly in Figure 10. In reality, pore structures often obey bimodal PSDs, i.e., with two peaks in the histogram ( Figure 11) and with much stronger impact on porosity. The influence of macro pores is highlighted in Figure 3.
As can be seen, for the monomodal PN in Figure 12, widespread invasion patterns with a high number of invasions is not achieved. This is in contrast to the bimodal PN in Figure 13. Monte Carlo simulations yielded an average gas saturation of 26% for the monomodal network and 38% for the bimodal network. This shows that widening of the PSD (Figure 9) and the introduction of macro pores (Figure 11),both, result in a change of the invasion process. This change is more significant in the second situation. While capillary fingering with narrow single gas branches is rather favored by narrow PSDs, widening of the invasion front with higher gas saturations is obtained by larger PSDs, and larger porosities. However, it can also be shown that the widening of the front can also be achieved when the porosity is kept constant and only the PSD is adjusted (Figure 13). The monomodal and bimodal networks used in simulations have a constant porosity of 71%. Figures 12  and 13 also show the saturation profiles of these simulations. They reveal the importance of the consideration of the PSD for characterization of the invasion process. These findings are in very good agreement with literature predictions on drainage invasion regimes [48,56,57]. According to [57], the invasion occurs in the capillary dominated regime, i.e. with rather low injection rate and negligible viscous forces in both fluids. The viscosity ratio of water/air is around 50, and the capillary number is calculated from the ratio of viscous over capillary forces: These findings are in very good agreement with literature predictions on drainage invasion regimes [48,56,57]. According to [57], the invasion occurs in the capillary dominated regime, i.e. with rather low injection rate and negligible viscous forces in both fluids. The viscosity ratio of water/air is around 50, and the capillary number is calculated from the ratio of viscous over capillary forces: These findings are in very good agreement with literature predictions on drainage invasion regimes [48,56,57]. According to [57], the invasion occurs in the capillary dominated regime, i.e., with rather low injection rate and negligible viscous forces in both fluids. The viscosity ratio of water/air is around 50, and the capillary number is calculated from the ratio of viscous over capillary forces: with liquid pressure difference ∆P l and capillary pressure gradient ∆P c . It depends on the viscosity of the liquid phase η, the wetting curvature by means of surface tension σ, and the velocity of the invading menisci v, as And Ref. [58] considering the PSD with mean throat radius r and standard deviation σ 0 . (In Equation (3), L t denotes the length of a throat). The conditions for capillary fingering are fulfilled, if the capillary number is below 10 −4 [48,59]. With the given values for viscosity, surface tension and PSD, this would be achieved if the invasion velocity is below 3.3 mm/min (with liquid viscosity of water η = 10 −3 Pa s and surface tension σ = 0.073 N/m at 20 • C, throat length L t = 50 µm and standard deviation = 1.5 µm). This is in good agreement with correlations of the invasion time and current density estimated in [46]. Transition to another invasion regime, e.g., stable displacement or viscous fingering [48] would occur if the invasion velocity would be increased or also if the standard deviation of the throat size would be decreased. Such effects are observed in the research of drying of PNs, where sufficiently narrow PSDs can result in stable, i.e., viscosity-dominated invasion [58]. When viscosity comes into play, the width of the invasion front scales with capillary number. However, in the absence of viscosity, the probability of the throat invasion depends on the PSD. In the limit of identical throat sizes, all throats would be equally selected for invasion. In the case of our simulation, where in a such a limit the selection would not be stochastically distributed but ordered by the throat number, the invasion would always occur in the throat with the lowest number (also refer to Figure 4) wherefore consequently this special situation (i.e., no distribution of throat and pore sizes) could not be simulated with the current algorithm. In contrast, when the throat size distribution becomes broader, the invasion follows the path of the least resistance, which results in a more random distribution of the phase patterns, provided that the throat sizes are randomly distributed. Following [60] and [61], the occupation probability decreases with growing ∆r.

Pore Network Simulations of Real Porous Structures
The above discussions are referred to rather artificial porous structures aiming at a more fundamental understanding of the influence of the pore structure on the invasion behavior. In the following, we would like to compare the results for realistic pore structures, extracted from felt, foam and sintered PTLs. The simulation results are compared to literature values from Arbabi et al. [20]. For this purpose, the PN simulations were repeated with a single-entry point (similar as in Figure 7). In more detail, the invasion starts from a single gas pore while all other surface pores are not connected to the oxygen supply. Simulations for each type of PTL were repeated five times with the PSDs shown in Figure 14. From these simulations, representative patterns, i.e., which were most frequently observed, are shown in Figure 15. It is seen that felt and foam PTL show more constricted patterns, while the sintered PTL allows relatively broad patterns with more gas invasions, which is explained well enough by the assumed PSDs for each type (Figure 14).
It is remarked that the PSDs are usually not provided in literature [20]. Due to this, we have selected a standard deviation of throat width so as to catch the gas-liquid phase distributions found for microfluidic experiments in the literature. The agreement of the simulated and the experimental invasion behavior is very well. Also, the results are in line with the discussion of the impact of PSDs above. Namely, the foam exhibits one capillary finger that invades the PN almost straightly. The felt PTL, with a slightly broader PSD (Figure 15), instead reveals a slight widening of the invasion front in horizontal direction. In contrast to that, the sintered PTL allows for a broadening of the invasion front and a significantly higher gas saturation at breakthrough.
In addition to that, it becomes clear from Figure 15 that the effect of liquid clustering occurs in all types of the PTL. The greatest number of isolated liquid clusters is observed in the sintered PTL where the number of invaded pores is higher and the invasion front appears more ramified. In the case of felt and foam PTL, isolated clusters are not seen on the experimental images. This could be because of the small size of the PN and the overall lower number of pore invasions. It is remarked that the PSDs are usually not provided in literature [20]. Due to this, we have selected a standard deviation of throat width so as to catch the gas-liquid phase distributions found for microfluidic experiments in the literature. The agreement of the simulated and the experimental invasion behavior is very well. Also, the results are in line with the discussion of the impact of PSDs above. Namely, the foam exhibits one capillary finger that invades the PN almost straightly. The felt PTL, with a slightly broader PSD (Figure 15), instead reveals a slight widening of the invasion front in horizontal direction. In contrast to that, the sintered PTL allows for a broadening of the invasion front and a significantly higher gas saturation at breakthrough.
In addition to that, it becomes clear from Figure 15 that the effect of liquid clustering occurs in all types of the PTL. The greatest number of isolated liquid clusters is observed in the sintered PTL where the number of invaded pores is higher and the invasion front appears more ramified. In the case of felt and foam PTL, isolated clusters are not seen on the experimental images. This could be because of the small size of the PN and the overall lower number of pore invasions.  It is remarked that the PSDs are usually not provided in literature [20]. Due to this, we have selected a standard deviation of throat width so as to catch the gas-liquid phase distributions found for microfluidic experiments in the literature. The agreement of the simulated and the experimental invasion behavior is very well. Also, the results are in line with the discussion of the impact of PSDs above. Namely, the foam exhibits one capillary finger that invades the PN almost straightly. The felt PTL, with a slightly broader PSD (Figure 15), instead reveals a slight widening of the invasion front in horizontal direction. In contrast to that, the sintered PTL allows for a broadening of the invasion front and a significantly higher gas saturation at breakthrough.
In addition to that, it becomes clear from Figure 15 that the effect of liquid clustering occurs in all types of the PTL. The greatest number of isolated liquid clusters is observed in the sintered PTL where the number of invaded pores is higher and the invasion front appears more ramified. In the case of felt and foam PTL, isolated clusters are not seen on the experimental images. This could be because of the small size of the PN and the overall lower number of pore invasions.

Summary and Conclusions
In this study, the applicability of PMNs for the simulation of water drainage from PTL was discussed. The PNM applies invasion percolation rules for hydrophilic drainage. It was shown that the PNM can reliably predict the invasion patterns of 2D microfluidic devices related to water electrolysis studies. In addition to that, the structure dependence of the invasion process was studied using various types of PTLs (foam, felt, sintered). In agreement with experimental data from literature, we could predict different invasion patterns for the three cases. Felt and foam PTLs showed narrow gas fingers while sintered PTL showed widespread invasion front patterns with a higher number of trapped clusters. This observation matches the theoretical investigations based on a Monte Carlo study of PSDs and invasion probability. Based on this, it could furthermore be shown, that the PSD affects the invasion patterns more significantly than the porosity. More clearly, it was revealed that the PSD could be adjusted independently of the porosity and that this resulted in different invasion patterns. The impact of PSD can also be extended towards porous media in fuel cells.
As a next step, the PN will be further enhanced to replicate the local coordination numbers in real PTLs with non-uniform distribution of pores and throats. The purpose would be to simulate mass transfer within the real structure of porous media rather than in an idealized network. These simulations would then be much closer to reality compared to the ones with a fixed coordination number in the network. The effective transport parameters, e.g., permeability, could also then be extracted by solving mass transfer equations pore by pore in a real structure. It would also be important to study the effect of local temperature changes in the system and liquid flow through corner films. Unsteady changes in the current density could also lead to pressure changes. For this reason, the application of imbibition rules along with drainage will become important.

Conflicts of Interest:
The authors declare no conflict of interest.