Computational Optimization of Porous Structures for Electrochemical Processes

Abstract

Porous structures are naturally involved in electrochemical processes. The specific architectures of the available porous materials, as well as their physical properties, crucially affect their applications, e.g., their use in fuel cells, batteries, or electrolysers. A key point is the correlation of transport properties (mass, heat, and charges) in the spatially—and in certain cases also temporally—distributed pore structure. In this paper, we use mathematical modeling to investigate the impact of the pore structure on the distribution of wetting and non-wetting phases in porous transport layers used in water electrolysis. We present and discuss the potential of pore network models and an upscaling strategy for the simulation of the saturation of the pore space with liquid and gas, as well as the computation of the relative permeabilities and oxygen dissolution and diffusion. It is studied how a change of structure, i.e., the spatial grading of the pore size distribution and porosity, change the transport properties. Several situations are investigated, including a vertical gradient ranging from small to large pore sizes and vice versa, as well as a dual-porosity network. The simulation results indicate that the specific porous structure has a significant impact on the spatial distribution of species and their respective relative permeabilities. In more detail, it is found that the continuous increase of pore sizes from the catalyst layer side towards the water inlet interface yields the best transport properties among the investigated pore networks. This outcome could be useful for the development of grading strategies, specifically for material optimization for improved transport kinetics in water electrolyser applications and for electrochemical processes in general.

1. Introduction

The control of process efficiency is in most cases based on the integral measurement of relevant parameters, such as the temperature, pressure, overall mass, concentration, and current density. Spatially and temporally resolved data are not available in most situations, either because of a lack of advanced equipment with high spatial and temporal resolution or because of data inaccessibility in the process. To achieve the required data accessibility in order to gain in situ information, it is often necessary to adjust the process for use with model experiments, which can result in the process deviating from reality. Furthermore, in most cases, having high spatial resolution compromises the temporal resolution. This especially restricts the research on structure interactions of porous media with overall process behavior. Examples include the application of electrodes in batteries, fuel cells, and water electrolysis [1,2,3,4,5]. Electrodes are composed of catalyst layers, gas diffusion or porous transport layers, as well as ion exchange membranes. Independent of the application, the involved porous structures are usually characterized by the ratio of the void fraction to the solid content, the average pore size and pore size distribution, the pore connectivity, the permeability, the effective diffusion coefficients, the effective thermal conductivity, and the wettability, etc. The structural material parameters are typically obtained from experimental measurements such as mercury porosimetry, BET (Brunauer-Emmett-Teller) measurement, scanning electron microscopy (SEM), microcomputed X-ray tomography (µ-CT), and even ultrasonic measurement and pulsed thermography on a larger scale [6,7,8,9,10]. From these methods, only microscopy and tomography are able to resolve the pore structure locally [11,12]. The other methods can only provide global information, such as information on the overall porosity, average pore size, and pore size distribution. In addition, transport parameters such as the relative permeability, diffusivity, thermal conductivity, and wettability are very often also determined as effective parameters that are averaged at the macro scale [13,14]. In particular situations, detailed information is required [15,16,17,18,19]. For example, transport losses in fuel cells or water electrolysis are mainly associated with the mass transfer resistances within the porous electrodes, which strongly depend on the local structure or wettability [20,21,22,23]. Various other examples are found in other research fields, such as hydrology or process engineering [24,25]. In many cases, several different substances are transported through the porous medium at the same time, and very often two-phase flow is a major issue. Especially in these situations, the distribution of species between larger and smaller pores is very relevant. While the liquid (wetting phase) and gas (non-wetting phase) are rather stagnant inside the smaller pores, pathways can develop through interconnected larger pores (Figure 1). Again, independent of the application, it is of interest to clearly identify how pores of different sizes are interconnected and how the overall architecture of the porous medium can be broken down into the local pore structure.

It is possible to analyze the impact of the distributed pore space on mass (and also heat) transfer using pore-scale modeling. Pore network models have been successfully applied to simulate water invasion in GDL microstructures used in proton exchange membrane fuel cells [26] and to simulate reactive transport [25].

In pore network modeling, two major routes are usually taken. (i) In the first option, the complex problem is simplified to single capillary experiments and simulations, as well as regular lattices [27,28,29]. More recently, with the development of more precise measurement techniques and more powerful computational tools, a second option is gaining more attention [18,30,31,32,33]. (ii) In this option, the structure of the porous medium is scanned using X-ray tomography, neutron tomography, or more recently FIB-TEM (Focused Ion Beam-Transmission Electron Microscopy) or FIB-SEM (Focused Ion Beam-Scanning Electron Microscopy) [34]. The images are then prepared for use for mass transfer simulations. Route (ii) certainly has an advantage in that the processes inside the porous medium can be visualized more realistically. However, the reliability of such models always depends on assumptions concerning the application of image filters, mesh precision, and applied thresholds in the image analysis and meshing of the computational domain. In this sense, route (i) allows more accurate physical correlation between the idealized structure and invasion.

In this paper, we sketch a general route for the estimation of transport properties based on pore network models. The application of pore network models can be useful in clarifying open questions related to the transport properties of differently structured materials, which are either estimated experimentally [14,35,36] or by macroscopic modeling [37]. Some examples are the more precise correlation of the structure and transport properties, the study of the impact of liquid cluster formation, and the role of the wetting of liquid films. As a first step, we use idealized 2D and 3D regular lattices, firstly to draw an overall picture of the impact of the pore structure on the distribution of wetting and non-wetting phases during a drainage invasion process, and secondly to develop a generic upscaling strategy using the transport coefficients and morphological data from the pore network model in a macroscopic continuum model. While we have demonstrated the relevance of pore size distribution in a recent paper related to water electrolysis [38], here the focus is on the defined gradients of the pore structure. The background for this investigation involves two factors. On the one hand, porous materials applied in electrochemical processes are very often produced with spatial gradients in the pore size [39]. In a wide range of production methods, the structure development occurs rather randomly under uncontrolled conditions [40]. Secondly, the softer materials, such as felt or foam, are very often compressed inside the technical equipment, resulting in a gradient of pore sizes [41,42]. On the other hand, it is of interest to assess how controlled grading of the pore structure can influentially impact the improvement of the overall transport kinetics of porous materials. This can be useful for classical electrochemical processes [43]. It is, therefore, important to better understand the correlation between structure and transport properties.

We introduce here a general approach that might be adapted for several situations. As an example, we use invasion of a non-wetting phase (oxygen) into a pore network saturated with a wetting phase (water) under hydrophilic conditions, as found in water electrolysis [44], and discuss the impact of the pore structure on the specific permeabilities of the partially saturated pore network. The effective parameters are applied in a continuous model to show the dynamic spatial distribution of oxygen inside the porous transport layer (PTL).

The paper is organized as follows. At first, the pore network modeling based on regular pore and throat lattices is briefly recalled in Section 2, with a focus on displacement processes occurring in the PTL during water electrolysis. In Section 3, the simulation results of the invasion of the non-wetting phase in graded pore networks are presented and discussed. They are further exploited for parameterization of a 1D macroscopic continuum model in Section 4. For this purpose, the morphological parameters obtained from the 2D pore network simulations in Section 3 will be merged into a 1D property coordinate.

2. Materials and Methods

2.1. Pore Network Model

The pore network model used in this study was developed in [45] and validated by comparison with reference experiments from the literature in [38]. The model is comprised of two tools: (i) the structure computation tool and (ii) the invasion simulation tool. At first, the structure is described extensively to give a detailed overview of the studied gradients. The invasion part is shortly summarized below—more details are given in [38,45].

Figure 2 illustrates a very frequent example of electrochemical applications. It shows compressed titanium felt with a graded pore structure, as used in water electrolysis (Figure 2a). The mean pore sizes and porosities of the titanium felt were estimated using GeoDict (Math2Market GmbH, Kaiserslautern, Germany) (Figure 2b). As can be seen, the pore structure differs between the top and the bottom sides. More clearly, the pore size increases from the catalyst layer (applied at the top side) to the water supply channel (connected to the bottom side). The variation is mainly induced by the manufacturer and is partly a result of compression of the felt (by 10–15%) during operation.

We use the example in Figure 2 as a starting point for the simulations. In more detail, we computed several pore networks with vertical gradients in pore size and porosity (note that the porosities of the regular pore network considered here are small compared to irregular porous media). In addition to the low-to-high (LTH) gradient shown in Figure 2, we also study the impact of a high-to-low (HTL) gradient. For the LTH gradient, two cases are studied. In the first case, a constant decrease in the average pore size from top to bottom is realized (Figure 3); in the second case, an abrupt change in porosity similarly as in Figure 2 is imposed on the pore structure.

2.2. Structure Tool

The mean pore and throat sizes $2\cdot \overline{r}=\overline{d}$ used in the pore network simulations are based on the example shown in Figure 2. In more detail, we used cylindrical throats with diameter $d_{ij}$ (or radius $r_{ij}$) and length L_ij = 50 μm = const and applied Gaussian-distributed throat sizes in the range of $\overline{r}\pm std\left(r\right)$ (

Table 1. Overview of simulation parameters.

	2D	3D
Network size (columns and rows)	80 $\times$ 10	10 $\times$ 10 $\times$ 7
Number of pore network slices	8	5
Number of pores	800	700
Number of throats	1360	1600
Average pore diameter range	(30 to 33.5) µm ± 4 µm	(30 to 32) µm ± 4 µm
Average throat diameter range	(22 to 29) µm ± 6 µm	(22 to 26) µm ± 6 µm
Lattice spacing	50 µm	50 µm
Thickness of network	500 µm	350 µm
Maximum porosity range	45–65%	30–50%

). For the pores, we used a spherical geometry with the mean and standard deviation given in Table 1.

Besides the constant gradient of throat and pore sizes imposed on $\overline{r}$, with $d\overline{r}/dz$ = $\pm$0.01 µm/µm for throat sizes and $d\overline{r}/dz$ = $\pm$0.005 µm/µm for pore sizes (Figure 3), a dual-porosity network was realized (Figure 4 and Figure 5). Both examples result in overlapping of the throat and pore sizes of each layer, as illustrated in Figure 4 and Figure 5.

In summary, the situations studied in this work are: (1) a reference 2D network without a gradient; (2) a 2D network with a continuous LTH gradient; (3) a 3D network with a gradient identical to case 2; (4) a 3D network with the same gradient as in case 2 but in the opposite direction (i.e., HTL); (5) a 3D network with two different porosity regions (denoted as a dual-porosity pore network, similarly to Figure 2). The spatial distribution of the porosity for all cases is provided in Figure 5. As can be seen, the variation of pore sizes has a direct impact on porosity. This is explained by the constant lattice spacing, i.e., the distance between the pores, as well as the constant overall network height (width and depth in 3D).

The simulation parameters are summarized in Table 1. The pore and throat sizes, as well as the porosities, indicate the selective parameters of the investigated gradients, whereas the wetting angle ($\theta$ = 60°) and temperature ($T$ = 80 °C) are kept constant. Note that the wetting angle was chosen so as to imitate the situation of water invasion inside titanium felt. The pressure is successively increased from 1 bar to $P_g\left(r_{ij}\right)$ in order to enable invasion of the pore network.

2.3. Invasion Tool

The pore network simulation is based on the drainage invasion algorithm with liquid clustering presented in [38,45]. The pores and throats are initially completely filled with liquid water. Oxygen then invades from the top side by stepwise increase of the gas pressure above the critical invasion pressure thresholds of the liquid filled pores. We assume high current densities and the formation of gas fingers with plug flow of the gas phase [46]. The order of water drainage from the throats and pores is computed from the capillary pressure curve based on the Young–Laplace equation:

$$P_l=P_g-\frac{2\sigma \cos \theta }{r_t}$$

(1)

where $P_l$ is the liquid pressure, $P_g$ is the gas pressure, $\sigma$ = 0.0627 kg s⁻² is the surface tension of water, $\theta$ is the wetting angle, and $r_t$ is the throat radius of throat $ij$. The pore invasion is computed accordingly in competition with the throat invasion.

The invasion process is assumed as quasi-steady, capillary dominated, non-viscous flow. It is assumed that steady gas–liquid distributions are obtained at breakthrough of the gas phase to the bottom of the pore network and at disconnection of the water supply route from the bottom side to the top by further increase of the gas pressure. The model assumptions (i.e., stepwise invasion, steady-state invasion patterns) are made in accordance with experimental results reported in [46] for low capillary numbers and incompressible gas flow, as well as on the simulation results for gas permeation in [47].

Once the steady invasion profiles are obtained, gas and liquid flow rates through the partially drained pore networks are computed from the Hagen–Poiseuille equation:

$${\dot{M}}_{l,g}=\frac{{\rho }_{l,g}\pi r_t^4}{8{\eta }_{l,g}L}\left(P_{l,g,1}-P_{l,g,2}\right),$$

(2)

Equation (2) takes into account that the transport paths for gas and liquid are hindered by isolated liquid clusters. The subscripts 1 and 2 in Equation (2) are the two pore neighbors of a throat $ij$ between which the pressure gradient is considered for the computation of flow. For this, a wetting liquid and an incompressible, ideal, isothermal gas are assumed as a first step. The discrete solution of the Hagen–Poiseuille equation can be used for the computation of the relative permeabilities of gas and liquid phases at the Darcy scale, assuming laminar flow:

$$K\cdot k_{rel}=\frac{{\dot{M}}_{l,g}}{A_{l,g}}\frac{\Delta z}{\Delta P_{l,g}}.$$

(3)

The value of absolute permeability $K$ is found for the same computations but with a completely saturated pore network with gas or liquid [45].

3. Results and Discussion

3.1. General Results

The 2D simulation results can only provide an orientation for the basic invasion mechanisms. In 3D, the pore connectivity is higher and different permeabilities are computed. However, the invasion profiles are very similar in 2D and 3D simulations. We, therefore, start here with the analysis of the 2D profiles before discussion of the 3D pore networks.

The simulated invasion patterns are very sensitive to the throat size variation in the direction of the invasion front, independent of being 2D or 3D. It has already been found that a change in the gradient can change the resulting invasion patterns. For example, gradients of $d\overline{r}/dz$ = −2e⁻³ did result in invasion patterns similar to invasion without a gradient, whereas for $d\overline{r}/dz$ = −0.02 an overlapping of throat sizes occurred and had a severe impact on the invasion patterns. This is illustrated in Figure 6 in the example of 2D simulations with HTL gradients; a similar finding was also made for 3D simulations using the same gradients. In Figure 6, the Gaussian-distributed throat sizes of selected layers of the 2D pore network are shown together with the invasion patterns. The invasion behavior changes from capillary fingering in Figure 6a to a flat front in Figure 6b. The simulations were repeated 15 times; the given invasion patterns are representative patterns. The example in Figure 6 clearly illustrates the relevance of the height of the gradient for the final saturation with gas and liquid, as well as the structure of the gas–liquid distribution.

Figure 7 summarizes the invasion profiles for 4 different gradients at breakthrough of the gas phase for both 2D and 3D simulations based on the parameters given in the previous section. In the LTH case, the larger pores are located at the bottom side. In this situation, invasion is favored towards the bottom of the network because of the lower invasion pressure thresholds that can be computed from Equation (1). Compared to the HTL simulation, the LTH is characterized by an invasion pattern such as viscous fingering [48] in both 2D and 3D simulations. Breakthrough of the gas phase occurs with a significantly lower number of invasions steps. This means that the final saturation with the gas phase at the breakthrough point is lower than compared to the HTL simulation, where the larger voids are located at the top side. In contrast to LTH, in HTL the top side is favorably invaded. Thus, in the case of HTL, the invasion results in the formation of a gas–liquid distribution, similarly as for stable displacement with a flat invasion front [48]. This is again observed in 2D and 3D simulations. However, the stabilization is more significant in 2D due to the lower pore connectivity.

3.2. Drainage Simulation

In order to obtain further insights into the impact of the gradient on the final distributions of gas and liquid and the resulting transport properties, a structure without a gradient is compared to a structure with an LTH gradient to show the effect of varying porosity. Besides the phase distributions in Figure 8, the relative gas and liquid permeabilities and saturation profiles are compared in Figure 9 and Figure 10.

The striking differences between the two simulations are the final saturations with liquid and gas phases at breakthrough, as well as the invasion behavior after breakthrough, if it is assumed that the gas phase can be forced to further invade the porous structure by increasing the gas pressure, e.g., provoked by the very low gas permeability in the studied case of laminar, incompressible gas flow (Figure 10) or based on diffusion resistances [49]. The pore network simulations indicate that the gas saturation at the top side can be reduced from around 50% to around 30% by application of a pore size gradient (Figure 9).

In the pore network without a gradient, the final liquid distribution is almost uniform across the height of the network (Figure 9; the low saturation at the bottom is explained by invasion rules and determination of liquid cluster labels [45]). In contrast, in the LTH pore network, a distinct watershed evolves in the center of the pore network, separating the gas and liquid phases. This is explained by the throat size distribution and the favorable invasion of larger throats in hydrophilic drainage.

The overall liquid saturation at breakthrough is 83% inside the pore network without a gradient and 93% in the LTH network, based on the averages from 15 simulations. The final saturation with liquid is on average 35% in the pore network without a gradient and 41% in the LTH network. Thus, more liquid is contained inside the pore network with the gradient. A similar finding is reported in [46] for the application of a microporous layer. Figure 8a,c and Figure 9 reveal that the LTH gradient is very favorable for viscous fingering. As a consequence, the liquid phase has higher connectivity and less isolated liquid clusters evolve in the simulation with LTH. As shown in Figure 10, the higher liquid connectivity results in an overall higher relative liquid permeability in the pore network with the gradient (between $S$ = 0.8 and $S$ = 0.5). The absolute permeabilities computed for the empty networks are $K$ = 1.65 × 10⁻¹² m² (network without gradient) and $K$ = 1.04 × 10⁻¹² m² (LTH). In [50], similar absolute permeabilities in the range of 10⁻¹² m² to 10⁻¹⁰ m² are reported for non-graded, fiber-based Ti PTLs. In detail, for pore sizes in the range of 32.4 ± 0.6 µm and porosity of 75%, the computed permeability of the non-graded structure is around 17–20 × 10⁻¹² m² in [50].

The two counter-current effects of liquid and gas transport induce an optimization problem in the structure in order to satisfy the optimal transport properties for both phases. From Figure 10, it can be concluded that the optimum gas–liquid configuration to obtain high liquid permeability is between the breakthrough point and S = 0.5 in case of the LTH gradient. In this range the computed gas permeability is also comparably high for the LTH gradient. This is explained by the saturation profiles, as discussed above.

Whether this issue depends on the gradient only or also on the dimension of the pore network (and thus pore connectivity) will be analyzed further using the 3D simulation results in the next section.

3.3. Drainage Simulation

The simulation results obtained from 3D pore networks are very similar to the 2D simulation results. As shown in Figure 11 for LTH and HTL gradients, higher overall gas saturations are obtained for the HTL gradient at breakthrough, whereas in case of the LTH gradient distinct gas fingers penetrate the pore network from top to bottom. When the gas pressure is further increased after breakthrough of the gas phase, a water–gas barrier is formed in the 3D pore networks.

Repetition of the simulations reveals the significance of the invasion behavior. The results are summarized in Figure 12. As can be seen, the liquid saturation at breakthrough of the gas phase is almost uniform in the LTH case (Figure 12a). The highest gas saturation is at the top side, which is around 20%. This is much lower than in the 2D case, where it is around 50%. Even when the simulation is continued to further invade the pore network until disruption of the water connectivity between the top and bottom (denoted as water disconnection), the gas saturation at the top side remains low (at around 40%), while gas invasion occurs mainly at the bottom side (where the liquid saturation is lower). The behavior is opposed in the case of the HTL, where the top side always has a higher gas saturation (and lower liquid saturation, Figure 12b) than the bottom side, where the pores and throats are smaller and the resulting invasion pressures are higher. However, interestingly, the LTH and HTL invasion profiles are not perfectly inversed. Instead, significant differences occur, especially when regarding the breakthrough point, while after breakthrough the saturation profiles tend to become almost opposed (black lines in Figure 12). As shown in Figure 13 for the example of the pore–throat configuration, this strongly affects the permeabilities of the two pore networks. Overall higher relative liquid permeabilities are also found in the LTH network after breakthrough. At the same time, the relative gas permeability is also comparably higher at higher liquid saturation levels. Note that the computed gas permeabilities are generally at a low level due to the assumptions of incompressible ideal gas and slow, laminar flow. Higher values are expected for higher overall gas pressures. In summary, this indicates that the LTH pore network generally has better transport properties than the HTL network.

In Figure 14 and Figure 15, the invasion of a dual-porosity pore network is studied, with similar structure as the titanium felt in Figure 2. Figure 14a reveals distinct gas fingers form inside the low-porosity region, whereas more pore invasions occur inside the high-porosity region than in the LTH case studied above. This observation is similar to the experimental results with microporous layers in [46].

The liquid saturation inside the top layer is almost similar to the LTH gradient (Figure 12a and Figure 14b). However, in the high-porosity region, significantly more oxygen is contained when the water disconnection is regarded. This results in a higher overall gas saturation and the formation of more isolated single-liquid clusters at breakthrough and at water disconnection. As can be found in Figure 15, the liquid permeabilities of the dual-porosity network are lower than in the LTH situation. The gas permeabilities are significantly lower and remain at a constant low level because the low-porosity region is not further invaded, even if the gas pressure is increased after breakthrough.

If one deduced a grading strategy from this comparison, a pore network with an LTH gradient would be recommended over a dual-porosity pore network if low gas saturations and higher liquid permeabilities were realized. Further studies in this direction incorporating simulations of the pore structure obtained from µ-CT measurements and artificial irregular pore networks with different pore size gradients are already in progress and can be used as a base for the development of a grading strategy in the future. Validation of the outcomes using Ti–felt PTLs could be based on X-ray synchrotron imaging or neutron imaging [49,50,51], as well as on experimental data, such as from impedance spectroscopy [52,53,54]. Apart from the optimization of mass transfer properties, the change of the ohmic resistances for charge transfer must also be considered to achieve the optimal structural design [35]. Experimental studies in the literature indicate that the realization of the connectivity of the PTL with the catalyst layer and the degradation of material properties play major roles in performance enhancement [35,49,55], while other studies reveal the impacts of temperature and pressure [47].

3.4. Upscaling to Continuum Model

For the formulation and parameterization of macroscopic continuum models, such as those used to simulate fuel cells or water electrolysis on a larger scale, experimental data is usually required. Here, we would like to discuss how pore network modeling can alternatively be applied to parameterize macroscopic models. In the following, the diffusive mass transfer $\dot{M}$ though the porous medium shall be regarded. This can generally be described through the correlation of the driving force, which can be a pressure gradient $\Delta P$ over distance $z$, multiplied by a transport coefficient $k$ that is related to a resistance factor $R$:

$${\dot{M}}_{eff}=A\frac{k}{R}\frac{\Delta P}{z}$$

(4)

The resistance factor basically indicates how much of the mass transfer through a porous medium is reduced compared to diffusion across a free surface with an identical cross-section. This depends on the porosity $\epsilon$ of the porous medium:

$$\epsilon =\frac{V_0}{V_0+V_s}$$

(5)

as well as on the tortuosity $\tau$, which can in the simplest form be given as the ratio of the length $l_p$ of a straight pore with an equivalent diameter and the same transport resistance as the porous medium, divided by the thickness of the porous medium $L$:

$$\tau =\frac{l_p}{L},$$

(6)

Or by the ratio of the geodesic length to the Euclidian length [56]. By definition, $0\le \epsilon \le 1$ and $\tau \ge 1$. In more advanced approaches, tortuosity is computed as a function of porosity, i.e., $\tau =f\left(\epsilon \right)$. One example is the Bruggeman relationship [57], where:

$$\tau ={\epsilon }^{1-\alpha }$$

(7)

with Bruggeman exponent $\alpha$. However, the complex dependency of the tortuosity on the pore structure (i.e., the pore size distribution, pore length, and connectivity) makes it almost impossible to identify $\tau$ correctly in most situations. Thus, in most cases $\tau$ remains a fitting parameter. In the cases studied here, $\tau$ = 1, as $l_p=L$ in the pore networks with regular lattices [57]. This means that the resistance factor

$$R=\frac{\tau }{\epsilon }.$$

(8)

becomes simply

$$R=\frac{1}{\epsilon }.$$

(9)

Although more complex situations (with $\tau \ne$ 1) might be studied using pore networks based on structure estimations from microscopy or tomography scans, in a first step our focus here is on porosity gradients.

In a homogeneous porous medium with equally distributed void space across the thickness, $R=1/\epsilon$ can be constant. In such situations, the effective mass flow rate could be related to the mass flow rate ${\dot{M}}_0$ through the empty box as:

$${\dot{M}}_{eff}=\frac{{\dot{M}}_0}{R}$$

(10)

In a regular pore network with $\tau$ = 1, Equation (10) reads as:

$${\dot{M}}_{eff}=\epsilon \cdot {\dot{M}}_0$$

(11)

However, in the situations studied in this paper, the porosity varies across the material. Ergo, the transport resistance $R$ also varies with height. As a consequence, as an alternative route we describe here the possibility to use the results from pore network modeling for the simulation of mass transfer through a porous medium with a pore size gradient instead of experimental estimation of spatially distributed parameters. As an example, we will discuss here the transport of dissolved oxygen in water through a partially liquid-saturated pore network. In a first step, we consider the 2D regular pore networks with the LTH gradient described above. In these networks, $\tau$ = 1, and thus $R$ is defined from Equation (9). The parameter of interest here is, thus, the spatially distributed porosity $\epsilon$. In addition, in the partially saturated pore network, the ratio of the liquid phase also has to be taken into account, as the diffusion of oxygen can only occur through the interconnected water clusters. Both parameters, i.e., the spatially distributed porosity (Figure 5a) and liquid saturation (Figure 9a), were calculated for each slice of the corresponding pore networks and used in the macroscopic model, as discussed in the following section.

3.5. Macroscopic Continuum Model

As already mentioned before, we consider here the spatial and temporal distribution of the concentrations of dissolved species in the liquid phase within the partially saturated pore network. In the case of water electrolysis, the transport of dissolved oxygen through the connected water clusters is of interest. Oxygen is produced inside the catalyst layer on top of the porous transport layer (not regarded in the simulations) and then invades the PTL in the form of gas fingers. Depending on the absorption coefficient ${\alpha }_{O_2-H_{2}O}$ for oxygen in water at the given temperature and pressure, part of the oxygen is dissolved inside the liquid phase. The dissolved oxygen is then transported along the concentration gradient, i.e., counter-currently to the water, by diffusion. The transport direction is, thus, from top to bottom (i.e., against the direction of the z-axis).

In order to simulate oxygen diffusion inside the liquid phase, a lumped one-dimensional macroscopic model is used. The model averages transport coefficients over the cross-section of the network and describes gradients along the z-axis (thickness) of the network. In general, the concentration of a species $A$ in the liquid phase at a specific point in the network z and at time point t is given as $c_A\left(t,z\right)$. The dynamics of this concentration can be characterized by the following partial differential equation

$$\frac{\partial {\varphi }_{\mathrm{l}}\left(\mathrm{z}\right)c_A\left(t,z\right)}{\partial t}+\frac{\partial }{\partial z}\left\{v_l\left(z\right){\varphi }_l\left(z\right)c_A\left(t,z\right)-D_A{\varphi }_l\left(z\right)\frac{\partial c_A\left(t,z\right)}{\partial z}\right\}=R\left(t,z,c_A\right) .$$

(12)

Therein, ${\varphi }_l\left(z\right)$ denotes the fraction of the liquid phase in the network, which can be computed from the porosity and liquid saturation at each position $z$:

$${\varphi }_l\left(z\right)=1-{\varphi }_s\left(z\right)-{\varphi }_g\left(z\right)=ϵ\left(z\right)Sa{t}_l\left(z\right)$$

(13)

which were condensed from the results of the 2D pore network simulation with the LTH pore size gradient at gas breakthrough, as presented in the previous section. Furthermore, the convective velocity of the liquid was obtained from mass conservation of the liquid flow through the network by:

$$v_l\left(z\right)=-\frac{{\dot{M}}_{eff}}{A{\varphi }_l\left(z\right)}$$

(14)

with the total surface area of the network computed from the values given in Table 1. The transport of the species against the water transport direction is described by Fick’s law of diffusion, with diffusion coefficient $D_A$. The source and sinks given by the absorption equilibrium of oxygen in water (which is temperature- and pressure-dependent) along the phase boundary in the vertical direction could principally be accounted for by the right-hand side term $R\left(t,z,c_A\right)$. For simplicity, we have neglected this term in this first approach. This means that oxygen is supplied only at the top side and the oxygen transfer over gas–liquid phase boundaries within the pore network is neglected.

For simulation of the introduced partial differential equation, a finite volume method was applied [58,59]. The resulting set of ordinary differential equations was solved in MATLAB using ode15s [60]. The calculations are based on dimensionless variables. For this purpose, it was assumed that $c_{O_2}\left(t=0,z\right)$ = 0 (initial condition) and $c_{O_2}\left(t,z=z_{max}\right)$ = 1, as well as $c_{O_2}\left(t,z=0\right)$ = 0 (boundary conditions). The results of the simulation are illustrated in Figure 16 for the example of the LTH gradient. As can be seen, a linear concentration profile is obtained, with a fast transition from the initial concentration profile to the steady state.

It can be deduced that the transport of dissolved oxygen in the liquid phase plays a minor role in the situation studied here, where the oxygen production rate is high enough to form gas fingers through the liquid-saturated pore network. In other words, with steady-state liquid and gas distribution, the gradient of dissolved oxygen along the height of the pore network is not high enough to transport a significant amount of additional oxygen by diffusion through the liquid phase. The high-pressure gradients over the liquid phase and the resulting high water flow rates result in significant viscous flow contributions and a low impact of oxygen diffusion. However, it is expected that oxygen diffusion could play a major role if the transport velocity of water counter-current to the oxygen diffusion direction is low.

The presented method of combined pore network simulation and macroscopic modeling is of specific interest in cases of changes of the pore network structure or changes of the liquid and gas distributions during dynamic invasion. It is a generic method that could also be applied for situations where pore and throat sizes vary over time, e.g., as a result of fouling effects or pore blocking resulting from contaminations in the liquid phase [25,61,62,63]. In such situations, a recursive procedure of alternating simulations is proposed.

4. Conclusions

In this paper, we have shown based on pore network simulations how gradients in pore sizes influence the invasion of gas into a liquid-saturated porous medium, as found in a PTL in water electrolysis. The different situations that were studied are gradients with increasing pore size in the invasion direction of the gas phase (LTH), the opposing structure (HTL), as well as a dual-porosity pore network. The simulation results were discussed for 2D and 3D pore networks. The results were compared to a pore network without a pore size gradient in a 2D simulation. It was shown that the final saturation with the gas phase decisively depends on the direction of the gradient, as well as on the design of the gradient, as the dual-porosity pore network revealed a different invasion behavior than the LTH network. Higher overall gas saturations at breakthrough of the gas phase were obtained in the case without a gradient, as well as in the pore networks with the HTL gradient, in both 2D and 3D simulations. Only when the gas pressure was increased beyond the breakthrough point did the formation of a water–gas barrier result in higher gas saturation in case of the LTH gradient. Interestingly, the obtained saturation profiles for LTH and HTL were not opposed when the breakthrough point was considered. The HTL gradient revealed a significantly lower slice saturation with liquid at the top side as compared to the results obtained for the LTH case. In the dual-porosity pore network, a low gas saturation (high liquid saturation) was observed at the low-porosity side. However, inside the high-porosity region, gas invasion resulted in the formation of more isolated liquid clusters, and thus in lower overall liquid permeability than in the LTH network, leaving some open questions for future work.

In summary, it could be demonstrated that optimization of transport properties in both phases (i.e., with high-gas and high-liquid permeabilities) is a bilateral problem, because higher gas permeabilities always result in the reduction of liquid permeability. This is especially problematic if the liquid phase is split up into a high number of isolated liquid clusters that hinder the transport of liquid and gas. This was mainly observed in the simulations without a pore size gradient, as well as in simulations with the HTL gradient and the dual-porosity pore network. In contrast, the simulation results indicate the overall better transport properties of the LTH gradient, where distinct gas fingers provide pathways for the gas phase, while at the same time maintaining high liquid permeabilities. Further systematic studies of structure optimization strategies also based on µ-CT images of different PTLs are already in progress.

Secondly, we have introduced an upscaling strategy based on combined pore network and macroscopic modeling. We showed how the distributed 2D morphological parameters porosity and liquid saturation obtained from the discrete pore network model can be lumped into a 1D macroscopic model. As an example, we discussed here the relevance of the transport of dissolved oxygen through the partially liquid-saturated pore network with the assumption of a high current density, and thus higher water flow rates. It was shown that in case of the assumption of steady-state liquid and gas distributions obtained at breakthrough of the gas phase, a constant linear oxygen profile was obtained within a few seconds. The diffusion flow rates were negligible compared to the oxygen flow through the gas fingers. However, they are expected to become much more relevant in cases of low current densities and in the transition of the transport regime from gas fingering towards bubble flow. This example reveals the possibility to couple pore network modeling and macroscopic modeling in situations where the structural effective parameters and measurement data are not available from experiments.

Symbols and Abbreviations

Symbols
Abase	Base area	m2
$c$	Normalized concentration	-
$c_A$	Concentration	kg/m3
$d$	Diameter	µm
$D$	Diffusion coefficient	m2/s
$ij$	Throat index	-
$k$	Relative permeability	-
$K$	Absolute permeability	m2
$L$	Distance between nodes	µm
$L_t$	Length of throat	µm
$\dot{M}$	Mass flow rate	kg/s
$P_c$	Capillary pressure	Pa
$P_l$	Liquid pressure	Pa
$r_p$	Pore radius	µm
$r_t$	Throat radius	µm
$\overline{r}$	Mean radius	µm
$R$	Resistance factor	kg/m3∙s
$t$	Time	s
$v_l$	Velocity	m/s
$V_p$	Volume of pore	m3
$V_s$	Volume of solid	m3
$V_t$	Volume of throat	m3
$V_0$	Volume of voids	m3
$z$	Thickness of PTL	m
$\epsilon$	Porosity	%
$\eta$	Kinematic viscosity	m2/s
${\rho }_l$	Density of liquid	kg/m3
$\theta$	Contact angle
$\sigma$	Surface tension	N/m
$\varphi$	Liquid fraction	-
$\tau$	Tortuosity	-
Abbreviations
µ-CT	Microcomputed tomography
FIB	Focused ion beam
HTL	High-to-low
LTH	Low-to-high
PSD	Pore size distribution
PTL	Porous transport layer
SEM	Scanning electron microscopy
TEM	Transmission electron microscopy