**Electrolysis Processes**

Special Issue Editor **Tanja Vidakovic-Koch**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Special Issue Editor* Tanja Vidakovic-Koch Max Planck Institute for Dynamics of Complex Technical Systems Sandtorstraße Germany

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Processes* (ISSN 2227-9717) (available at: https://www.mdpi.com/journal/processes/special issues/ electrolysis processes).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Article Number*, Page Range.

**ISBN 978-3-03936-386-5 (Pbk) ISBN 978-3-03936-387-2 (PDF)**

c 2020 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**



### **About the Special Issue Editor**

**Tanja Vidakovic-Koch** studied chemical engineering at University of Belgrade, Serbia, where she also obtained her Magister of Science degree in electrochemistry. She pursued further doctoral studies in process and systems engineering at Max Planck Institute for Dynamics of Complex Technical Systems (MPI-DCTS), Magdeburg, Germany, where she obtained her Dr.-Ing. degree as well as lecturer (Habilitation) qualification for Electrochemical Process Engineering ("venia legendi") from Otto von Guericke University (OvGU) in Magdeburg, Germany. Currently she is Head of the Electrochemical Energy Conversion group at the MPI-DCTS, Magdeburg, Germany, and she lectures at OvGU. Her research focuses on the model-based analysis of linear and nonlinear dynamics of electrochemical energy conversion in polymer electrolyte membrane (PEM) fuel cells, PEM electrolyzers, and enzymatic fuel cells.

### *Editorial* **Editorial on Special Issue Electrolysis Processes**

### **Tanja Vidakovi´c-Koch**

Electrochemical Energy Conversion, Max Planck Institute for Dynamics of Complex Technical Systems, D-39106 Magdeburg, Germany; vidakovic@mpi-magdeburg.mpg.de

Received: 7 May 2020; Accepted: 7 May 2020; Published: 14 May 2020

Renewable energies such as solar, hydro or wind power are in principal abundant but subjected to strong fluctuations. Therefore, development of new technologies for storage of these renewable energies is of special interest. Electrochemical technologies are ideal candidates for the use of excess current, and consequently an increased electrification of chemical processes is expected. In this respect, there are different pathways to utilize excess current electrochemically. Intermediate energy storages in (a) chemical energy carriers like hydrogen via water electrolysis or (b) electrochemical energy storage devices like batteries are perhaps most accepted and discussed solutions. Additionally, excess current can be used with the main goal not to be stored for later use, but to solve some environmental issue or for construction purposes. Possible applications are waste water treatment and electromachining. The article collection in this special issue consists of one review paper and nine original research papers and discusses these topics in more detail. As a Guest Editor of this special issue, I thank all authors for their contributions and wish all readers interesting insights and new inspirations for their works.

### **Electrolysis Processes for Intermediate Energy Storage in Chemical Energy Carriers Like Hydrogen–Water Electrolysis**

There are three main technologies for water electrolysis: alkaline water electrolysis (AEL), proton exchange membrane (or polymer electrolyte membrane) electrolysis (PEMEL), and solid oxide electrolysis (SOEL), with two of them (AEL and PEMEL) at high technical readiness level. Despite these facts and intensive discussions on water electrolysis as a key technology for generation of pure hydrogen using renewable electricity, currently less than 4% of hydrogen originates from electrolysis, with the main part originating not from water but from from chlor-alkali electrolysis where hydrogen is a by-product of chlorine production [1]. Broad introduction of cost competitive and preferably zero carbon routes for hydrogen production is urgently required. In the review paper of Brauns and Turek [1], some of the main challenges hindering broader penetration of water electrolysis are discussed with a focus on AEL. As the authors write, AEL is a key technology for large-scale hydrogen production powered by renewable energy. However, conventional electrolyzers are designed for operation at fixed process conditions, therefore, the implementation of fluctuating and highly intermittent renewable energy is challenging. Their system analysis enabled important insights and a roadmap for more energy efficient systems. According to these authors, in order to be competitive with the conventional path based on fossil energy sources, each component of a hydrogen energy system needs to be optimized to increase the operation time and system efficiency. They stress that by combining AEL with hydrogen storage tanks and fuel cells, power grid stabilization can be achieved. As a consequence, the conventional spinning reserve can be reduced, which additionally lowers the carbon dioxide emissions.

Water electrolysis was also in the focus of Vorhauer et al. [2] and Altaf et al. [3]. One of the reasons for lowering energy efficiency of water electrolysis at high current densities is mass transport limitation in porous transport layers (PTL) of water electrolyzers. This issue is intrinsic for both low temperature water technologies (AEL and PEMEL) and is likely caused by the counter current transport of oxygen gas produced at the anode to the educt (water or OH- ions). In two publications by Vorhauer et al. [2] and Altaf et al. [3], this issue was studied with the help of porous network theory and

with special emphasis on PEMEL. Pore network models (PNM) are powerful tools to simulate invasion and transport processes in different porous media with applications across different disciplines like geology, chemical engineering as well as electrochemical engineering (e.g., fuel cell applications). In their pioneering contribution, Vorhauer et al. [2] described a first application of a PNM of drainage for the prediction of the oxygen and water invasion process inside the anodic PTL at high current densities. According to the authors, in the simulation with narrow pore size distribution, the volumetric ratio of the liquid transporting clusters connected between the catalyst layer and the water supply channel is only around 3% of the total liquid volume contained inside the pore network at the moment when the water supply route through the pore network is interrupted; whereas around 40% of the volume is occupied by the continuous gas phase. The majority of liquid clusters are disconnected from the water supply routes through the pore network if liquid films along the walls of the porous transport layer are disregarded. Moreover, these clusters hinder the countercurrent oxygen transport. They also based on these sketches a new route for the extraction of transport parameters from Monte Carlo simulations, incorporating pore scale flow computations and Darcy flow.

In the publication by Altaf et al. [3], results from Monte Carlo pore network simulations are shown and compared qualitatively to microfluidic experiments from literature. The literature-based experimental invasion patterns of different types of PTLs, such as felt, foam or sintered, have shown good agreement with pore network simulations. Additionally, the impact of pore size distribution on the phase patterns of oxygen and water inside the pore network was studied. These very promising results further supported pore network modeling as a valuable tool for gaining new insights in the transport processes in porous layers during water electrolysis.

Reliable models of gas-liquid equilibria in aqueous electrolytes are of significant importance for proper description of many electrochemical processes including gases as products or reactants (major examples are water, brine, CO2 electrolysis, but also reactions in polymer electrolyte fuel cells). Nabipour et al. [4] proposed a novel solubility estimation tool for gases in aqueous electrolyte solutions based on an extreme learning machine (ELM) algorithm. Although the presented study cases are not directly relevant for water electrolysis, the here developed novel methodology has a potential to be transferred to more relevant examples for electrochemical applications.

### **Electrolysis Processes for Intermediate Energy Storage in Electrochemical Energy Storage Devices**

Different types of electrochemical energy storage devices are discussed in the framework of intermediate energy storage from renewables. Common examples are rechargeable batteries like Li-ion or redox flow batteries. In this special issue, two less discussed options are described, a so-called acid-base flow battery [5] and supercapacitors [6]. An acid-base flow battery is proposed by Xia et al. [5]. This very interesting solution is based on reverse electrodialalysis with bi-polar membranes. During charging, the system operates in electrolysis mode; hydrogen and oxygen evolution reactions take place at the cathode and anode, respectively; and acid and base solutions are regenerated in corresponding compartments. Alike normal water electrolysis, the two half-reactions take place at different pH values (acidic and alkaline conditions). During discharge, neutralization of acid and base produces electricity in the process of reverse electrodialysis with bipolar membranes in an apparently fairly overlooked flow battery concept. The authors demonstrated this concept with stack experiments, consisting of 5 to 20 repeating cell units at lab scale. The first results are very promising and the studied system shows already energy density in the similar range of redox flow batteries. The challenges and measures to further increase energy efficiency of this new type of acid-base flow battery are discussed.

Due to the high power fluctuations that are inherent to renewable energy sources, the dynamics of the storage media is of great importance when designing storage concepts for renewable energy. Electrochemical storage systems showing even better dynamics than batteries are so-called supercapacitors. These devices can be quickly charged and discharged, but have lower energy density than batteries. Still, due to their favorable dynamics, they can be a valuable addition to batteries

in the framework of intermediate energy storage from renewables. The storage capacity of these devices depends largely on employed materials. In this special issue, development of novel electrode material for supercapacitor application based on pseudocapacitance is discussed. Guragain et al. [6] developed a large-surface-area MCO2O4 material in which a tubular microstructure leads to a noticeable pseudocapacitive property with the excellent specific capacitance value exceeding 407.2 F/g at 2 mV/s scan rate. In addition, a Coulombic efficiency ~100% and excellent cycling stability with 100% capacitance retention was noted even after 5000 cycles. These tubular MCO2O4 microstructures display peak power density exceeding 7000 W/kg. Based on these authors, the superior performance of the tubular MCO2O4 microstructure electrode is attributed to their high surface area, adequate pore volume distribution, and active carbon matrix, which allows effective redox reaction and diffusion of hydrated ions.

### **Electrolysis Processes Aiming to Solve Environmental Issues or for Construction Purposes**

In addition to clean energy, the requirement of clean water is of the upmost importance for prosperity and healthiness of the world population. With respect to water, one issue is water scarcity, causing 1.2 billion people to suffer from a lack of water [7]. Another issue is industry or domestic-based water pollution [7–9]. Due to the widespread use of chlorinated compounds such as HCl, NaCl, and MgCl2 in the industry, the content of chloride ion in wastewater is increasing [8,9]. If it is discharged directly to the environment, diverse environmental issues will be created with consequences for water resources, soil and human health. At present the most widely used method for Cl− waste water treatment is chemical precipitation. Electrochemical oxidation processes have a high potential for efficient removal of even trace amounts of Cl− due to high degradation efficiency, no secondary pollution, modularity, flexibility, and use of renewable electricity, which can contribute to stabilization of the energy grid. The degradation efficiency and degradation products of the electrochemical oxidation process change with the anode material, therefore, two articles by Xu et al. [8] and Liu et al. [9] discuss the development of new anode materials for electrochemical Cl− removal. Xu et al. [8] studied porous Ti/SnO2-Sb2O3/PbO2 electrodes for electrocatalytic oxidation of chloride ions by varying different parameters. The results have shown that Ti/SnO2-Sb2O3/PbO2 electrodes with 150 μm pore size had the best removal effect on chloride ions with removal ratios amounting up to 98.5% when the initial concentration of chloride ion was 10 g L−1. The advanced electrode structure minimized oxygen evolution as a side reaction, which further increased the removal effect of chloride ions. In the further publication, Liu et al. [9] studied different material combinations. The porous Ti/Sb–SnO2/Ni–Ce–PbO2 electrode was prepared by using a porous Ti plate as a substrate, an Sb-doped SnO2 as an intermediate, and a PbO2 doped with Ni and Ce as an active layer. The authors studied also the mechanism of electrochemical dechlorinating. They found out that the increase of OH− inhibits the degradation of Cl−.

In addition to inorganic catalysts, oxidation of mainly organic pollutants in the waste water treatment can be performed with the help of biological catalysts, for example microorganisms. These systems can operate either as microbial fuel cells (MFC) or in electrolysis mode. A contribution by Muddemann et al. [7] discuss current challenges in the scale up of these systems, with special emphasis on oxygen reducing cathode. The authors demonstrated a strong increase in the MFC performance and long term stability upon improving catalyst coating quality.

Finally Liu et al. [10] describe electrochemical discharge machining (ECDM) as an effective way to fabricate micro structures in non-conductive materials, such as quartz glass and ceramics. This has significant importance for the development of micro electromechanical systems (MEMS), such as micro reactors and micro medical devices, which often consist of the micro structures of nonconductive materials, such as glass, ceramics and silicon nitride and which are difficult to fabricate using the traditional machining methods.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The author declares no conflict of interest.

### **References**


© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Review* **Alkaline Water Electrolysis Powered by Renewable Energy: A Review**

### **Jörn Brauns \* and Thomas Turek**

Institute of Chemical and Electrochemical Process Engineering, Clausthal University of Technology, Leibnizstr. 17, 38678 Clausthal-Zellerfeld, Germany; turek@icvt.tu-clausthal.de

**\*** Correspondence: brauns@icvt.tu-clausthal.de; Tel.: +49-5323-72-2473

Received: 23 January 2020; Accepted: 18 February 2020; Published: 21 February 2020

**Abstract:** Alkaline water electrolysis is a key technology for large-scale hydrogen production powered by renewable energy. As conventional electrolyzers are designed for operation at fixed process conditions, the implementation of fluctuating and highly intermittent renewable energy is challenging. This contribution shows the recent state of system descriptions for alkaline water electrolysis and renewable energies, such as solar and wind power. Each component of a hydrogen energy system needs to be optimized to increase the operation time and system efficiency. Only in this way can hydrogen produced by electrolysis processes be competitive with the conventional path based on fossil energy sources. Conventional alkaline water electrolyzers show a limited part-load range due to an increased gas impurity at low power availability. As explosive mixtures of hydrogen and oxygen must be prevented, a safety shutdown is performed when reaching specific gas contamination. Furthermore, the cell voltage should be optimized to maintain a high efficiency. While photovoltaic panels can be directly coupled to alkaline water electrolyzers, wind turbines require suitable converters with additional losses. By combining alkaline water electrolysis with hydrogen storage tanks and fuel cells, power grid stabilization can be performed. As a consequence, the conventional spinning reserve can be reduced, which additionally lowers the carbon dioxide emissions.

**Keywords:** alkaline water electrolysis; hydrogen; renewable energy; sustainable; dynamic; fluctuations; wind; solar; photovoltaic; limitations

### **1. Introduction**

Hydrogen is considered a promising energy carrier for a sustainable future when it is produced by utilizing renewable energy [1]. Today, less than 4% of hydrogen production is based on electrolysis processes, of which the main part is hydrogen as a by-product of chlorine production. Hence, the major share of the needed hydrogen depends on the fossil path through the steam reforming of natural gas [2]. This situation is caused by the higher production costs of electrolysis processes compared to the conventional fossil sources, due to high electricity costs and interfering laws [3]. To reduce CO2 emissions and to become independent of fossil energy carriers, the share of hydrogen produced using renewable power sources needs to be increased significantly in the next few decades. Therefore, water electrolysis is a key technology for splitting water into hydrogen and oxygen by using renewable energy. After drying and removing oxygen impurities, the hydrogen purity is higher than 99.9%, and the hydrogen can be directly used in the following processes or in the transport sector [4]. Solar and wind energy are the preferred renewable power sources for hydrogen production, as their distribution is the most widespread [5,6]. Hydropower, biomass, and geothermal energy are alternatives, and are often utilized for the base load [7]. The main problem with using renewable energy is the unevenly distributed and intermittent local availability [6]. With a higher share of renewable energy from wind turbines or solar photovoltaic panels and fair CO2 emission costs, hydrogen production by water

electrolysis will become more attractive. The combination of water electrolysis with renewable energy is particularly advantageous, as excess electrical energy can be chemically stored in hydrogen to balance the discrepancy between energy demand and production [6]. For large-scale applications, the hydrogen can be stored in salt caverns, storage tanks, or the gas grid [8–12]. Smaller hydrogen quantities can also be stored in metal hydrides [13,14].

For water electrolysis, there are three technologies available: Alkaline water electrolysis (AEL), proton exchange membrane (or polymer electrolyte membrane) electrolysis (PEMEL), and solid oxide electrolysis (SOEL) [15–18]. While the low-temperature technologies, AEL and PEMEL, both provide high technology readiness levels, the high-temperature SOEL technology is still in the development stage [19]. Alkaline water electrolysis uses concentrated lye as an electrolyte and requires a gas-impermeable separator to prevent the product gases from mixing. The electrodes consist of non-noble metals like nickel with an electrocatalytic coating. PEMEL uses a humidified polymer membrane as the electrolyte and noble metals like platinum and iridium oxide as the electrocatalysts. Both technologies are operated at temperatures from 50 to 80 ◦C and allow operation pressures of up to 30 bar. The nominal stack efficiency of both technologies is around 70% [18,20]. SOEL is also known as high-temperature (HTEL) or steam electrolysis, as gaseous water is converted into hydrogen and oxygen at temperatures between 700 and 900 ◦C. Theoretically, stack efficiencies near 100% are possible due to positive thermodynamic effects on power consumption at higher temperatures. However, the increased thermal demand requires a suitable waste heat source from the chemical, metallurgical, or thermal power generation industry for economical operation. Moreover, the corrosive environment demands further material development [6,20,21]. As a consequence, SOEL provides only small stack capacities below 10 kW, compared to 6 MW for AEL and 2 MW for PEMEL [20]. Hence, the investment costs and the lifetime determine whether AEL or PEMEL is the most favorable system design for a large-scale application. Today, the investment costs for AEL are from 800 to 1500 € kW−<sup>1</sup> and for PEMEL from 1400 to 2100 € kW<sup>−</sup>1. Furthermore, the lifetime of alkaline water electrolyzers is higher and the annual maintenance costs are lower compared to a PEMEL system [15,20,22,23]. Often, PEMEL systems are preferred for dynamic operation due to the short start-up time and a broad load flexibility range. The shortcomings of AEL are gradually being overcome by further development [24]. Therefore, this review focuses on alkaline water electrolysis powered by renewable energy. To ensure safety and high efficiency, alkaline water electrolyzers must be optimized for dynamic operation. Hence, the process needs to be analyzed for how the dynamics will affect the system performance and what aspects should be considered when fluctuating renewable energy is used instead of a constant load [25]. Thus, this contribution shows model descriptions for alkaline water electrolysis, photovoltaic panels, and wind turbines to identify the limitations when combining all components into a hydrogen energy system. Furthermore, theoretical models can help to solve the existing problems using intelligent system design and suitable operation strategies.

This study mainly contains literature that was obtained with the keywords alkaline electrolyzer (or electrolyser or electrolysis) in combination with one of the following words: Renewable, sustainable, green, dynamic, fluctuation, intermittent, solar, photovoltaic, wind, and power to gas. For a broader overview, additional literature is also included. Figure 1 shows the number of annual publications that are listed in the Web of Science Database for the given keywords from 1990 to 2019. Additionally, the keyword alkaline is replaced by other water electrolysis technologies to show the share of technology-specific publications [26]. Around 2010, the number of annual publications started to increase persistently due to the discussion about the energy turnaround, especially in Germany and other European countries [9,27]. Furthermore, the topic is often discussed technology-independently, as the number of technology-specific publications is small compared to publications with unspecified water electrolysis technologies. While the low-temperature technologies, AEL and PEMEL, show an equal share of technology-specific publications, the high-temperature technology SOEL is mentioned less. This distribution reflects the recent considerations of which

technology may be favored for sustainable hydrogen production. Particularly, alkaline water electrolysis is considered as the most reliable method for large-scale hydrogen production [5,21].

**Figure 1.** The number of publications per year from 1990 to 2019 containing the specified keywords. Around 2010, the publication rate increases due to greater interest in the energy turnaround. While the topic is often discussed technology-independently (unspecified), more publications for low-temperature technologies, like alkaline water electrolysis (AEL) and proton exchange membrane electrolysis (PEMEL), are available than for the high-temperature technology solid oxide electrolysis (SOEL) [26].

### **2. Alkaline Water Electrolysis**

Alkaline water electrolysis is used to split water into the gases hydrogen and oxygen using electric energy. The chemical reactions are given in the Equations (1)–(3). At the cathode, water molecules are reduced by electrons to hydrogen and negatively charged hydroxide ions. At the anode, hydroxide ions are oxidized to oxygen and water while releasing electrons. Overall, a water molecule reacts to hydrogen and oxygen in the ratio of 2:1.

$$\text{Cathode:}\qquad 2\,\text{H}\_2\text{O}\_{\text{(l)}} + 2\,\text{e}^- \longrightarrow \text{H}\_{2\text{(g)}} + 2\,\text{OH}^-\_{\text{(aq)}}\tag{1}$$

$$\text{Anode:}\qquad 2\,\text{OH}^-\_{\text{(aq)}} \longrightarrow 0.5\,\text{O}\_{2(g)} + \text{H}\_2\text{O}\_{(l)} + 2\,\text{e}^-\tag{2}$$

$$\text{Overall reaction:}\qquad \mathrm{H}\_{2}\mathrm{O}\_{\text{(l)}} \longrightarrow \mathrm{H}\_{2\text{(g)}} + 0.5\,\mathrm{O}\_{2\text{(g)}}\tag{3}$$

The required cell voltage for this electrochemical reaction can be determined by thermodynamics. The free reaction enthalpy ΔR*G* in (4) can be calculated with the reaction enthalpy ΔR*H*, the temperature *T*, and the reaction entropy ΔR*S*.

$$
\Delta\_{\mathbb{R}}G = \Delta\_{\mathbb{R}}H - T \cdot \Delta\_{\mathbb{R}}S \tag{4}
$$

The reversible cell voltage *U*rev in (5) is determined by the ratio of the free reaction enthalpy ΔR*G* to the product of the number of exchanged electrons *z* = 2 and the Faraday constant *F* (96,485 C mol<sup>−</sup>1) [28].

$$\mathcal{U}\_{\text{rev}} = -\frac{\Delta\_{\text{R}}G}{z \cdot F} \tag{5}$$

At a temperature of 25 ◦C and an ambient pressure of 1 bar (standard conditions), the free reaction enthalpy for the water splitting reaction is ΔR*G* = 237 kJ mol<sup>−</sup>1, which leads to a reversible cell voltage of *U*rev = −1.23 V. As the free reaction enthalpy is positive at standard conditions, the water splitting is a non-spontaneous reaction [28]. Due to irreversibilities, the actual cell voltage needs to be higher than the reversible cell voltage for the water splitting reaction. The thermoneutral voltage *U*th in Equation (6) depends on the reaction enthalpy ΔR*H*, which is composed of the free reaction enthalpy ΔR*G* and irreversible thermal losses *T* · ΔR*S*.

$$\mathcal{U}\_{\text{th}} = -\frac{\Delta\_{\text{R}}H}{z \cdot F} \tag{6}$$

At standard conditions, the reaction enthalpy for water electrolysis is ΔR*H* = 286 kJ mol<sup>−</sup>1. Hence, the thermoneutral voltage is *U*th = −1.48 V [28].

### **3. System**

A schematic flow diagram of alkaline water electrolysis is shown in Figure 2. The electrolyte is pumped through the electrolysis stack, where the product gases are formed. While natural convection can be a cost-efficient alternative, gas coverage of the electrode surface can raise the required cell voltage and therefore increase the operational costs [29]. Additionally, most alkaline water electrolyzer systems provide a temperature control for the electrolyte to maintain an optimal temperature range.

**Figure 2.** A schematic flow diagram of an alkaline water electrolyzer. The electrolyte is pumped through the electrolysis cell where the gas evolution takes place. Adjacent gas separators split both phases, and the liquid phase flows back to the electrolysis stack. Heat exchangers ensure that the optimal temperature is maintained, and the product gases can be purified afterward.

The two-phase mixtures of liquid electrolyte and product gas leave the electrolysis cell and enter subsequent gas separators. Mostly, the phase separation is realized with a high residence time in large tanks. The product gas is demisted and dried before it is purified to the desired level [30]. The liquid electrolyte leaves the gas separator and is pumped back to the electrolysis stack. As the product gases are soluble in the electrolyte solution, the mixing of both electrolyte cycles causes losses and higher gas impurities. An alternative can be to use partly separated electrolyte cycles with an equalization line for liquid level balancing of both vessels [31,32]. With separated electrolyte cycles, the electrolyte concentration will increase on the cathodic side due to water consumption and decrease on the anodic side due to water production. Therefore, the electrolyte requires mixing, on occasion, to maintain an optimal electrolyte conductivity.

### **4. Cell Design and Cell Voltage**

The design of the electrolysis stack depends on the manufacturer; however, some general similarities can be observed. Two variants of cell designs are shown in Figure 3. Earlier alkaline water electrolyzers used a conventional assembly with a defined distance between both electrodes. Later, this concept was replaced by the zero-gap assembly, where the electrodes are directly pressed onto the separator to minimize ohmic losses due to the electrolyte. Porous materials like Zirfon™ Perl UTP 500 (AGFA) or dense anion exchange membranes can be used as the separator [33–37].

**Figure 3.** Different cell designs for alkaline water electrolysis. Whereas (**a**) shows a conventional assembly with a defined distance between both electrodes, (**b**) depicts a zero-gap assembly where the electrodes are directly pressed onto the separator [38].

During operation, the required cell voltage is always higher than the reversible cell voltage due to different effects. A calculated cell voltage profile is displayed in Figure 4. In addition to the ohmic losses, *I* · *R*ohm, there are activation overvoltages of the electrodes, *η*act. The ohmic resistance of the cell design is affected by the electronic conductivity of the electrode material, the specific conductivity of the electrolyte, the ionic conductivity of the separator material, and gas bubble effects.

**Figure 4.** The calculated cell voltage of an atmospheric alkaline water electrolyzer at a temperature of 60 ◦C according to Equation (8). The overall cell voltage consists of the reversible cell voltage *U*rev, ohmic losses *I* · *R*ohm, and activation overvoltages *η*act [39,40].

*Processes* **2020**, *8*, 248

The zero-gap design tries to eliminate the electrolyte losses by minimizing the electrode distance. There is still a minimal gap between both electrodes, which can increase the cell voltage. The activation overvoltages are defined by the electrode materials. Whereas nickel is the most-used electrode material, it provides very high overvoltages for the oxygen and hydrogen evolution reactions [41–44]. Hence, electrocatalytic materials are added to the electrodes. Iron is a cost-efficient catalyst for the oxygen evolution reaction [41,42,45]. Molybdenum decreases the overvoltage for the evolution of hydrogen at the cathode [44,46,47].

Several authors have proposed correlations for the modeling of cell voltage. Equation (7) considers the operation temperature *ϑ* and the current density *j* by describing the dependencies with empirical parameters. While the parameters *ri* reflect ohmic losses, *s* and *ti* stand for the activation overvoltages of the hydrogen and oxygen evolution reactions [28].

$$\mathcal{U}\_{\rm cell} = \mathcal{U}\_{\rm rev} + \left(r\_1 + r\_2 \cdot \vartheta\right) \cdot j + s \cdot \log\left[\left(t\_1 + \frac{t\_2}{\vartheta} + \frac{t\_3}{\vartheta^2}\right) \cdot j + 1\right] \tag{7}$$

This correlation can be extended with the effects of the operation pressure *p* in (8) by adding the empirical parameters *di*, which specify the additional losses owing to pressurized operation [39]. In general, the reversible cell voltage increases with the pressure; however, the ohmic resistance caused by the gas bubbles decreases as the bubble diameter becomes smaller. Hence, both effects equalize each other and only small differences can be observed [48].

$$\mathcal{U}\_{\rm cell} = \mathcal{U}\_{\rm rev} + \left[ (r\_1 + d\_1) + r\_2 \cdot \theta + d\_2 \cdot p \right] \cdot j + s \cdot \log \left[ \left( t\_1 + \frac{t\_2}{\theta} + \frac{t\_3}{\theta^2} \right) \cdot j + 1 \right] \tag{8}$$

The correlations (7) and (8) are empirical and therefore only valid for the actual system to which they are adjusted. The correlation parameters and a suitable equation for the reversible cell voltage under atmospheric conditions can be found in the Appendix A in Table A1 and Equation (A1). Other authors have proposed physically reasonable models based on actual dimensions and properties of the system rather than on empirical correlations.

An example for such an approach is Equation (9), in which the terms are split into experimentally determinable parts [49].

$$\mathcal{U}\_{\rm cell} = \mathcal{U}\_{\rm rev} + \eta\_{\rm act}^{\rm c} + \eta\_{\rm act}^{\rm a} + I \cdot (\mathcal{R}\_{\rm c} + \mathcal{R}\_{\rm a} + \mathcal{R}\_{\rm elec} + \mathcal{R}\_{\rm mmem}) \tag{9}$$

The cell voltage *U*cell is calculated with the reversible voltage *U*rev, the activation overvoltages *η*act, and the ohmic resistances. Whereas *R*<sup>c</sup> and *R*<sup>a</sup> represent the reciprocal electronic conductivity of the electrode materials, *R*ele stands for the ohmic loss caused by the electrolyte conductivity. Additionally, the ohmic resistance *R*mem of the separator material is taken into account. The activation overvoltages *η*act can be calculated with the Butler–Volmer equation. In most cases, the simplified Tafel equation is sufficient to describe the resulting overpotentials [40]. The required Tafel slope and exchange current density can be extracted from experimental data. Hence, those parameters are only valid for the actual system design; however, they can be easily replaced by other data when needed. As the ohmic resistances of the electrodes (*R*<sup>c</sup> and *R*a) only depend on the electronic conductivity and electrode dimensions, both values are known. In most cases, the ohmic resistance of the electrode is comparably small and can be neglected. The electrolyte resistance *R*ele is determined by the specific electrolyte conductivity and the cell design. Whereas the electrolyte gap is minimal in zero-gap designs, conventional setups maintain a defined distance between both electrodes. As the specific conductivity of the electrolyte gap is affected by gas bubbles, there is an optimal electrode distance for conventional designs [50]. If the electrode distance is too small, the gas bubbles accumulate in the gap and lower the conductivity. With increasing distance, the bubble detachment is enhanced and the specific conductivity increases. It is a trade-off between a small electrolyte gap—as the ohmic resistance increases linearly with this parameter—and a better conductivity of the space between both electrodes. In addition to the decreasing electrolyte conductivity with higher amounts of gas bubbles, the active electrode surface can be blocked by gaseous compounds, which leads to additional losses [49]. As this phenomenon depends on the cell design and operation concept, there are difficulties in describing it properly. Therefore, it is often neglected, or empirical correlations referring to the gas hold-up are utilized [49].

Furthermore, the installed separator material also has significant ohmic losses. While the porous separator Zirfon™ Perl UTP 500 is often used, anion exchange membranes are promising alternatives. For Zirfon™-based materials, experimental data of the resistance at a fixed electrolyte concentration for different temperatures are available [51].

The most-used electrolyte for alkaline water electrolysis is an aqueous solution of potassium hydroxide (KOH) with 20 to 30 wt.% KOH, as the specific conductivity is optimal at the typical temperature range from 50 to 80 ◦C [25]. A cheaper alternative would be a diluted sodium hydroxide solution (NaOH), which has a lower conductivity [52]. Calculated specific electrolyte conductivities for both electrolyte solutions at different temperatures are shown in Figure 5. While KOH provides a specific conductivity around 95 S m−<sup>1</sup> at 50 ◦C, NaOH reaches a value around 65 S m<sup>−</sup>1. At a temperature of 25 ◦C, a similar effect can be seen. The conductivity of KOH is around 40 to 50% higher than the conductivity of a NaOH solution at the optimal weight percentage. Another aspect is the solubility of the product gases inside the electrolyte, as this influences the resulting product gas purity. In general, the gas solubility decreases with an increasing electrolyte concentration due to the salting-out behavior [53]. NaOH also shows a slightly higher salting-out effect than that of KOH. Hence, the product gas solubility is higher in a KOH solution [54–56].

**Figure 5.** The calculated specific electrolyte conductivity as a function of the electrolyte concentrations of sodium hydroxide (NaOH) and potassium hydroxide (KOH) solutions at different temperatures obtained by Equations (A2) and (A3). The correlation parameters can be found in Table A2 [52,57].

Another approach is to use ionic liquids (ILs) as the electrolyte or as an additive, owing to their remarkable properties [5,6,21]. Ionic liquids are organic substances which are liquid at room temperature and are electrically conductive [58]. A negligible vapor pressure, non-inflammability, and thermal stability are promising arguments for their utilization in water electrolysis. Furthermore, ILs can be used over a wide electrochemical window [59]. The absorption and separation of gases is an additional area of application [60,61]. However, the toxicity of ILs is a current field of research, and the viscosity is comparably high, which should be taken into account before any large-scale implementation [6,58,59]. In addition to providing high-efficiency water electrolysis at low temperatures, ILs are chemically inert and therefore do not require expensive electrode materials [62].

### **5. Gas Purity**

Gas purity is an important criterion of alkaline water electrolysis. While the produced hydrogen typically has a purity higher than 99.9 vol.% (without additional purification), the gas purity of oxygen is in the range of 99.0 to 99.5 vol.% [48]. As both product gases can form explosive mixtures in the range of approximately 4 to 96 vol.% of foreign gas contamination, technical safety limits for an emergency shutdown of the whole electrolyzer system are at a level of 2 vol.% [31,63]. Therefore, the product gas impurity needs to be below this limit during operation to ensure continuous production. Experimentally determined anodic gas impurities for alkaline water electrolysis are presented in Figure 6 for different operation modes. The current densities are in the range from 0.05 to 0.7 A cm−<sup>2</sup> and the system pressures range from 1 to 20 bar [64].

**Figure 6.** Anodic gas impurity (H2 in O2) in relation to the current density at different pressure levels for (**a**) separated and (**b**) mixed electrolyte cycles, at a temperature of 60 ◦C, with an electrolyte concentration of approximately 32 wt.% and an electrolyte volume flow of 0.35 L min−<sup>1</sup> [64].

While the gas impurities with separated electrolyte cycles are below 0.7 vol.% for all tested current densities and pressure levels, mixing of the electrolyte cycles increases the gas impurity significantly. Furthermore, two similarities can be seen. The gas impurity lowers with increasing current density and increases at higher pressure levels. Both effects are physically explainable. While the contamination flux stays constant with varying current densities, the amount of produced gas becomes lower in a linear relationship. Hence, at a higher current density, the contamination is more diluted than at a lower current density [32,64]. As a consequence, the operation in the part-load range is more critical due to the higher gas impurity. The amount of dissolved product gas increases with pressure; thus, high concentration gradients for the diffusion through the separator material are available, and more dissolved foreign gas reaches the other half-cell when mixing [64]. However, operation at slightly elevated pressures is favorable, as the costly first mechanical compression level can be avoided by the direct compression inside the electrolyzer system [65]. With mixed electrolyte cycles, the gas impurity reaches critical values even at higher current densities during pressurized operation. While at atmospheric pressure, the gas impurity is only at a current density of 0.05 A cm<sup>−</sup>2, slightly above the safety limit of 2 vol.% H2 in O2, this limit is already reached at 0.5 A cm−<sup>2</sup> for a system pressure of 10 bar. At 20 bar, no sufficient gas purity could be measured, as even a current density of 0.7 A cm−<sup>2</sup> results in a gas impurity of 2.5 vol.%.

### **6. Periphery**

The operation of alkaline water electrolyzers is also affected by the installed periphery, including power supplies. The output signal may contain a specific number of ripples, which directly influences the process performance [66]. Power supplies with a high ripple propensity lower the overall efficiency and, therefore, signal smoothing is necessary. The ripple formation is avoidable, but the component costs will be higher [67]. In general, thyristor-based power supplies tend to deliver a higher degree of fluctuation, and transistor-based systems output a smoother signal. Additionally, a higher ripple frequency does not affect the system performance as much as the occurrence of low-frequency ripples [68]. Furthermore, a coherence between the ripple behavior of a power supply and the product gas quality of alkaline water electrolysis can be observed [69].

### **7. Renewable Energy**

The combination of alkaline water electrolysis with renewable energy is essential for sustainable hydrogen production without significant carbon dioxide emissions. While solar and wind energy are often favored due to their wide availability, other renewable energies, such as hydropower, biomass, and geothermal energy, are frequently utilized for the base load [7]. The direct usage of renewable energy in the power grid is difficult due to the mismatch between energy demand and production and the limited storage possibilities for electricity. Hence, excess electric energy should be chemically stored in hydrogen for later usage [6]. Due to the fluctuating and intermittent behavior of solar and wind power, alkaline water electrolyzers must be adapted to a dynamic operation. To evaluate the requirements, local weather data can be used to extract the amplitudes and frequencies of fluctuations.

Typical time-related profiles for solar radiation and wind velocity are shown in Figure 7. The data were measured by the weather station of the Clausthal University of Technology on the rooftop of a university building. Whereas the wind velocity shows a mean value of around 3.8 m s<sup>−</sup>1, the significant solar radiation is only available during the daytime. Hence, the averaged value over the whole day is 233 W m−<sup>2</sup> for a sunny day and only 29 W m−<sup>2</sup> for a cloudy day.

**Figure 7.** Typical time-related profiles for (**a**) solar radiation and (**b**) wind velocity, measured by the weather station of the Clausthal University of Technology. Though solar radiation peaks around noon, wind velocity shows sinusoidal oscillations.

The volume flow of the produced hydrogen directly follows the renewable energy profile used for operation [70]. Only a short delay is noticeable for the gas purity, which is defined by the system volume [71]. Due to the possibility of direct coupling of water electrolysis and photovoltaic panels, this technology is highly appropriate for renewable hydrogen production [29,72,73]. As photovoltaic

panels require high investment costs, wind power is often favored for large-scale hydrogen production. In comparison with photovoltaic power, wind power shows a higher degree of fluctuation and is very intermittent. Therefore, the dynamic operation of alkaline water electrolyzers is more challenging [4].

Hence, the dynamic behavior of an alkaline water electrolyzer can be used to develop suitable system designs and to operate existing systems safely and efficiently. As measurements of solar radiation and wind velocity are often available for a given location, the theoretically available renewable energy can be calculated and used as an input during the system design. Different approaches exist for the calculation of solar photovoltaic power and wind turbine power. While the current–voltage characteristics of photovoltaic panels can be expressed as a function of manufacturer data and solar radiation, the power of wind turbines is a fraction of the maximum available wind power, which is defined by the wind speed and a performance coefficient [72,74].

### *7.1. Solar Photovoltaic Power*

The behavior of photovoltaic panels can be described by single-diode and two-diode models with varying degrees of complexity. Often, the solution must be obtained iteratively or with numerical methods when very detailed models are utilized [75,76]. Simple models with analytical solutions are a recent research topic, as a short processing time can be needed for online characterization and the optimization of existing systems [75]. In Figure 8, the coupling possibilities of an alkaline water electrolyzer and solar photovoltaic panels are shown. Additional losses of a DC/DC transformer can be avoided when a direct coupling of the systems can be realized. Otherwise, the transformation ensures a fit of both systems by an indirect coupling [73,77,78].

**Figure 8.** Schematic of alkaline water electrolysis powered by solar energy. Photovoltaic panels convert the solar radiation into electricity, which can be used for the operation. The implementation of a DC/DC power converter is optional, as direct and indirect coupling is possible [70,78,79].

When a direct coupling of both systems is to be realized, the possible operation points can be determined by the intersection of the current–voltage curves. A typical current–voltage characteristic of an alkaline water electrolyzer is given by (8). The resulting current of a photovoltaic cell *I*PV at different solar radiation levels can be described by (10) with a suitable single-diode model as a function of the voltage *U*PV [29,72,73]. Therefore, specific data from the photovoltaic (PV) panel and the ambient conditions are required in order to calculate the photocurrent *I*ph, the reverse saturation current *I*s, and the thermal voltage *U*T. Furthermore, the serial *R*<sup>s</sup> and parallel *R*<sup>p</sup> resistance of the photovoltaic panel must be available.

$$I\text{pV} = I\_{\text{Ph}} - I\_{\text{s}} \cdot \left[ \exp\left(\frac{\text{LI}\_{\text{PV}} + I\_{\text{PV}} \cdot R\_{\text{s}}}{\text{LI}\_{\text{T}}}\right) - 1\right] - \frac{\text{LI}\_{\text{PV}} + I\_{\text{PV}} \cdot R\_{\text{s}}}{R\_{\text{p}}} \tag{10}$$

The photocurrent *I*ph is defined in (11), which shows a linear relationship with the solar radiation *E*sun absorbed by the photovoltaic cell. A higher cell temperature *T*c increases the photocurrent.

$$I\_{\rm ph} = \left(0.003 \,\mathrm{m}^2 \,\mathrm{V}^{-1} + 10^{-7} \mathrm{m}^2 \,\mathrm{V}^{-1} \,\mathrm{K}^{-1} \cdot T\_{\rm c}\right) \cdot E\_{\rm sun} \tag{11}$$

The reverse saturation current *I*<sup>s</sup> can be calculated by (12) with the short-circuit current *I*sc, the open-cell voltage *U*oc, and the thermal voltage *U*T. Whereas the short-circuit current and the

open-cell voltage are provided by the manufacturer, the thermal voltage depends on the physical properties.

$$I\_8 = \frac{I\_{\rm sc}}{\exp\left(\frac{UL\_{\rm sc}}{U\_{\rm T}}\right) - 1} \tag{12}$$

An equation for the thermal voltage is given in (13), which is based on the Boltzmann constant *<sup>k</sup>*<sup>B</sup> (1.3806·10−<sup>23</sup> J K<sup>−</sup>1) and the electron charge *<sup>e</sup>* (1.602 19·10−<sup>19</sup> C) [72]. Additionally, the number of serially connected cells, *n*s, and the cell temperature are required. Furthermore, the non-ideality factor *m* contains any deviations from the theoretical behavior.

$$
\Delta L\_{\rm T} = m \cdot \frac{n\_{\rm s} \cdot k\_{\rm B} \cdot T\_{\rm c}}{\varepsilon} \tag{13}
$$

In addition to these equations, the calculation of the resulting current of a photovoltaic cell requires knowledge of the serial (*R*s) and parallel (*R*p) resistance of the system. By adding parallel photovoltaic cells, the current multiplies by the amount of parallel paths *n*p. Suitable parameters of an existing photovoltaic cell setup are given in Table 1. For this exemplary calculation, a constant temperature of the photovoltaic cell is assumed. Otherwise, the cell temperature increases with the absorbed solar radiation. While simple linear approaches already result in a good agreement with experimental data, a complete energy balance is the best way to determine the temperature exactly [29,72].

**Table 1.** Parameters for the example calculation of the photovoltaic current using Equation (10). The number of serial *n*s and parallel *n*p connected photovoltaic cells, the short-circuit current *I*sc, the open-cell voltage *U*oc, the serial *R*s and parallel resistance *R*p, and the non-ideality factor *m* are setup-specific data. A constant cell temperature *T*c is assumed [29,72,73].


The results of the example calculation are shown in Figure 9. The current–voltage characteristics are given for different solar radiation levels from 200 to 1000 W m<sup>−</sup>2, in combination with a typical polarization curve of an alkaline water electrolyzer (10 cm2 electrode area) from (8) in Figure 9a. The power–voltage curves for the photovoltaic cell are shown in Figure 9b. The maximal power point (MPP) for each radiation level is marked with a dot in both diagrams.

In Figure 9a, the characteristics of the alkaline water electrolyzer deviate from the MPP curve. Therefore, the photovoltaic cell cannot deliver the maximal power, and the overall efficiency decreases. Hence, both systems should be optimized until the alkaline water electrolyzer performs close to the maximal power output [73,80]. The alternative would be an indirect coupling of both systems with the integration of a DC/DC converter, which also implies losses, with an efficiency of around 90% [81,82].

**Figure 9.** Example calculation results of the (**a**) current–voltage characteristics of a photovoltaic panel at different solar radiation levels and the corresponding (**b**) power–voltage curve. Additionally, a current–voltage characteristic of an alkaline water electrolyzer (AEL) is implemented. The intersections determine the possible operation points. For an efficient operation, the distance to the maximal power points (MPP) should be minimal [29,72,73].

### *7.2. Wind Power*

As the power from photovoltaic cells is only available during the daytime, wind power is another important energy source for the renewable production of hydrogen. The schematic concept is shown in Figure 10. For the implementation of conventional wind turbines, an AC/DC converter is essential. The efficiency of an AC/DC conversion is also approximately 90% [82,83].

**Figure 10.** Schematic of alkaline water electrolysis powered by wind energy. Wind turbines convert the available wind power into electricity, which can be used for the operation. The implementation of a suitable AC/DC converter is mandatory [74,79].

For the calculation of the wind turbine power, the exact wind velocity at the height of the turbine rotor should be known. Often, the wind velocity is measured at rooftops or special measurement facilities with a defined height of approximately 10 m, which is significantly lower than the height of a wind turbine, around 100 m [84]. Therefore, the measured data should be corrected to the desired height by (14).

$$\upsilon\_{\text{wind}} = \upsilon\_{\text{wind,ref}} \cdot \frac{\ln\left(\frac{z\_{\text{wind}}}{z\_0}\right)}{\ln\left(\frac{z\_{\text{wind,ref}}}{z\_0}\right)}\tag{14}$$

The wind velocity *v*wind at the height *z*wind can be determined from the measured wind velocity *v*wind,ref at the height *z*wind,ref in combination with the roughness of the terrain *z*<sup>0</sup> [48]. To obtain the output power of a wind turbine *P*turbine, first, the theoretical wind power *P*wind needs to be calculated using (15). Therefore, the air density *ρ* (from 1.22 to 1.3 kg m<sup>−</sup>3), the area spanned by the rotor blades *A*, and the wind velocity are needed [74,85].

$$P\_{\rm wind} = \frac{1}{2} \cdot \rho \cdot A \cdot v\_{\rm wind}^3 \tag{15}$$

The maximal wind power cannot be completely converted into wind turbine power. This circumstance is considered by the implementation of the performance coefficient *C*P, which lowers the maximal reachable power output. The actual wind turbine power results from the product of the wind power and the performance coefficient in (16).

$$P\_{\text{turbine}} = P\_{\text{wind}} \cdot \mathcal{C}\_{\text{P}} \tag{16}$$

The determination of the correct performance coefficient is a complete research topic in itself, which consists of empirical correlations and computational fluid dynamics (CFD) simulation studies. Often, experimental data are used to fit the correlations to the measurements [74]. An example equation for the performance coefficient is shown in (17) [74,79].

$$\mathcal{C}\_{\rm P} = 0.22 \cdot \left( \frac{116}{\lambda\_{\rm i}} - 0.4 \cdot \beta - 5 \right) \cdot \exp \left( -\frac{12.5}{\lambda\_{\rm i}} \right) \tag{17}$$

Therefore, the pitch angle of the turbine blades *β* has to be defined and the tip speed ratio *λ* needs to be calculated in (18) from the turbine blade radius *R*, the rotational speed *ω*, and the wind speed [74].

$$
\lambda = \frac{R \cdot \omega}{v\_{\text{wind}}} \tag{18}
$$

The calculation of the performance coefficient also requires the parameter *λ*i, which is described by (19) based on the tip speed ratio and the blade pitch angle [74].

$$\frac{1}{\lambda\_{\text{i}}} = \frac{1}{\lambda + 0.08 \cdot \beta} - \frac{0.035}{\beta^3 + 1} \tag{19}$$

For the blade radius, a value of 46.5 m is assumed, which is a typical blade length for a wind turbine with a rated power of 2 MW [74]. In Figure 11, the calculation results for the performance coefficient, depending on the tip speed ratio and the turbine power at different wind velocities, are shown. The performance coefficient of conventional wind turbines is limited at *C*<sup>p</sup> = 0.593 [74]. In this example, a maximal performance coefficient of approximately *C*<sup>p</sup> = 0.450 is reached for a blade pitch angle of *β* = 0°. With an increasing pitch angle, the maximum of the performance coefficient decreases and shifts towards smaller tip speed ratios.

**Figure 11.** Example calculation results of (**a**) the performance coefficient for various rotor blade pitch angles using Equation (17) and (**b**) the wind turbine power for different wind velocities using Equation (16). The maximum power point (MPP) trajectory is marked [74,79].

For the calculation of the turbine power in Figure 11b, a pitch angle of *β* = 6° is assumed. With increasing wind velocity, the value of the maximal power point (MPP) becomes higher and shifts towards faster rotational speeds. The rated wind speed of this exemplary wind turbine is at 11 m s−<sup>1</sup> with rotational speeds from 6 to 17 min<sup>−</sup>1. The cut-in wind speed is 3ms−<sup>1</sup> and the cut-out wind speed is 22 m s−<sup>1</sup> [74]. In comparison with the power characteristics of photovoltaic panels, the polarization curve of alkaline water electrolyzers can not be directly optimized towards the MPP trajectory, as the optimal operation point highly depends on the wind turbine design and weather conditions. Therefore, an efficient AC/DC converter is the best option for maintaining an efficient operation of an alkaline water electrolyzer [82].

### **8. Hydrogen Energy System and Power Grid Stabilization**

An exemplary process scheme for a hydrogen energy system is provided in Figure 12. Photovoltaic panels and wind turbines are connected with suitable converters to a DC bus, from which alkaline water electrolyzers are powered. The produced hydrogen can be stored for later application in fuel cells. To raise the fuel cell efficiency, the produced oxygen can be used instead of air. Therefore, an additional storage tank must be available, which incurs further costs [86].

The fuel cells are also connected to the DC bus, and the power can be used by the electricity grid with DC/AC converters. At lower energy demands, hydrogen can be produced and converted back into energy when it is needed. As conventional alkaline water electrolyzers are designed for operation at constant conditions, occurring fluctuations may be damped by additional energy storage devices like batteries, supercapacitors, or flywheels [25,28,82]. When excess energy is available, this energy storage can be charged to be fully available when needed. The damping quantity is limited to a certain degree of fluctuation, as the energy storage amount is also restricted to the capacity of all installed devices. Additionally, the produced hydrogen can also be used for the decarbonization of industrial processes or as a fuel in the transport sector [87–89]. To raise the overall efficiency, some DC/DC converters could be neglected by optimized system designs by lowering the system flexibility. Furthermore, when the alkaline water electrolyzers are able to operate under dynamic conditions, additional energy storage devices are not required or, at least, the number of such devices could be lowered. There are still some challenges for electrolyzer manufacturers to overcome before this possibility becomes available.

**Figure 12.** The schematic process scheme of a hydrogen energy system. Photovoltaic panels and wind turbines generate renewable energy to power alkaline water electrolyzers, and stored hydrogen can be converted back into electricity by fuel cells. Therefore, either oxygen or air can be utilized. Additional energy storage devices can damp fluctuations, and the complete hydrogen energy system can be used for power grid stabilization [25,28,82,87].

With an increasing share of renewable energies in the power grid, it is difficult to maintain a constant power frequency. Such hydrogen energy systems or alkaline water electrolyzers can be used to stabilize the power frequency by damping the fluctuations. An additional benefit would be the reduction of the conventional spinning reserve, which reduces costs and CO2 emissions [87,90]. A predictive control can be used for stable and efficient operation. Pressurized alkaline electrolyzers are more suitable for damping fast fluctuations, whereas atmospheric units can handle the slow fluctuations [87].

### **9. Limitations and Solution Approaches**

The implementation of a hydrogen energy system into the existing power grid is a challenging task with some limitations which must be overcome in order to guarantee high system availability. The main problem of an alkaline water electrolyzer powered by renewable energy is the high gas impurity in the part-load range, which can cause a safety shutdown when reaching a foreign gas contamination of 2 vol.% [31,91]. Hence, the annual operation time is limited to the time spans with sufficient renewable energy [91].

### *9.1. Limited Operation Time*

The limited operation time leads to a high number of startup and shutdown cycles, which can exceed the maximal start/stop count defined by the manufacturer and, therefore, can lower the expected system lifetime or warranty agreements. Mainly, the electrodes are affected by the repetitive start/stop behavior and the electrode degradation is accelerated [48,82]. Nickel electrodes are known to degrade significantly after 5000 to 10,000 start/stop cycles. When operating with photovoltaic power only, 7000 to 11,000 cycles are already reached in the period of 20 to 30 years. The fluctuating nature of renewable energy amplifies the electrode degradation, as this phenomenon acts partly as a start/stop process [92]. This issue can be solved by the development of stable electrode compositions or self-repairing electrode surfaces [92].

To circumvent the drawbacks of having only one renewable power source, such as in the daytime-limited operation with solar power, the combination of several energy sources enhances the overall efficiency. While the operation with only PV shows a faradaic efficiency of approximately 40%, wind power leads to a faradaic efficiency of around 80%. The combination of both technologies enhances the faradaic efficiency above 85% [79].

To hinder the gas impurity from reaching the lower explosion limit, the part-load range of most alkaline electrolyzers is limited to 10 to 25% of their nominal load [82,91]. Fluctuations below the minimal load can be balanced out with the implementation of energy storage devices, as shown in Figure 12; however, in some scenarios, the available energy storage will not be sufficient. When the gas impurity is still in a tolerable region, short periods without an electrode polarization can be allowed. The cathode starts to degrade noticeably below a voltage of around 0.25 V [82]. Thus, the complete shutdown can be held until reaching this voltage limit. The available time depends on the electrode composition, as the electrochemical double layer acts as a capacitor and delays the voltage breakdown after a power loss. Experimentally, a time span of around 10 min has been reported [82].

#### *9.2. Optimal System Design and Operation Strategies*

To mitigate the rise of gas impurities during low power availability, an optimal system design can allow enough time until sufficient energy is available again. While the gas volume inside the system acts as a buffer tank and dilutes the gas contamination, the liquid and the solid volume of an electrolyzer buffers the system temperature during part-load operation [25,71].

Furthermore, to maintain an efficient operation, the system temperature has to be in an optimal range of 50 to 80 ◦C for an electrolyte solution with 20 to 30 wt.% KOH [25]. As most renewable-energy-powered alkaline water electrolyzers will not provide a separate heating unit, the temperature needs to be reached and maintained only by the heat of the reaction [4]. Temperatures above 80 ◦C should be avoided with a suitable cooling system to prevent high degradation rates. An alternative would be the operation at low temperatures to damp electrode degradation, but then, very active electrocatalysts are needed to reach a sufficient efficiency [86].

More experimental and theoretical work is needed to fully understand the dynamic behavior of alkaline water electrolyzers powered by renewable energy [25]. In addition to an optimal system design, suitable dynamic operation strategies can be beneficial for lowering the gas impurity. While low gas impurities occur with separated electrolyte cycles, high gas impurities result in combined mode. The measured stationary gas impurities in Figure 6 are reached after a specific duration. When the electrolyzer is able to switch between both operation modes automatically, this can be used to switch to the separated mode when the gas impurity is too high, and then combine again when a sufficient gas production rate is available. Experimental work shows that the gas impurity can be almost halved by this approach [31].

The primary reason for high gas contamination is the continuous operation at low current densities. This circumstance can be prevented by reducing the overall cell area (overloading) or by subdividing the system into several smaller blocks [91]. While the implementation of electrolyzers with smaller electrode areas also limits the maximal load compared to larger systems, partial system operation is a more elegant method. During low power availability, single stacks or compartments of a system with multiple stacks can be powered off, which lowers the available electrode area and therefore results in higher current densities [93]. Obviously, this strategy causes problems in maintaining the optimal system temperature due to the disabled components. An alternative method to prevent adverse process states is the use of predictive control systems. For example, when low renewable power availability is forecasted, the system can change the temperature, pressure, or operation mode to a more suitable state before negative effects occur [87].

### **10. Conclusions**

The combination of alkaline water electrolysis and renewable energy for sustainable hydrogen production is an essential step towards the decarbonization of industrial processes and the transport sector [87–89]. To determine the most relevant limitations and to propose suitable solution approaches, the technologies have to be fully understood [25]. Whereas the process of alkaline water electrolysis can be defined by current–voltage characteristics and the resulting gas impurity, photovoltaic panels and wind turbines should be operated at the maximal power point [73,74,79]. Therefore, the influencing parameters must be known. Different model approaches exist, out of which the most suitable one should be chosen. While empirical correlations are often only valid for the specific experimental setup, physically reasonable models can be used in a more general way to develop new solutions. For alkaline water electrolysis, many experimental and theoretical data are available to calculate and analyze the cell voltage under operation conditions. As the actual system design and cell arrangement differ for every electrolyzer, certain parameters have to be determined experimentally to use the proposed models for another system. Mainly, this issue exists for electrode compositions and separator materials. To describe the gas purity of hydrogen and oxygen mathematically, only models and correlations on an empirical basis are currently available due to the high number of influencing variables [31,32]. As the gas impurity mainly determines the system availability of an alkaline water electrolyzer, more research for the development of physically-based models is needed. The dynamic system behavior should be analyzed, as optimized dynamic operation strategies can be beneficial for the overall system efficiency. Many models with different complexity levels are available for the description of the current–voltage characteristics of photovoltaic panels. Most models rely on physical principles and manufacturer data [75]. Thus, proper modeling for different systems is possible. The power conversion by wind turbines can be described by system properties and suitable correlations for the performance coefficient [74]. As this variable is influenced by many parameters, including the design of the turbine blades, the correlation should only be used for very similar wind turbines, or the parameters must be determined experimentally or by simulation.

To conclude, there are appropriate models available for all components of a hydrogen energy system. However, some descriptions need further improvement to be applicable to a variety of different system designs. With this knowledge and with experimental studies, many researchers have already examined the limitations of renewable-powered alkaline water electrolyzers [48,79,82]. The central prospect is to increase the operation time through intelligent system designs and advantageous operational concepts. While the implementation of conventional energy storage devices to damp the dynamics is a first logical step, alkaline water electrolyzers should be enabled to handle all dynamics directly to reduce costs and to enhance the efficiency [25]. As the hydrogen production from fossil energy carriers is less expensive than hydrogen from electrolysis processes, only optimized systems with the use of excess renewable energy can be competitive.

**Author Contributions:** Conceptualization, methodology, software, validation, formal analysis, investigation, data curation, writing—original draft preparation, visualization, J.B.; resources, writing—review and editing, J.B. and T.T.; supervision, project administration, funding acquisition, T.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Project numbers: 290019031; 391348959.

**Acknowledgments:** The authors thank the Institute of Electrical Information Technology (IEI) of the Clausthal University of Technology for providing the weather data.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

### **Abbreviations**


#### **Appendix A. Correlations and Parameters**

A correlation for the reversible cell voltage *U*rev of alkaline water electrolysis is given in (A1). The obtained value can be used for the calculation of the cell voltage in (7) or (8) at atmospheric conditions. For a pressurized system, extended correlations are required, as the reversible cell voltage increases at higher pressures [40]. The empirical correlation parameters for the calculation of cell voltage by (7) and (8) are given in Table A1.

$$\mathrm{U}\_{\mathrm{rev}} = 1.503 \,\mathrm{42\,V} - 9.956 \cdot 10^{-4} \,\mathrm{V} \cdot \left(\frac{T}{K}\right) + 2.5 \cdot 10^{-7} \,\mathrm{V} \cdot \left(\frac{T}{K}\right)^2\tag{A1}$$



The correlations for the calculation of specific electrolyte conductivity for KOH and NaOH can be found in (A2) and (A3). The required correlation parameters are listed in Table A2. The validity range for (A2) is a temperature *T* from 258.15 to 373.15 K and KOH mass fractions *w*KOH between 0.15 and 0.45. Equation (A3) is valid for temperatures *ϑ* between 25 and 50 ◦C and NaOH mass fractions *w*NaOH from 0.08 to 0.25 [52,57].

$$\begin{aligned} \boldsymbol{\sigma\_{\rm KOH}} &= \mathbf{K\_{1}} \cdot (100 \cdot \boldsymbol{w\_{\rm KOH}}) + \mathbf{K\_{2}} \cdot \boldsymbol{T} + \mathbf{K\_{3}} \cdot \boldsymbol{T^{2}} + \mathbf{K\_{4}} \cdot \boldsymbol{T} \cdot (100 \cdot \boldsymbol{w\_{\rm KOH}}) \\ &+ \mathbf{K\_{5}} \cdot \boldsymbol{T^{2}} \cdot (100 \cdot \boldsymbol{w\_{\rm KOH}})^{\mathbf{K\_{6}}} + \mathbf{K\_{7}} \cdot \frac{\boldsymbol{T}}{(100 \cdot \boldsymbol{w\_{\rm KOH}})} + \mathbf{K\_{8}} \cdot \frac{(100 \cdot \boldsymbol{w\_{\rm KOH}})}{\boldsymbol{T}} \end{aligned} \tag{A2}$$

$$\sigma\_{\text{NaOH}} = K\_1 + K\_2 \cdot \theta + K\_3 \cdot w\_{\text{NaOH}}^3 + K\_4 \cdot w\_{\text{NaOH}}^2 + K\_5 \cdot w\_{\text{NaOH}} \tag{A3}$$


**Table A2.** Parameters for the calculation of the specific electrolyte conductivities of KOH and NaOH solutions by Equations (A2) and (A3) [52,57].

### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Pore Network Simulation of Gas-Liquid Distribution in Porous Transport Layers**

### **Nicole Vorhauer 1,\*, Haashir Altaf 1, Evangelos Tsotsas <sup>1</sup> and Tanja Vidakovic-Koch <sup>2</sup>**


Received: 13 July 2019; Accepted: 20 August 2019; Published: 23 August 2019

**Abstract:** Pore network models are powerful tools to simulate invasion and transport processes in porous media. They are widely applied in the field of geology and the drying of porous media, and have recently also received attention in fuel cell applications. Here we want to describe and discuss how pore network models can be used as a prescriptive tool for future water electrolysis technologies. In detail, we suggest in a first approach a pore network model of drainage for the prediction of the oxygen and water invasion process inside the anodic porous transport layer at high current densities. We neglect wetting liquid films and show that, in this situation, numerous isolated liquid clusters develop when oxygen invades the pore network. In the simulation with narrow pore size distribution, the volumetric ratio of the liquid transporting clusters connected between the catalyst layer and the water supply channel is only around 3% of the total liquid volume contained inside the pore network at the moment when the water supply route through the pore network is interrupted; whereas around 40% of the volume is occupied by the continuous gas phase. The majority of liquid clusters are disconnected from the water supply routes through the pore network if liquid films along the walls of the porous transport layer are disregarded. Moreover, these clusters hinder the countercurrent oxygen transport. A higher ratio of liquid transporting clusters was obtained for greater pore size distribution. Based on the results of pore network drainage simulations, we sketch a new route for the extraction of transport parameters from Monte Carlo simulations, incorporating pore scale flow computations and Darcy flow.

**Keywords:** pore network model; Monte Carlo simulation; drainage invasion; porous transport layer; clustering effect; water electrolysis

### **1. Introduction**

Pore network models (PNMs) are discrete mathematical models that are basically used to simulate and predict pore scale processes. Different types of PNMs are generally distinguished. These are (i) PNMs for quasi-static drainage invasion processes [1–4], (ii) PNMs for quasi-static imbibition invasion processes [1,3–5], (iii) PNMs of drainage with phase transition and diffusion of the vapor (especially applied in drying research) [6–8], and (iv) PNMs for the computation of dynamic pore scale fluid flow [9–13]. While the first three approaches usually assume quasi-static invasion of the pore space, the fourth approach considers viscous flow of the liquid phase and dynamic invasion of the pore space. There are several examples that combine (i) and (ii) [3,4,14]. There also exists a number of PNMs that additionally incorporate pore scale mechanisms, such as liquid film flow [15–17] or crystallization [18–20]. Only a few models are available that take into account coupled heat and mass transfer [21,22] or that at least consider the invasion and transport processes under non-isothermal conditions [23–25].

A wide range of PNM applications related to fuel cells is already available in literature, e.g., [26–33]. Major focus of the PN studies related to fuel cells is essentially on liquid water distribution in

dependence of the pore structure [28], liquid clustering [29], the thickness of the gas diffusion layer [30], hydrophobicity [27], local variation of wetting properties, condensation-induced liquid water formation inside the pore structure [31,33], ice formation, model validation by X-ray tomography [32], and neutron tomography [34,35].

So far, PNMs are only rarely applied in water electrolysis research. Experimental studies incorporating microfluidic platforms of pore networks (PNs) are presented by Bazylak et al. [36–39]. Their investigations are based on the correlation of the oxygen flow rate with current density. The microfluidic platforms employed in their studies are based on the geometric information obtained from micro tomography measurements of different titanium porous transport layers (PTLs). They showed that the developing gas bubbles penetrate the porous structure of the model PTLs and form continuous fractal gas branches that cover the PN at the breakthrough point. The network covering gas fingers at the breakthrough point are assumed as stable as long as the volume flow rates correlating with the current density are constant (In [38,39], the air flow rates were 1, 5, and 10 μL min−<sup>1</sup> with current densities of 1.4, 7.0, and 14 A cm<sup>−</sup>2, and the liquid flow rates were 5 and 10 μL min<sup>−</sup>1. The same liquid flow rates and air flow rate of 1 μL min−<sup>1</sup> were applied in the numerical simulations in [37]. In [36], higher liquid flow rates of 505 μL min−<sup>1</sup> and gas flow rates of 10 to 300 μL min−<sup>1</sup> were applied.) It was illustrated that the shape of the penetrating gas fingers and the final saturation of the microfluidic PN with oxygen are dictated by the pore structure when the invasion occurs in a capillary regime. Higher gas saturation values were obtained for lower porosity of the PN and smaller pore and throat sizes (sintered PTL in [39]). However, an explanation for this outcome is not given in the paper. Additionally, only 2-dimensional (2D) microfluidic experiments are presented; as will be discussed in this paper, the relationship between pore size distribution and saturation is different in 3-dimensional (3D) PNs where essentially the interconnectivity of the pore space is higher. However, it is surmised from PN modeling that the pore size distribution (PSD), especially the standard deviation of pore and throat sizes, were different in the three cases studied in [39] (felt, sintered, and foam PTL). The presence of large and small pores results in the competitive invasion of the PTL. This can especially be illustrated for strongly heterogeneous pore structures, e.g., bimodal PSD, [40]. Principally, the invasion process follows the path of least resistance. In the hydrophilic micromodels used in experiments in [39], this is the path that follows the lowest invasion pressure thresholds of the neighboring pores; it is thus the route of the largest neighbor pores. In more detail, invasion of the PTL at constant current densities is a process of quasi-static drainage of water. This process occurs in distinct steps. This is designated as the mechanism 'one throat at a time' in [36] and can be simulated with PNMs. To invade any liquid filled pore or throat with water, the pressure inside the gas phase has to overcome a critical invasion pressure threshold that depends on the curvature of the gas-liquid interface. If the pore and throat sizes are distinct, the invasion events occur at distinct pressures. According to this, the pressure curve follows a trend of alternating pressure increase (to achieve the critical invasion pressure thresholds), during which the saturation of the microfluidic network remains constant, and the pressure drops during the sudden invasion (saturation decrease) when the critical invasion pressure is achieved [41]. In [36], it is also shown that the invasion of the PN can continue after the breakthrough of one gas branch if the oxygen flow rate, and thus the gas pressures are increased accordingly.

In this paper we focus on water electrolysis cells with in the sandwich coordinate countercurrent flow of water and oxygen (perpendicular to the catalyst layer and the proton exchange membrane) [42]. Additionally, we consider the situation of high oxygen production rates (current density = 0.6–6 A/cm<sup>2</sup> [43]; oxygen flux < 4 μL s−<sup>1</sup> mm<sup>−</sup>2). Based on this, we assume plug flow of oxygen and water. In the first approach presented in this paper, we use pore network Monte Carlo simulations (PNMCSs) for the prediction of the invasion of an initially fully water saturated PN with oxygen. We furthermore assume that a stationary distribution of oxygen and water inside the PTL is achieved when the oxygen flow paths connect the catalyst layer, where the oxygen is produced, with the water supply channel, while the water connectivity between both sides is maintained, presuming constant oxygen production rates and constant pressure and temperature along the PTL. As a result of the PNMCSs, we

show and discuss the invasion patterns in the moment when the water supply route is disconnected, as well as the impact of isolated single liquid clusters on the relative permeability of oxygen and water. Following parameter estimation concepts, e.g., in drainage of soils [1] or drying [44,45], we propose to estimate relative permeability from the stationary final gas-liquid distribution obtained from the PN drainage simulation.

### **2. Pore Network Model**

A schematic sketch of the PNM is given in Figure 1. Here we consider idealized 3D PNs with PSDs in the range of typical PTLs usually applied for water electrolysis (Figure 2). The application of idealized lattice structures is commonly practiced if the physical mechanisms and pore scale effects are analyzed. More advanced studies base the PN simulations on the real structure of the porous medium (e.g., [36,46]). This is also foreseen in our studies in the next step. However, in this paper we present and discuss the basics of the method and the basic outcomes of the drainage simulation and thus disregard the more expensive (concerning time and computational effort) method.

**Figure 1.** (**a**) Schematic illustration of the pore network model (PNM) of drainage with pore and throat numbering (pore numbers in blue; throat numbers in black). (**b**) Breakthrough path of oxygen from the catalyst side to the water supply channel. Note the periodic boundary conditions that allow connection of edge throats and pores with each other.

As can be seen in Figure 1a, the PN consists of circular pores, i.e., the sites and cylindrical throats, i.e., the bonds. Both pores and throats can in principal be occupied by either liquid water or gaseous oxygen. The liquid occupied pores are shown in black in Figure 1a; the gas occupied pores are shown in red. Similarly, the liquid occupied throats are shown in blue and the empty throats are shown in red. The pores at the bottom side of the network can be interpreted as water sources as they represent the connection links to the water supply channel. All of these pores are initially saturated with water (black pores in Figure 1). In contrast, pores at the top side of the PN are connected to the oxygen sources inside the catalyst layer at the anode side of the water electrolysis cell and directly supplied with gaseous oxygen (red pores in Figure 1). The water supply route from the water supply channel connected to the bottom of the PN towards the catalyst layer (at the top side) is given by the interconnected blue throats, while the red pores and throats provide distinct routes for the gas phase. Figure 1b shows the path of oxygen through the PN. As can be seen, the oxygen path covers the PN from the catalyst side to the water supply channel in this example. This situation is called a breakthrough. The breakthrough of the gas phase is achieved when one of the bottom pores is occupied with oxygen. On the other hand, the

water supply is interrupted once the blue cluster is either completely split up into numerous smaller single clusters or disconnected from either side.

The geometric arrangement of pores and throats and the neighbor relations are specified by the different matrices pnp, tnt, tnp, and pnt. Following the example in Figure 1a, the neighbor relations of the 3D PN as illustrated here are:

$$\text{pnp}(\text{pore}10) = \begin{bmatrix} \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & 11 & 12 & 13 & 16 & 19\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \end{bmatrix} \tag{1}$$

$$\text{trp}(\text{throat1}) = \begin{bmatrix} 10 & 11 \\ \vdots & \vdots \end{bmatrix} \tag{2}$$

$$\text{pnt}(\text{pore}10) = \begin{bmatrix} \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & 3 & 28 & 34 & 55 & 64\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \end{bmatrix} \tag{3}$$

$$\text{trut}(\text{throat1}) = \begin{bmatrix} 2 & 3 & 28 & 29 & 34 & 35 & 55 & 56 & 64 & 65 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \end{bmatrix} \tag{4}$$

From this, it follows that the coordination number of the PN is 6. Note that we apply periodic boundary conditions. This means that the pores at all lateral sides are connected to each other in order to reduce confinement effects. While the radius distribution of the throats and pores stochastically vary in the PNMCSs, the other geometrical parameter, such as throat length and the neighbor relations, are kept constant (Table 1). The PSDs of pores and throats in the PN simulations presented below are illustrated in Figure 2.

**Figure 2.** Representative pore size distribution (PSD) of the pore networks (PNs) studied in the pore network Monte Carlo simulations (PNMCSs) in Section 3.1. The PN with low standard deviation (STD) of pore and throat sizes is denoted STD 0.5 μm, and the PN with a higher standard variation is denoted STD 1 μm (refer to discussions below in Section 3.1). The first peaks correspond to the sizes of throats and the second peaks to the pore sizes, respectively. Note that the overlap of pore and throat sizes is greater for a higher standard deviation. The porosity of both PNs is 21%.


**Table 1.** Simulation parameters.

We assume that the oxygen is homogeneously distributed along the surface of the PN and neglect spatial and temporal pressure fluctuations. Instead, we assume a constant oxygen supply rate. Based on this, we anticipate plug flow of oxygen and water and compute the invasion of the PN based on the Young–Laplace equation, following the concepts in [8,47,48]:

$$P\_{\rm l,t,p} = P - \frac{2\sigma \cos \theta}{\mathbf{r}\_{\rm t,p}} \tag{5}$$

The order of invasion thus follows the radius distribution of the throats (with index t) and pores (with index p). Viscous flow is disregarded for the drainage invasion computation. The liquid clusters are labeled based on the throat saturation, taking into account the saturation in neighboring pores. This means that any neighbor throats that contain water belong to the same liquid cluster if the pore between them is also saturated with water (Figure 1). The gas phase is not labeled as it is continuous. Note that liquid clusters that are disconnected from the cluster spanning the PTL from the catalyst side to the water channel side (all blue throats in Figure 1a) are not further invaded. These clusters remain in their original size and can be regarded as transport barriers for oxygen flow (Figure 3). In regard to an efficient PTL mass transfer, it would be affordable to avoid clustering of the liquid phase; relevant issues are also intensively studied e.g., in hydrodynamics of porous geological structures and soils [49]. However, in regard to the optimization of the PTL, several aspects interact with each other, including mass transfer, heat transfer, and electrical conductivity [50].

The PNMCSs presented here are restricted to the computation of the point when the fragmented clusters are not further invaded. In most cases, the fragmented clusters do not span the network as they are either connected to only one of its open sides (top or bottom) or isolated in the center of the PN. The simulation can result in several gas branches penetrating the network from top to bottom. (Note that the gas phase always forms a continuum). This is in good agreement with experimental findings reported in [36].

isolated clusters (**a**) (**b**)

continuous liquid cluster

**Figure 3.** Liquid clustering in a 3-dimensional (3D) PN on the example of a small 10 × 10 × 10 network with homogeneous PSD and the parameters specified in Table 1. The transport barrier clusters (isolated and discrete) are shown in grey (**a**,**b**); the liquid transporting clusters are shown in blue (**a**,**b**); and the continuous gas phase is shown in magenta (**a**). Pores are not shown for reasons of readability.

The computation of the quasi-static invasion profile does not incorporate the solution of mass transfer equations because invasion percolation with trapping is assumed here. Instead, the pore scale fluid transport equations are set-up and solved in the second step based on the stationary invasion patterns from step 1 (Figure 4). For the computation in step 2, the breakthrough invasion patterns, i.e., at disconnection of water transport routes, are used [29,30,37–39]. Only the liquid cluster spanning the PN from top to bottom in the moment before disconnection and the continuous gas phase are considered. The mass transfer through the spanning clusters and the gas phase is computed pore-to-pore, based on the Hagen–Poiseuille equation:

$$\dot{\mathbf{M}}\_{\rm l,g} = \frac{\rho\_{\rm l,g} \pi \mathbf{r}\_{\rm l}{}^{4}}{8 \eta\_{\rm l,g} \rm L} \Big(\mathbf{P}\_{\rm l,g,1} - \mathbf{P}\_{\rm l,g,2}\big) . \tag{6}$$

for the liquid phase l and the gas phase g. Incompressible, Newtonian viscous flow is assumed (the compression factor of oxygen at operation conditions of P = 10 bar and T = 50 ◦C is roughly 1). Note that the assumption of the Hagen–Poiseuille flow is a strong model simplification for the pore flow because L d. A more advanced approach would account for the radial velocity of the flow [51]. In Equation (6), Pl,1 and Pg,1 are the liquid and vapor pressures in neighbor pore 1 of a throat (compare with Equation (2)) and Pl,2 and Pg,2 are the liquid and vapor pressures in neighbor pore 2 of that throat, respectively. The resulting set of linear equations is then transferred into the matrix notation:

$$P\_{\rm l,g} = \mathbf{A}\_{\rm l,g} \Big(\mathbf{g}\_{\rm l,g}\big) \big| \mathbf{b}\_{\rm l,g} \tag{7}$$

In these equations, **A**<sup>l</sup> and **A**<sup>g</sup> represent the matrices of liquid and vapor conductivities of the throats (e.g., refer to [52] for further details), with

$$\mathbf{g}\_{\rm l,g} = \frac{\rho\_{\rm l,g} \pi \mathbf{r}\_{\rm l}^4}{8 \eta\_{\rm l,g} \rm L}. \tag{8}$$

Following discussions in [36], conductivities inside pores are not computed as the pores are not interpreted as hydraulic conductors. Validation of this assumption could be further studied by Lattice–Boltzmann simulation, e.g., [53,54].

The boundary conditions for each throat are given in bl and bg, which are specified for the pore neighbors 1 and 2 of a throat. Generally:

$$\mathbf{b}\_{\mathbf{l},\mathfrak{g}} = \mathbf{g}\_{\mathbf{l},\mathfrak{g}} \mathbf{P}\_{\mathbf{l},\mathfrak{g}} \tag{9}$$

The boundary conditions for the PN simulation are P1 at the bottom side of the PN and P2 at the top side (Figure 4). Solving Equation (7) yields the vapor and liquid pressure fields in the PN. With this, Equation (6) can be solved for each throat. Once the liquid and vapor flow rates through the liquid throats in spanning liquid clusters and the gas throats are known, the overall mass flow rates are computed in step 3 (Figure 4c). From this follows the relative permeability for either the liquid (l) or gas phase (g): .

$$\mathbf{K} \cdot \mathbf{k}\_{\mathrm{l,g}} = \frac{\eta\_{\mathrm{l,g}}}{\varrho\_{\mathrm{l,g}}} \frac{\mathbf{M}\_{\mathrm{l,g}}}{\mathrm{A}} \frac{\Delta \mathbf{z}}{\Delta \mathbf{P}\_{\mathrm{l,g}}}.\tag{10}$$

The absolute permeability K is obtained for the same computations but a totally empty (gas permeability) or totally saturated (liquid permeability) network.

**Figure 4.** Extraction of efficient transport parameters based on PNMCSs. (**a**) Step 1: Computation of the steady state invasion patterns at the disconnection of the water supply route by PNMCS. (**b**) Step 2: Computation of the pore scale fluid flow based on the patterns obtained from step 1. (**c**) From the flow rates in step 2, the relative and absolute permeabilities of the PN are computed in step 3 on the Darcy scale.

Note that high pressures up to P = 30 bar are postulated as usual operating conditions of water electrolysis cells. It is remarked that the evaporation of water, even at prevailing temperature of 50–80 ◦C, can be disregarded at such high pressures. Additionally, in the simulations presented here, we generally assumed hydrophilic conditions with cosθ = 1 and constant temperature and pressure as well.

### **3. PN Simulations and Results**

### *3.1. 3D PNMCS of Drainage*

We present here the simulation of one realization of a 3D PN with square lattice and 30 × 30 × 10 pores (Figure 5). The simulation parameters are summarized in Table 1. The overall pressure, as well as the temperature, were kept constant. The lattice spacing between the pores was L = 50 μm and the throat length was kept constant with Lt = 27 μm. The pore and throat sizes varied in the range given in Figure 2, with standard deviations 0.5 μm and 1 μm. The simulations were repeated 20 times with randomly distributed pore and throat radii within the range given in Figure 2.

According to [55] or [36,46], porosity of the PTL usually varies between 54% and 85%, depending on the kind of the material. Usually, sintered powder has a lower porosity than felt PTLs or foam PTLs. However, the porosity in our simulations was only around 21%. Larger porosities would only be achieved in our PN with greater variation of the pore and throat sizes (i.e., by higher standard deviation) and much shorter throat lengths. To reach the high porosities given in the literature, overlapping of the pore volumes and negligible throat length would have been necessary. This would lead to different invasion effects than currently underlying in the proposed PNM, namely site invasion instead of bond invasion and different flow regimes than those postulated above. (Exemplarily, the porosity could be increased to 43% for throat lengths of 7 μm. Lt < rt would then require a revision for the assumption of a developed Poiseuille flow inside the throats). Instead, we expect to achieve higher porosities in PNs based on the real porous structure, with various heterogeneities of the pore space (see discussions in Section 4 below). Additionally, due to the relatively long throat length applied in our PN, the thickness is greater than the given values in literature [36,55] where the 10 pore rows correspond to a thickness between 170 μm and 300 μm. Independent of this, the PN simulations presented below nicely illustrate the invasion mechanisms that drive the drainage process. The presented method can be easily adapted to the simulation of the real porous structure.

**Figure 5.** 3D PN under study. Pores are shown in black and throats in blue. The top and bottom pore rows are only vertically connected to their neighbor throats. All pores and throats are initially saturated with liquid.

In these simulations, pores were invaded independently of their throat neighbors. Note that the drainage process is principally a bond invasion process wherefore pores could be invaded together with the throat neighbors [4]. This is also revealed by the experiments presented in [36]. However, as shown in Figure 2 and also represented in Figure 6, the pores are comparably large in our PN with invasion pressure thresholds in the range of the throats. The curves of pore and throat capillary pressures partly overlap and thus indicate the competitive invasion. Apart from that, the invasion pressure thresholds vary linearly with the size of pores and throats if all other parameter are kept constant (Equation (5), Table 1). Figure 6 indicates that the capillary pressures are rather small in the given range of pore sizes, which results in the decrease of liquid pressure (against the gas pressure) of only a few millibars and thus a concave gas-liquid interface.

**Figure 6.** Capillary pressure variation in the studied range of pore sizes (standard deviation of 0.5 μm, Figure 2). Due to the overlapping of the curves the pore invasion is computed independently of the throat invasion.

The simulation results of one PN simulation (with standard deviation 0.5 μm, Figure 2) are summarized below and in Figures 7–9. At first, the capillary pressure curves of the PN simulation are shown in Figure 7. In the hydrophilic drainage process studied here, the available meniscus pore or throat with the lowest capillary pressure is first invaded. Since all surface pores are available with liquid menisci at the start of the simulation, the blue curve initially follows the radius distribution of the surface pores. Once all surface pores are invaded, the capillary curve passes into the trend of the throats, whereas the pore invasion becomes random depending on the distribution of the available meniscus pores. Note that more and more pores become available for invasion when the drainage front widens. However, the blue cloud only represents the pores in liquid clusters spanning the PN, since all other clusters are stationary and not invaded. This holds for the throats analogously. The black curve follows the invasion of the largest throats in the spanning clusters, while the randomly distributed points below the curve correspond again to the PSD of throats in the PN.

**Figure 7.** Capillary pressure curves of the drainage invasion process initiated in all surface pores. Note that the small pores and throats with higher capillary pressures (Figure 6) are not invaded.

Figure 8 summarizes the gas-liquid distributions of the PN simulation. At first, the situation a few moments before loss of the liquid connectivity between the vertical throats in the top and bottom row is illustrated in Figure 8a. As can be seen, the liquid phase is multiply disconnected and penetrated by gas branches (in white). The gas phase forms a fractal structure that moves downwards towards the water supply channel. The invasion of the gas phase leads to the formation of isolated liquid clusters

that remain behind the front. The size of these clusters is stable as they are not further invaded. It is remarked that the invasion process occurs in all three dimensions in the homogeneous pore structure underlying this simulation (isotropic invasion). This is explained with the equal distribution of liquid pressure (associated with PSDs) throughout the PN, constant wettability, temperature, and pressure conditions. However, it is noted that the horizontal breakthrough occurs much earlier than the vertical breakthrough. This is explained with the oxygen accessibility of the surface pores. In detail, all of the surface pores are initially connected to one liquid-filled throat (downwards) and the gas bulk phase (upwards). Essentially, based on the invasion percolation algorithm and due to the fact that the peak of the PSD of pores is shifted towards larger radii, all pores empty initially before any throat is invaded. This results in the formation of an oxygen front that spans the PN horizontally (the size is roughly NixNj) (Figure 8c,d). According to this, the liquid connectivity between the boundary pores at the top side and the bottom side is already lost at the start of the process. A different situation would be expected if the surface pores were smaller and the invasion pressure were higher.

**Figure 8.** (**a**) PN shortly before loss of connectivity of the liquid supply route through the blue cluster. (**b**) Only the liquid transporting cluster (main cluster) is shown here. Note that due to the applied periodic boundary conditions, the cluster appears at two sides of the PN, while it is disconnected in the center. (**c**) Empty surface pores at the top side (catalyst side). Pores connected to a liquid-filled throat of the spanning cluster are shown in light blue. (**d**) Surface pores at the bottom side (water channel side). Pores connected to the spanning (main) cluster in blue; all other pores in red.

At the moment when the liquid supply route is interrupted (or shortly before that moment), the overall saturation with liquid is S = 0.5996. However, as shown in Figure 8, there exists only one liquid transporting cluster (main cluster). The other liquid-filled throats are shown in grey and are either single throats or part of isolated clusters. The number of clusters at the interruption of liquid connectivity between the surface throats (first and last row of vertical throats) is 933. The maximum cluster size is 624 throats; this is the blue cluster in Figure 8a,b. The mean size of the other clusters is 9.5 throats, thus approximately two orders of magnitude smaller. The liquid transporting cluster contains approximately 5.7417 <sup>×</sup> 10−<sup>3</sup> <sup>μ</sup>L of water. Its volume referred to the total liquid volume contained in the PN at the interruption of the liquid path is only 3%. This reveals the relevance of the pore scale information about the liquid connectivity for the calculation of the liquid permeability. A different situation is expected in presence of wetting liquid films, as will be discussed below. It is also noted that higher permeabilities would be obtained if the drainage simulation would be interrupted already at a higher overall saturation. However, in the simulation presented in Figures 7–9,

a breakthrough of the gas phase occurred shortly before disconnection of the main cluster from the top side.

Furthermore, the saturation of the exchange interfaces plays a vital role. The surface and bottom saturations are highlighted in Figure 8c,d. Note that the light blue pores at the top side are already empty. These pores are connected to liquid-filled vertical throats right below the surface. The bottom pores shown in blue are liquid filled. The red pores in these images are either empty (filled with oxygen) or belong to isolated clusters. The ratio of the wetted surface (taking into account the vertical throats inside the first and last layer of the PN) is around 2% at the top surface and around 1% at the bottom side, thus very low. This results in low liquid permeabilities, as will be shown below.

From the distributions shown in Figure 8, the overall liquid saturation Sl of the pore network can be calculated from:

$$\mathbf{S}\_{\mathbf{l}} = \frac{\sum\_{\mathbf{l}=1}^{22500} \mathbf{V}\_{\mathbf{l}} \cdot \mathbf{S}\_{\mathbf{l}} + \sum\_{\mathbf{N}\_{\mathbf{P}}=1}^{9000} \mathbf{V}\_{\mathbf{P}} \cdot \mathbf{S}\_{\mathbf{P}}}{\sum \mathbf{V}\_{\mathbf{t}} + \sum \mathbf{V}\_{\mathbf{P}}},\tag{11}$$

with pore and throat saturation Sp and St. The overall gas saturation Sg is:

$$\mathbf{S\_{\%}} = \mathbf{1} - \mathbf{S\_{\text{lv}}} \tag{12}$$

accordingly. The liquid and gas saturations are calculated for each network slice (Nk = 1:10) and illustrated as a function of the vertical position ((Nk − 1)·L, see Figure 1) in Figure 9a,b. Figure 9a,b shows the transient saturation profiles of one PN simulation. It reveals that the PN is basically invaded by oxygen from the top side. If one disregards the top and bottom surface layers of the PN (as they have a different overall volume, Figure 1), the saturation of liquid is evenly distributed. This means that isolated liquid clusters homogeneously cover the PN. The liquid confined in the liquid transporting main cluster spanning the PN is illustrated in Figure 9c. The liquid saturation of the main cluster is calculated analog to Equation (11) by only taking into account the liquid inside it. Two different options are compared with the liquid saturation profiles in Figure 9c. These are option 1:

$$\mathbf{S\_{MC,V\_{tot}}} = \frac{\sum\_{\mathbf{N\_t(N\_k)}} \mathbf{V\_t \cdot S\_{t,MC}} + \sum\_{\mathbf{N\_P(N\_k)}} \mathbf{V\_P \cdot S\_{p,MC}}}{\sum\_{\mathbf{N\_t(N\_k)}} \mathbf{V\_t + \sum\_{N\_P(N\_k)}} \mathbf{V\_P}},\tag{13}$$

saturation on the base of the total void volume contained in slice Nk; and option 2:

$$\mathbf{S\_{MC,V\_l}} = \frac{\sum\_{\mathbf{N\_t(N\_k)}} \mathbf{V\_t \cdot S\_{t,MC}} + \sum\_{\mathbf{N\_t(N\_k)}} \mathbf{V\_P \cdot S\_{p,MC}}}{\sum\_{\mathbf{N\_t(N\_k)}} \mathbf{V\_t \cdot S\_t} + \sum\_{\mathbf{N\_p(N\_k)}} \mathbf{V\_P \cdot S\_P}},\tag{14}$$

saturation on the base of the total liquid volume contained in that slice. Note that only the profiles for S = 1 till S = 0.62 are shown here, because the main cluster splits into two clusters at S = 0.62. The agreement of the curves in Figure 9a,c is good at the start of drainage (when the overall saturation S is high), indicating that initially nearly all the liquid is contained in one cluster. At later stages of the invasion process, a gradient develops. Comparison of the blue and black lines in Figure 9c with the red curves reveals that the center of the volume of the main cluster is located at the bottom side of the PN, thus closer to the water supply channel, whereas connectivity to the top side is lost already at the very beginning of the drainage (when all surface pores are invaded). The difference between the black curves and the red curves is related to the liquid volume contained in single clusters. As can be seen, more single clusters are found inside the upper region of the pore network (i.e., at higher values of (Nk − 1)·L). This reveals that the main cluster is traveling through the PN, leaving the isolated clusters behind its front. From this it can be expected that the relative permeability of the

liquid phase continuously decreases with progress of invasion, while the oxygen permeability is 0 as long as the breakthrough of the gas phase is not observed (Section 3.2). This outcome might be an indicator for the relevance of the thickness of the PTL for water exchange processes in water electrolysis cells [56]. According to this, thinner PTLs might be more convenient in terms of the water supply routes. However, further simulations are necessary to proof the consistency of this assumption. This finding is essentially an explanation for the low permeability of the liquid phase, as discussed below.

**Figure 9.** (**a**) Transient liquid saturation profiles of the one PN simulation shown in Figures 7 and 8. (**b**) The according gas saturation profiles. The catalyst side is found on the right (corresponding to the top of PN) and the water supply channel is found on the left (corresponding to the bottom of PN). The saturation profiles are shown for overall liquid saturation S = [1:0.61] (the final saturation is S = 0.6). (**c**) Saturation profiles of the liquid transporting clusters for S = [1:0.62].

It is remarked that different simulation results are obtained for a higher standard deviation (1 μm) (Figure 10). The liquid saturation is still homogeneously distributed throughout the network, but the level of the final slice saturation is lower, namely S = 0.575 (in Figure 9 S = 0.6). Moreover, the main liquid cluster is larger and appears more dense. It contains 8% of the total liquid volume. This reflects the importance of studying the impact of the PSD rather than the impact of porosity, which is around 21% in both PNs [57].

**Figure 10.** (**a**) Transient saturation profiles of PN with a standard deviation of 1 μm. The saturation profiles are shown for overall liquid saturation S = [1:0.58] (the final saturation is S = 0.575). (**b**) Main cluster.

Consistency of these results is proven by each 20 repetitions of the PN simulation (PNMCS). The results are illustrated in Figure 11, which shows the steady state saturation profiles at the moment of interruption of liquid transport. It is highlighted that the average saturation in the center of the PN is slightly higher in the network with lower standard deviation. This means that more liquid is contained in isolated clusters in this case, as previously observed.

**Figure 11.** (**a**) Steady state saturation profiles of 20 PN realizations with standard deviation 0.5 μm. (**b**) Steady state saturation profiles of 20 PN realizations with standard deviation 1 μm.

Similar simulation results are obtained if the PN surface is invaded at only one sight. In the simulation shown in Figure 12, the one open pore at the surface is highlighted in light blue. Note that the PN is identical to the one shown in Figures 7–9. The invasion follows a similar route as before and the main cluster spanning the PN at breakthrough is almost identical. The saturation profiles in Figure 12 are obviously different because the surface pore row basically remains saturated. However, the average saturation in the center of the PN is again 0.6 at the end of the process. Isotropy of the invasion process can be shown for such a situation, i.e., when only one pore in the center of the surface area is accessible for oxygen (Figure 13).

**Figure 12.** (**a**) PN invasion from one surface pore (shown is the complete PN with gas phase in white, disconnected clusters in grey, and the liquid transporting cluster at the breakthrough of the gas phase in blue). (**b**) Saturation profiles for S = [1:0.64].

Figure 13a plots the label of that in the invasion event invaded vertical (Nk-direction) or horizontal (Ni and Nj-direction) throat (also refer to Figure 1). The three clouds at the start of the process indicate that throats are invaded in each direction (red circles). The breakthrough point is achieved when the throat with the lowest label in each horde is invaded. This is highlighted in Figure 13a. The plot furthermore reveals that, afterwards, horizontal and vertical throats are equally invaded. Figure 13b shows the invasion velocity:

$$\frac{\text{dI}'}{\text{dI}\_{\text{tot}}} = \frac{\text{N}\_{\text{ver,hor}}}{\text{N}\_{\text{tot}}},\tag{15}$$

i.e., the number of invasions in the horizontal and vertical throats (Nhor and Nver, respectively) I' related to the total number of invasions (pores and throats) Itot. According to Figure 13b, more vertical throats are initially invaded, whereas, at later stages of the drainage process, the invasion occurs more often in horizontal throats. This is explained with the PN geometry in Figure 5. It is the main reason for the clustering effect discussed above. In order to prevent the clustering effect, it would be more convenient if the invasion process was anisotropic and occurred mainly in vertical throats. Taking the invasion pressure curve in Figure 6 as a reference, this could be achieved by a gradient in throat sizes (small at the top side, large at the bottom side).

**Figure 13.** (**a**) Order of invasion on the example of the PN discussed above. According to Figure 1, higher throat labels are associated with the top side of the PN. (**b**) Invasion velocities in vertical and horizontal direction.

It is highlighted that a different invasion behavior would occur if also temperature, pressure, or wettability would spatially vary along the vertical or the horizontal direction. Depending on the situation, anisotropic invasion can also be expected. We plan to illustrate such situations in a forthcoming paper.

### *3.2. Estimation of Relative Permeabilities*

The final gas-liquid phase distribution of the PN drainage simulation is employed for the calculation of gas and liquid permeability. Note that due to the trapping of clusters, the relative permeabilities of the liquid and gas phases cannot be related to each other (kl 1 − kg).

For calculation of the absolute and relative permeabilities, we have considered the liquid flow through pores in Nk = 2 and Nk = 9 (Figures 1 and 5) because all surface pores are initially invaded. The pressure gradient is imposed between the pores in these rows. The relative permeability of the liquid phase is related to the total surface area, although only 1–2% of the interfaces are wetted by the liquid spanning cluster. For the PN simulation presented in Figures 7–9, we found an absolute permeability of K = 1.9639 <sup>×</sup> 10−<sup>12</sup> m2 (liquid phase) and K = 6.5612 <sup>×</sup> 10−<sup>12</sup> m<sup>2</sup> (gas phase) (in [50], similar values are obtained for a porosity of 50%). The relative permeabilities are kl = 7.8439 <sup>×</sup> 10−<sup>4</sup> for the liquid phase and kg = 0.0068 for the gas phase. The relatively low oxygen permeability is associated to the hindrance by isolated liquid clusters. Higher liquid permeabilities are obtained if the calculations are repeated for higher liquid saturations of the PN, i.e., if the PN is not fully

drained (Figure 14). (Note that the gas phase is not penetrating the PN for liquid saturations shown in Figure 14). These values are in good agreement with previous results in [45] and they clearly illustrate the different relative permeabilities of 2D and 3D PNs, which are related to the higher connectivity of the 3D PN and the associated clustering effect. This situation can change if wetting liquid films are considered in the PNM [57]; however, their development and extension depend on the wettability of the surface, temperature, pressure, and process conditions. (In [57] it is mentioned that the wettability alteration of titanium PTLs occurs due to the oxidation of its surface.) Additionally, a heterogeneous pore structure with interpenetrating large and small pores can lead to higher relative permeabilities if the breakthrough occurs in the large pores, whereas the small pores remain saturated and provide the pathway for water transport. The relative permeabilities are thus a function of PSD [57]. The absolute permeability depends on the ratio of pore space to solid, i.e., on the porosity (e.g., [1,50]). However, according to Equation (10), higher absolute permeability does not necessarily increase relative permeabilities.

**Figure 14.** Permeability variation for simulations with low and high standard deviation of PSD.

### *3.3. Oxygen Production Rate*

The decomposition of water occurs if a voltage higher than the thermoneutral voltage is supplied to the electrolysis cell, e.g., U = 1.5–2.2 V [43]. It is anticipated that, at low voltages, oxygen can be solved inside the water that is transported through the PTL towards the catalyst layer. However, if higher oxygen production rates are achieved, the gas can invade the PTL if the pressure in the gas phase is high enough. In [58,59], the different transport regimes for oxygen through the PTL that correlate with the current density of the electrolyzer are shown. In this, dispersed bubbly flow (#1), plug flow (#2), slug flow (#3), churn flow (#4), and annular flow (#5) are distinguished. The latter flow regime is based on wetting liquid films forming along the surface of the pores and sustaining liquid transport between the water flow channel and catalyst layer. For simplicity, we assumed that the current density is always high enough to allow for a plug flow of the gas phase without formation of liquid films; this is option #6 (Figure 15a,b). This assumption is reasonable in the face of the total volume of the pore network of only a few μL (see below). However, if good wettability of the PTL surface is anticipated (hydrophilic properties), liquid films are likely to form along the solid pore surface. These films can enhance liquid transport through invaded pores as they can connect the water supply channel with liquid clusters at the top side of the PTL [15,16,60]. Basically, in presence of liquid films, the previously (and above discussed) isolated clusters would not occur, but rather the complete liquid phase in the PN would be interconnected and participate in liquid transport. This could enhance liquid transport properties significantly. (Exemplarily, in drying processes, liquid films can reduce the overall drying time by orders of magnitude, e.g., [15,17]). The presence of these films strongly depends on the surface roughness, the geometry of pore corners, wettability, temperature, and pressure, as previously mentioned.

**Figure 15.** (**a**) Plug flow invasion of the gas phase (invasion regime #6) completely separating the gas transporting from the liquid transporting pore space. (**b**) Annular invasion of the gas phase according to [58] (invasion regime #5), leading to annularly wetted pore space transporting gas towards the water supply channel in the center and liquid films transporting water towards the catalyst layer, as well as fully saturated pores and throats.

The drainage algorithm presented above in Section 2 works independently of the current oxygen production rate. We postulate that the oxygen production rate is always high enough to (i) exceed the solubility threshold of oxygen in water (around 40 mg/kg water at 25 ◦C and 1 bar up to 1 g/kg water at the same temperature and around 25 bar [61]) and to (ii) simultaneously flood the drained channels with oxygen. If this is fulfilled, the invasion process can be seen as fast and is expected to take place in a very short time period.

The oxygen production rate can be calculated using Faraday's law, assuming that the current efficiency of oxygen is equal to 1 and stochiometry

$$2\text{H}\_2\text{O} \rightleftharpoons 4\text{H}^+ + 4\text{e}^- + \text{O}\_2.\tag{16}$$

Based on Faraday's law:

$$\mathbf{Q} = \mathbf{I} \cdot \mathbf{t} = \mathbf{F} \cdot \mathbf{z} \cdot \mathbf{N}\_{\mathbf{O}\_{2'}} \tag{17}$$

where z is the number of exchanged electrons (z = 4) and F is Faraday's constant (F <sup>=</sup> 9.64853 <sup>×</sup> <sup>10</sup><sup>4</sup> A s mol<sup>−</sup>1) the moles of electrolyzed water can be correlated with the current of the electrolyzer. From this follows:

$$\frac{d\mathbf{dN\_{O\_2}}}{dt} = \dot{\mathbf{N\_{O\_2}}} = \frac{\mathbf{I}}{\mathbf{F} \cdot \mathbf{z'}}\tag{18}$$

the formed molar amount of oxygen per unit time. The volume flux of oxygen can be further calculated from: .

$$
\dot{\mathbf{v}}\_{\rm O\_2} \left[ \frac{\rm m^3}{\rm m^2 s} \right] = \frac{\rm N\_{\rm O\_2}}{\rm A} \cdot \frac{\overline{\rm R} \cdot \rm T}{\rm P} = \frac{\rm J}{\rm F \cdot z} \cdot \frac{\overline{\rm R} \cdot \rm T}{\rm P}, \tag{19}
$$

with current density

$$\mathbf{J} = \frac{\mathbf{I}}{\mathbf{A}} \tag{20}$$

in A m−2. The relationship in Equation (19) is illustrated in Figure 16a; the linear dependence is in agreement with values reported in [39]. Taking the PN from Section 3.1 as a basis, the surface area connected to the catalyst layer is roughly 0.22 mm2. If it is furthermore assumed that the complete surface area is electrolytically reactive, the oxygen flux in Figure 16a can be converted into a volumetric flow rate: .

$$\dot{\mathbf{V}}\_{\rm O\_2} \left[ \mu \text{L min}^{-1} \right] = 0.22 \left[ \text{mm}^2 \right] \dot{\mathbf{v}}\_{\rm O\_2} \left[ \text{mm min}^{-1} \right] \cdot 10^3. \tag{21}$$

**Figure 16.** (**a**) Linear dependence of the oxygen volume flux on current density. (**b**) Required invasion time for the constant PN volume VPN (slope) in dependence of the reversed volumetric production flow rate.

Figure 16b illustrates the time dependence of the PN invasion on the reversed volumetric flow 1/ . VO2 <sup>=</sup> 1/ . vO2A [sμL<sup>−</sup>1], with the slope being the total volume of the PN in Section 3.1,

$$\mathbf{t} = \mathbf{V}\_{\text{PN}} \frac{1}{\dot{\mathbf{v}}\_{\text{O}\_2} \mathbf{A}} \mathbf{A}' \tag{22}$$

and A being the total cross section of the PN surface. The total volume of the PN discussed in Section 3.1 is approximately VPN = 0.19 μL; the pore volume is 0.0234 μL and the throat volume is 0.17 μL. In detail, at a current density of 2 A cm<sup>−</sup>2, the complete PN could be invaded within 0.0675 s, whereas reduction of the current density by factor 1000 increases the invasion time to roughly 1 min. This is in agreement with experiments reported in [36], where, at postulated current densities of 7 A cm−2*,* the according volume flow rates of the gas phase resulted in invasion times <1 min. Note that always only the dry part of the PN is available for the gas phase. Since this region is much smaller than the total volume, the breakthrough of the gas phase can occur in a shorter time. Additionally, the available cross section for gas invasion is usually lower than the given value in the presence of some liquid pores at the surface (Figure 8).

### **4. Summary and Discussions**

We have presented a method to study the pore scale transport of oxygen and water through the PTL of water electrolysis cells. The method is based on Monte Carlo pore network simulations (PNMCS). It allows simulation of the distribution of the two fluids inside the pore space on a physical base, i.e., without incorporating effective parameters. We have shown that this can be a useful and reliable tool to numerically estimate the pore scale distribution of gas and liquid and the transport coefficients of PTLs, such as relative and absolute permeabilities. More clearly, if these are estimated more accurately based on the proposed method, transport characteristics of PTLs could be predicted more precisely in the future. Moreover, it is highlighted that the discrete method allows correlation of specific transport characteristics with the individual structure of the pore space. The proposed method, therefore, opens up a new route for designing powerful PTLs on customized demand in the future.

The developed PNM basically incorporates invasion percolation rules for hydrophilic drainage with trapping of liquid clusters. Pores and throats are separately invaded in our PN because of their similar invasion pressure thresholds. We assumed plug flow of the liquid and gas phases and negligible film flow. We also assumed that the oxygen production rate is always high enough to homogeneously invade the PN through all surface pores; pressure gradients were disregarded. The correlation with Faraday's law showed that the small PN can be invaded during less than one minute if the current density is higher than 2 mA cm<sup>−</sup>2. The invasion process was simulated until the disconnection of the spanning liquid clusters, i.e., when the transport route for water from the top (assumed to be connected to the catalyst layer) to the bottom (assumed to be connected to water supply channel) was interrupted. Due to the discrete character of the model, it was possible to distinguish between isolated clusters, which do not contribute to liquid transport but rather hinder the oxygen flow, and the liquid transporting clusters, which span the network. The possibility to separate single liquid pores and throats depending on their designation as transporting or non-transporting elements particularly reveals the strength of PN modeling and the perspectives for future applications. Different situations were studied. At first, we illustrated relevance of the pore size distribution (PSD). Secondly, we changed the invasion rules of surface pores. For the latter, we allowed invasion of all surface pores in the first situation and restricted invasion to only one pore in the center of the PN in the second situation. We have shown that the final saturation of the PNs with liquid as well as the liquid volume contained in the spanning liquid cluster depends primarily on the PSD; the porosity was kept constant in the situations studied in this paper. Additionally, our simulations revealed that the liquid transporting clusters cover only a very small percentage of the total volume of the PN, whereas most of the liquid phase is disconnected in isolated single clusters if the drainage simulation is continued until disconnection of the water supply route. Additionally, the spanning cluster was travelling through the PN during the drainage invasion and had its center of mass at the bottom side at the moment when liquid connectivity was interrupted. The clustering effect and the travelling spanning cluster are major outcomes of the PN drainage simulation with a significant impact on the liquid and gas permeabilities, which were very low in the studied cases. Furthermore, we found that the spanning liquid clusters providing the transport route for water cover less than 2% of both surfaces (top, bottom). Though this low value can be attributed to the idealized PN structure, the associated low porosity, and the postulated constant boundary conditions, new PN simulations shall further enlighten the major characteristics of mass transfer in PTLs. Therefore, we plan to incorporate micro tomography scans of the PTL in the PNM in the next step. The instrument to be used is a Procon CT alpha with maximum resolution of 0.6 μm, equipped with image processing software, Mavi, from the Fraunhofer Institute for Industrial Mathematics ITWM Kaiserslautern. Following concepts, e.g., described in [1,2,62], we aim to run the PN simulation using the extracted real structure of the porous medium. This will allow more realistic prediction of the gas-liquid distribution for real PTLs as well as study of the role of heterogeneous pore structures, such as bimodal PSDs. Furthermore, microfluidic experiments will be necessary to validate the assumptions of the model. Current open questions concern the limits of the flow regimes (based on discussion in [58,59]), the relevance of liquid flow through corner films [57], the impact of local temperature and wettability variations, as well as the dynamic invasion in the presence of current density fluctuations. Based on [36], we plan to illustrate in more detail the impact of flow regimes (in terms of higher capillary numbers), as well as unsteady operation conditions related to the fluctuation of the current density. In the latter case, we expect multiple redistribution of liquid due to pressure (and eventually also temperature) variations. This implies the application of imbibition invasion rules additionally to drainage rules. In this context, it will also be worth studying if the disconnected liquid clusters can be reconnected and open up more routes for water transport, especially in the situation where liquid films support liquid flow through corners.

**Author Contributions:** Conceptualization, N.V., E.T., and T.V.-K.; methodology, N.V., E.T., and T.V.-K.; software, N.V. and H.A.; validation, N.V., H.A., and T.V.-K.; formal analysis, N.V. and T.V.-K.; investigation, N.V.; resources, N.V. and E.T.; data curation, N.V. and H.A.; writing—original draft preparation, N.V.; writing—review and editing, T.V.-K.; visualization, N.V.; supervision, N.V.; project administration, N.V., E.T., and T.V-K.; funding acquisition, N.V., E.T., and T.V.-K.

**Funding:** This research received no external funding.

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

#### *Processes* **2019**, *7*, 558

### **Abbreviations**


### *Processes* **2019**, *7*, 558

### **Abbreviations**


### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

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

### **Haashir Altaf 1, Nicole Vorhauer 1,\*, Evangelos Tsotsas <sup>1</sup> and Tanja Vidakovi´c-Koch <sup>2</sup>**


Received: 12 February 2020; Accepted: 19 March 2020; Published: 21 March 2020

**Abstract:** 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 flows counter currently to the educt through the PTL. When it invades the water-filled 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 film flows 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 microfluidic experiments from literature. The invasion patterns of different 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.

**Keywords:** pore network model; drainage invasion; pore size distribution; porous transport layer; electrolysis

### **1. 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–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–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–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/cm2, the hydrogen production rates were very low and at values above 1 A/cm2 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–16,20,26–28] and modeling approaches [29–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<sup>2</sup> 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.


**Table 1.** Study of the interrelation of mass transfer limitations and porous transport layers (PTL) structure.

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–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.

**Figure 1.** (**a**) Schematic representation of different porosity regions in PTL as studied in [33]. (**b**) Pore size gradient within PTL as in the study from [39].

### **2. 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–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–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:

$$P\_l = P\_{atm} - \frac{2\sigma\cos\theta}{r} \tag{1}$$

where *Pl* is the liquid pressure, *Patm* is the atmospheric pressure, σ in kg s−<sup>2</sup> 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.

**Figure 2.** Different throat size distributions with same porosity but different mean throat diameter and standard deviation of the throat size. Solid in white and liquid saturated void space (i.e., pores and throats) in blue.

**Figure 3.** Different saturations for constant porosity and constant mean throat diameter but different organization of the pore network (PN). Solid and empty pores and throats in white and liquid saturated pores and throats in blue.

### **3. 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.

**Figure 4.** Scheme of the algorithm.

**Figure 5.** Pore and throat numbering in a 2D PN. The dashed lines illustrate the periodicity at the lateral boundaries.

### *3.1. Network Generation*

Pore and throat radii (*rp*,*rt*) 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.

**Figure 6.** Geometric information about pores and throats.

#### *3.2. 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.

### *3.3. 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.

### *3.4. Model Assumptions*


### **4. 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.

**Figure 7.** Microfluidic drainage experiment from Arbabi et al. [20] (Reprinted with the permission from Elsevier, 2014). Solids and gas invaded area in black, liquid in white. The PN is identified by the 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.

**Figure 8.** (**a**) PNM simulation result, (**b**) invasion pattern comparison of simulation and experiment result. Liquid-filled throats are shown in black, invaded pores in red, and invaded throats in blue. Liquid-filled pores are not shown.

### **5. Monte Carlo Simulations**

### *5.1. 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.

**Figure 9.** Histograms of pore size distribution (PSD) with varying standard deviation in μm as indicated in the legend.

**Figure 10.** Gas saturation at breakthrough of the gas phase for different PSDs. The mean value of throat sizes is 17 μm.


**Table 3.** Total void volumes and porosities associated with the gas saturation at breakthrough.

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.

### *5.2. 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.** 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.

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.

**(a) Figure 12.** *Cont*.

**Figure 12.** (**a**) Exemplary invasion patterns of the monomodal PN with porosity 71%. Liquid in blue, gas in white and solid in gray. The arrow indicates the direction of gas invasion. (**b**) Saturation profiles for different overall number of invaded throats and pores achieved during one drainage simulation of the PN with randomly distributed pore and throat sizes.

**Figure 13.** (**a**) Exemplary invasion patterns of the bimodal PN with porosity 71%. Liquid in blue, gas in white and solid in gray. Macro-pores are represented by thicker lines. The arrow indicates the direction of gas invasion. (**b**) Saturation profiles for different overall number of invaded throats and pores from one drainage simulation.

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:

$$\text{Ca} = \frac{\Delta P\_I}{\Delta P\_c} \tag{2}$$

with liquid pressure difference Δ*Pl* and capillary pressure gradient Δ*Pc*. 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

$$
\Delta P\_l = \frac{8\eta v L\_t}{\tilde{r}^2} \tag{3}
$$

And

$$
\Delta P\_{\varepsilon} = \frac{2\sigma\sigma\_{0}}{\overline{r}^{2}}\tag{4}
$$

Ref. [58] considering the PSD with mean throat radius *r* and standard deviation σ0. (In Equation (3), *Lt* denotes the length of a throat). The conditions for capillary fingering are fulfilled, if the capillary number is below 10−<sup>4</sup> [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−<sup>3</sup> Pa s and surface tension σ = 0.073 N/m at 20 ◦C, throat length *Lt* = 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*.

### *5.3. 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.

**Figure 14.** PSDs used for different PTL types with standard deviations: 1.0 μm for foam, 1.7 μm for felt and 2.5 μm for sintered.

**Figure 15.** Invasion pattern comparison: (**a**) foam PTL, (**b**) felt PTL, (**c**) sintered PTL (Experimental images from [20] are reprinted with the permission from Elsevier, 2014).

### **6. 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.

**Author Contributions:** Conceptualization, N.V. and T.V.-K.; methodology, N.V.; software, H.A.; validation, H.A.; formal analysis, N.V. and T.V.-K.; investigation, H.A.; resources, N.V. and E.T.; data curation, H.A.; writing—original draft preparation, H.A.; writing—review and editing, N.V.; visualization, H.A.; supervision, N.V. and T.V.-K.; project administration, N.V., T.V.-K. and E.T.; funding acquisition, N.V., T.V.-K. and E.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

### **Symbols**


### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Extreme Learning Machine-Based Model for Solubility Estimation of Hydrocarbon Gases in Electrolyte Solutions**

### **Narjes Nabipour 1, Amir Mosavi 2,3,4,5, Alireza Baghban 6, Shahaboddin Shamshirband 7,8,\* and Imre Felde <sup>9</sup>**


Received: 26 November 2019; Accepted: 7 January 2020; Published: 9 January 2020

**Abstract:** Calculating hydrocarbon components solubility of natural gases is known as one of the important issues for operational works in petroleum and chemical engineering. In this work, a novel solubility estimation tool has been proposed for hydrocarbon gases—including methane, ethane, propane, and butane—in aqueous electrolyte solutions based on extreme learning machine (ELM) algorithm. Comparing the ELM outputs with a comprehensive real databank which has 1175 solubility points yielded R-squared values of 0.985 and 0.987 for training and testing phases respectively. Furthermore, the visual comparison of estimated and actual hydrocarbon solubility led to confirm the ability of proposed solubility model. Additionally, sensitivity analysis has been employed on the input variables of model to identify their impacts on hydrocarbon solubility. Such a comprehensive and reliable study can help engineers and scientists to successfully determine the important thermodynamic properties, which are key factors in optimizing and designing different industrial units such as refineries and petrochemical plants.

**Keywords:** hydrocarbon gases; solubility; natural gas; extreme learning machines; electrolyte solution; prediction model; big data; data science; deep learning; chemical process model; machine learning

### **1. Introduction**

Solubility of hydrocarbon and nonhydrocarbon gases—i.e., mixtures of methane, ethane, propane, CO2, and N2 in aqueous phases—is known as one of the important practical and theoretical challenges in petroleum, geochemical, and chemical engineering. This property has an effective role in different processes, such as achieving optimum conditions for oil and gas transportation, gas hydrate formation, designing thermal separation processes, gas sequestration for protecting environment, and coal gasification. Petroleum reservoirs normally have some natural gases with aqueous solution at

high-pressure and high-temperature conditions so that the solubility of gas becomes attractive for engineers [1–8]. In production and transportation of hydrocarbons, it is possible that water content of gas undergoes an alteration in phase from vapor to ice and gas hydrates. The crystalline solid phases called gas hydrates are created when small-sized gas molecules are trapped in lattice of water molecules. Creation of hydrates can cause major flow assurance problems during production and transportation of hydrocarbons steps such as pipeline blockage, corrosion, and many other issues resulted from the two-phase flow [1,9–11].

In the recent years, investigations on CO2 solubility in aqueous electrolyte solutions have grown significantly as well as they are related to CO2 capture and storage. It is a clear fact that the dominant cause of global warming is emission of CO2 gas generated from fossil fuels so its sequestration and disposal in the ocean have been known as a reasonable choice to overcome global warming problems [12–14]. Simulation of enhanced oil recovery, design of supercritical extraction, and optimization of CO2 dissolution in the ocean need a comprehensive knowledge about carbon dioxide solubility in aqueous electrolytes solutions [13–15].

Investigation of natural gas phase behavior in aqueous solutions in different operational conditions is known one of the important issues in the industry, which has wide applications for avoiding problems in designing and optimization of gas processing. In the literature, there are different solubility datasets for various gas–liquid systems. These datasets mostly include hydrocarbons' dissolution in water/brine systems [1,4,5,9,16–20] and non-hydrocarbons such as CO2 and N2 dissolution in water/brine systems [7,12–14,18,21–24]. A brief summary of the hydrocarbon systems datasets is shown in Table 1 for hydrocarbons. The experimental data of water content of hydrocarbons and non-hydrocarbons are limited because of difficulties in measurement of the low water content gases at high pressure and low temperature. Mohammadi and coworkers expressed that an accurate estimation of water content can be obtained by gas solubility data, therefore, they overcame the complexities of experimental determination of the water content in natural gases [1]. Due to limited number of measurement data, wide attempts have been made to model and describe the gas-liquid equilibrium in aqueous electrolyte solutions. There are several thermodynamic models which uses the Henry's constant, activity coefficient, and cubic equations of state to obtain more information about the equilibrium conditions. The changes of Henry's constant for the pressure lower than 5 MPa are negligible and it is dominantly affected by temperatures [19]. The high dependency on temperature is obvious at low temperature and also the nonlinear decreasing relationship is observed at high temperatures [25]. Furthermore, there is just a limited number of Henry's constants for hydrocarbon systems at low temperature. According to this fact, there are several drawbacks in applying the Henry's law, whereas it has great ability for accurate prediction of solubility. As an example, it is suitable for dilute solutions or near-ideal solutions [26]. Additionally, this method is correct for single compounds in no chemical reaction conditions for aqueous phase. Another method is cubic EOS which has several advantages such as small number of parameters, computational efficiency, and ease of performance [3,4,21]. The EOSs were proposed originally for pure fluids, after that, their applications were expanded for mixtures by combining the constants from different pure components. This extension can be done by different methods such as Dalton's law of additive partial pressures and Amagat's rule of additive volumes [5]. For complex compounds, there are some limitations in accuracy of EOS which highlight the importance of empirical adjustments by dealing with the binary interaction parameters. In order to determine these parameters, a reliable source of experimental data for vapor-liquid equilibrium is required which induces some uncertainty into EOSs [7].

Due to above discussions, development of an accurate and reliable approach for estimation of solubility of hydrocarbons and non-hydrocarbons in aqueous electrolyte solutions has been highlighted. Nowadays, machine learning approaches have shown extensive applications in different topics [27–35]. This work organizes a novel artificial intelligence method called extreme learning machine (ELM) to estimate solubility of hydrocarbons in aqueous electrolyte mixtures in terms of types of gas, mole fractions of gases, pressure, temperature, and ionic strength.


**Table 1.** Details of experimental hydrocarbons solubility in aqueous electrolyte solutions.

### **2. Methodology**

### *2.1. Experimental Dataset Collection*

In order to construct a highly accurate and comprehensive model capable of estimating the solubility of mixtures of hydrocarbons in aqueous electrolyte solutions, a comprehensive databank was provided based on existing experimental data in Table 1. This databank contains total number of 1175 solubility points for hydrocarbons (881 and 294 points for training and testing phases, respectively) (see Table S1 of data set in Supplementary Materials). According to the literature [1,4,5,9,16–20], the solubility of gases in these systems is highly function of aqueous solutions, pressure, temperature, and gaseous phase composition. The aqueous phase composition was change into ionic strength (*I*) from salt concentrations to reduce dimensions of modeling process. The following equation presents the relationships between ionic strength, valance of charged ions (*zi*), and molar concentration of each ion (*mi*).

$$I = \frac{1}{2} \sum m\_i |z\_i|^2 \tag{1}$$

In this study, the solubility of hydrocarbons is predicted in terms of concentration of components in gaseous mixture, ionic strength of solution, temperature, and pressure.

$$\eta\_{\rm li} = f(\mathbb{C}1, \mathbb{C}2, \mathbb{C}3, \mathbb{C}4, I, P, T, \text{idx}) \tag{2}$$

In which, η*<sup>h</sup>* represents the hydrocarbon solubility in aqueous phase; *C*(1–4) are known as the methane, ethane, propane, and butane mole fraction in gas phase (0–99.99); *I* denotes the ionic strength based on molarity (0–37.35); *T* denotes the temperature in terms of ◦C (1.4–245.15); *P* shows the pressure in MPa (0.3–100), and *idx* symbolizes the index of fraction whose solubility is to be determined (1,2,3,4).

### *2.2. Extreme Learning Machine*

Huang proposed a new intelligence method based on single-layer feedforward neural network (SLFFNN) called extreme learning machine to satisfy the drawbacks of gradient-based algorithms such low training speed and low learning rate. In the ELM algorithm, the hidden nodes are selected randomly and the weights of output of the SLFFNN are calculated by applying Moore–Penrose generalized inverse [36,37].

The scheme of ELM algorithm is demonstrated in Figure 1. By assuming *N* training sets such as (*xi*, *yi*) <sup>∈</sup> <sup>R</sup><sup>n</sup> <sup>×</sup> <sup>R</sup><sup>m</sup> for *<sup>L</sup>* hidden nodes, the SLFFNN algorithms can be written as

$$\sum\_{i=1}^{L} \beta\_i f\_i(\mathbf{x}\_j) = \sum\_{i=1}^{L} \beta\_i f\_i(a\_i.b\_i.\mathbf{x}) \tag{3}$$

In which, *ai* = [*ai*1, ... , *ain*] <sup>T</sup> points to input weights matrix which is related to hidden nodes, β*<sup>i</sup>* = [β*i*1, ... , β*im*] <sup>T</sup> represents the output weights matrix which is related to hidden nodes, and *bi* symbolizes the hidden layer bias.

$$\sum\_{i=1}^{L} \beta\_i f\_i(\mathbf{x}\_j) = H\beta \tag{4}$$

In which, β = [β1, ... , β*L*] and *h*(*x*) = [*h*1(*x*), ... *hL*(*x*)] are known as the hidden layer output matrix and the output weight matrix.

The first step of this model is the random calculation of input weight and the bias of hidden layer for the training phase. Then, for determining these values, the hidden layer matrix is obtained by utilization of input variables. Then, the SLFFNN training is changed to a least-square problem. The ELM algorithms implement regularization theory to define a target function as [38–40]

$$\min L\_{ELM} = \frac{1}{2} \left\| \beta \right\|^2 + \frac{c}{2} \left\| T - H\beta \right\|^2 \tag{5}$$

**Figure 1.** Structure of ELM algorithm.

### **3. Results and Discussion**

In this study, the solubility of hydrocarbons in the aqueous electrolyte phase is determined based on ELM algorithm. To this end, the sigmoid function is set as activation function and the input weights were initialized randomly in range of (−1, 1). Additionally, the number of nodes in the hidden layers was estimated as 30 based on the lowest value of RMSE as determined in Figure 2. As shown, after 30 nodes, by increasing complexity of model, the testing error increased so the optimum structure of the algorithm has 30 nodes to prohibit overfitting.

**Figure 2.** Obtaining optimum structure of proposed algorithm.

In the following, the statistical results of the estimation of hydrocarbon solubility are inserted in Table 2. The following equations are used to achieve this end:

$$\text{Mean relative error} \left( \text{MRE} \right) \\ = \frac{100}{N} \Sigma\_{i=1}^{N} (\frac{\mathbf{X}\_{i}^{\text{actual}} - \mathbf{X}\_{i}^{\text{predicted}}}{\mathbf{X}\_{i}^{\text{actual}}}) \\ \tag{6}$$

$$\text{Root mean square error (RMSE)} = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left( \left( \mathbf{X}\_i^{\text{actual}} - \mathbf{X}\_i^{\text{predicted}} \right)^2 \right)}\tag{7}$$

$$\text{Mean squared error (MSE)} = \frac{1}{N} \sum\_{i=1}^{N} \left( X\_i^{\text{actual}} - X\_i^{\text{predicted}} \right)^2 \tag{8}$$

$$\text{IR}-\text{squared}\ (\text{R2}) \ = 1 - \frac{\sum\_{i=1}^{N} \left( \mathbf{X}\_i^{\text{actual}} - \mathbf{X}\_i^{\text{predicted}} \right)^2}{\sum\_{i=1}^{N} \left( \mathbf{X}\_i^{\text{actual}} - \overline{\mathbf{X}^{\text{actual}}} \right)^2} \tag{9}$$

As shown in Table 2, the MRE, MSE, and RMSE are determined as 22.049, 1.33285 <sup>×</sup> <sup>10</sup><sup>−</sup>8, and 0.0001 for training phase respectively. Moreover, for testing phase, MRE <sup>=</sup> 22.054, MSE <sup>=</sup> 1.05351 <sup>×</sup> <sup>10</sup>−<sup>8</sup> and RMSE = 0.0001 are calculated. The estimated R<sup>2</sup> values are 0.985, 0.987, and 0.985 for training, testing, and overall datasets respectively. These results give the knowledge about the high degree of accuracy for proposed ELM algorithm.


**Table 2.** Statistical analyses of developed model.

On the one hand, the comparison between the estimated and real hydrocarbons solubility in aqueous electrolyte solutions are shown in Figure 3. This depiction demonstrates an excellent agreement between estimated and real solubility values. Figure 4 also represents the regression plot of actual hydrocarbons solubility versus estimated one. A light cloud of data near the 45◦ line expresses the validity and accuracy of ELM algorithm. Additionally, Figure 5 also shows the distribution of relative deviations between forecasted and actual hydrocarbons solubility in aqueous solutions. It can be seen that the ELM outputs deviate slightly from the real solubility and most of relative deviations are near to zero. Furthermore, Figure 6 shows the histograms of relative deviations for training and testing phases. In this demonstration, frequency diagram confirms that most of the error points are close to zero and also cumulative axis express the fact that range of deviation is very limited and the highest slope of the cumulative curve occurred near the zero point.

**Figure 3.** Comparison of actual and estimated solubility of hydrocarbons.

**Figure 4.** Cross plot of actual and estimated solubility of hydrocarbons.

**Figure 5.** Relative deviation between actual and estimated solubility of hydrocarbons.

**Figure 6.** Histogram diagram of relative deviations.

The ELM algorithm implemented in the current work shows an excellent ability in calculation of solubility of hydrocarbons in aqueous phases. One of the important factors which can influence the validation of model is degree of precision of utilized data. In order to clarify the accuracy of solubility databank, the leverage mathematical method is recruited. This method has some rules to identify the suspected solubility data so that a matrix which is known as hat matrix, should be constructed based on formulation [41–45]

$$H = \mathcal{U} \left(\mathcal{U}^T \mathcal{U}\right)^{-1} \mathcal{U}^T \tag{10}$$

In which, *U* symbolizes a matrix of *i* × *j* dimensional. *i* and *j* are known as the number of algorithm parameter and training points which are used for determination of critical leverage limit as

$$H^\* = \mathfrak{Z}(j+1)/i\tag{11}$$

In order to detect the reliable zone, there are two standard residual indexes (−3 and 3) which are used in the leverage method. As shown in Figure 7, the reliable area is bound by these two residual indexes and critical leverage limit. The critical straight lines are shown by red and green colors. This plot is known as William's plot. In this plot, normalized residual is depicted versus hat value which is determined from the main diagonal of aforementioned matrix. It is obvious that the major number of solubility data are located in this area which expresses validation of the hydrocarbon solubility databank.

**Figure 7.** Detection of suspected data for hydrocarbon solubility dataset.

In the most of parametric studies, it is a valuable attempt to identify the effectiveness of all inputs on the target. According to this fact, the sensitivity analysis is employed to investigate effect of concentration of components in gaseous mixture; ionic strength of solution; and temperature and pressure on the solubility of hydrocarbons in aqueous electrolyte systems. To this end, the relevancy factor should be determined as follows for each input parameter [46–54]:

$$r = \frac{\sum\_{i=1}^{n} (X\_{k,i} - \overline{X\_k}) \left(Y\_i - \overline{Y}\right)}{\sqrt{\sum\_{i=1}^{n} \left(X\_{k,i} - \overline{X\_k}\right)^2 \sum\_{i=1}^{n} \left(Y\_i - \overline{Y}\right)^2}} \tag{12}$$

In which *Yi* and *Y* denote the '*i*' th output and output average. *Xk*,*<sup>i</sup>* and *Xk* are known as '*k*' th of input and average of input. Figure 8 shows the relevancy factor for each effective variable of hydrocarbon solubility. It is necessary to explain that the relevancy factor lies in range of −1 to 1 so that the higher absolute value has more impact on hydrocarbon solubility. Furthermore, the positive relevancy factor shows the straight relationship between input and target. The relevancy factors for pressure, temperature, the index of fraction, ionic strength, methane, ethane, propane, and butane mole fraction in gas phase are 0.52, 0.20, −0.48, −0.16, 0.11, 0.06, −0.19, and −0.07 respectively. According to this explanation and results, as pressure, temperature, and mole fraction of methane and ethane increase, the solubility of investigated hydrocarbon increases. Moreover, pressure and mole fraction of ethane in gaseous phase are the most and least effective parameters for determination of solubility of hydrocarbons.

**Figure 8.** Sensitivity analysis for solubility of hydrocarbons.

### **4. Conclusions**

The hydrocarbon solubility in aqueous electrolyte phases at high temperature and pressure conditions is known as a major effective parameter in variety of applications for petroleum industries and chemical engineering. Numerous attempts have been made in the current study to suggest a highly accurate and comprehensive predicting tool on the basis of extreme learning machine to calculate hydrocarbons solubility in wide ranges of operational conditions. Comparing the ELM outputs with a comprehensive real databank which has 1175 solubility points concluded to R-squared values of 0.985 and 0.987 for training and testing phases respectively. The excellent agreements of ELM and real hydrocarbon solubility values express that the ELM algorithm is a valuable tool for design and optimization of various processes that are related to vapor-liquid equilibrium. Furthermore, this study gives more information about the intensity of each input parameter on solubility of hydrocarbons. Due to the aforementioned results, this work has potential application in commercial software packages such as CMG and ECLIPSE for simulation of fluid flow in porous media.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2227-9717/8/1/92/s1, Table S1: data set.

**Author Contributions:** Conceptualization, A.B. and S.S.; Methodology, N.N., A.M., and A.B.; Software, N.N., A.M., A.B., S.S., and I.F.; Validation, N.N., A.M., A.B., S.S., and I.F.; Formal analysis, N.N., A.M., A.B., S.S., and I.F.; Investigation, N.N., A.M., A.B., S.S., and I.F.; Resources, N.N., A.M., A.B., S.S., and I.F.; Data curation, N.N., A.M., A.B., S.S., and I.F.; Writing—Original draft preparation, N.N., A.M., and A.B.; Writing—Review and editing, N.N., A.M., A.B., S.S., and I.F.; Visualization, N.N., A.M., A.B., S.S., and I.F.; Project administration and submission, A.M., Supervision, I.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** We acknowledge the financial support of this work by the Hungarian State and the European Union under the EFOP-3.6.1-16-2016-00010 project.

**Conflicts of Interest:** The authors declare no conflicts of interest.

### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Acid-Base Flow Battery, Based on Reverse Electrodialysis with Bi-Polar Membranes: Stack Experiments**

### **Jiabing Xia, Gerhart Eigenberger, Heinrich Strathmann and Ulrich Nieken \***

Institute of Chemical Process Engineering, University of Stuttgart, Boeblinger Strasse 78, D-70199 Stuttgart, Germany; icvt@icvt.uni-stuttgart.de (J.X.); gerhart.eigenberger@icvt.uni-stuttgart.de (G.E.); heinrich.strathmann@icvt.uni-stuttgart.de (H.S.)

**\*** Correspondence: ulrich.nieken@icvt.uni-stuttgart.de

Received: 19 November 2019; Accepted: 9 January 2020; Published: 11 January 2020

**Abstract:** Neutralization of acid and base to produce electricity in the process of reverse electrodialysis with bipolar membranes (REDBP) presents an interesting but until now fairly overlooked flow battery concept. Previously, we presented single-cell experiments, which explain the principle and discuss the potential of this process. In this contribution, we discuss experiments with REDBP stacks at lab scale, consisting of 5 to 20 repeating cell units. They demonstrate that the single-cell results can be extrapolated to respective stacks, although additional losses have to be considered. As in other flow battery stacks, losses by shunt currents through the parallel electrolyte feed/exit lines increases with the number of connected cell units, whereas the relative importance of electrode losses decreases with increasing cell number. Experimental results are presented with 1 mole L−<sup>1</sup> acid (HCl) and base (NaOH) for open circuit as well as for charge and discharge with up to 18 mA/cm2 current density. Measures to further increase the efficiency of this novel flow battery concept are discussed.

**Keywords:** electrical energy storage; acid-base neutralization flow battery; reverse electrodialysis with bipolar membranes; stack test results

### **1. Introduction**

The efficient storage and recovery of electrical energy is a key issue in the transition of electric energy generation from fossil fuels towards sustainable sources like solar and wind power. Flow batteries present an attractive solution, in particular for stationary and decentralized applications, since the amount of energy stored is separated from energy conversion. Presently, vanadium redox flow batteries are the most important and best-studied type of flow batteries. In redox flow batteries, the electric potential is generated by the electrochemical reactions at the electrodes.

As an alternative, the reverse electrodialysis (RED) has been proposed, where an electrical potential is generated across ion exchange membranes, which separate salt solutions of different ionic concentrations. In most cases studied, seawater and river water have been used as the two solutions [1]. Such a RED-stack consists of a sequence of cells containing one compartment for seawater and one for river water, separated by cation and anion exchange membranes. However, the open voltage of such a cell is only about 0.15 V, which requires a large number of subsequent cells to obtain a reasonable stack voltage.

An interesting alternative RED-concept, which will be considered in the following, is based on the neutralization reaction of acid and base at the internal interface of a bipolar membrane, the reverse electrodialysis with bipolar membranes (REDBP).

In [2] we presented an overview of previous publications on REDBP [3–5], together with our experimental result for a single REDBP cell unit, showing the potential advantages and discussing the present shortcomings of this—so far fairly overlooked—concept.

In the present contribution we extend our experimental results to REDBP stacks, consisting of 5 to 20 cell units. Figure 1 from [2] is a sketch of the stack under charge (lower half) and discharge conditions (upper half). Only one of several repeating cell units between the electrode chambers is shown and only the required main ionic fluxes and the electrode reactions are depicted.

**Figure 1.** Schematic of main ionic fluxes and electrode reactions of the REDBP (reverse electrodialysis with bipolar membranes) stack used, in the lower part during charge and in the upper part during discharge. Only one of several repeating cell units is shown between the electrode compartments. The respective electrode reactions during charge and discharge are indicated at the electrodes. CM: cation exchange membrane, AM: anion exchange membrane, BM: bipolar membrane.

Each repeating cell unit consists of an acid (HCl) and a base chamber (NaOH) at both sides of the bipolar membrane (BM). A salt compartment (NaCl), separated from the acid or base chamber by an anion exchange membrane (AM) or a cation exchange membrane (CM), provides or accepts the respective Cl<sup>−</sup> and Na<sup>+</sup> ions to ensure electroneutrality. To form a REDBP stack, several repeating cell units have to be assembled, with an electrode compartment at each end of the stack.

The main element of each repeating cell unit is the bipolar membrane. Under charge (lower half of Figure 1), a sufficiently high electrical potential has to be applied over the electrodes such that water is split inside the bipolar membrane into H<sup>+</sup> and OH−. H<sup>+</sup> and OH<sup>−</sup> combine with Cl<sup>−</sup> and Na<sup>+</sup> from the adjacent salt chambers to produce HCl and NaOH, which can be stored in separate storage tanks. This corresponds to the well-established process of acid and base generation from salt solutions by electrodialysis with bipolar membranes [6]. The respective electrode reactions are water splitting under generation of OH<sup>−</sup> and of H2 at the negative electrode and with H<sup>+</sup> and O2 generation at the positive electrode. The applied DC power supply transports electrons e− from the positive to the negative electrode.

If the DC power supply is removed ("open circuit conditions") a potential will establish across each of the bipolar membranes in the stack. It prevents the further recombination of H<sup>+</sup> and OH<sup>−</sup> inside each bipolar membrane. The sum of all resulting membrane potentials can be measured between the electrodes of the stack as open circuit voltage (OCV).

If the electrodes are connected via an external load (upper part of Figure 1) the unit can be used as a battery. Now H<sup>+</sup> from HCl and OH<sup>−</sup> from NaOH recombine to water inside each bipolar membrane and generate electric energy; the ionic fluxes in each repeating cell unit reverse. The resulting electric potential across all repeating cell units drives the depicted electrode reactions. The electrode reactions transform the ionic into electronic current, which flows through the load of the external circuit. The reaction details inside the bipolar membrane during charge and discharge have been discussed in [2]. They depend on the permselectivities and on the fixed charge densities of the ion exchange membranes used.

Commercially available membranes currently limit acid and base concentrations to 2 M, but long-term stability allows only for operation up to 1 M [2]. Using 1 M acid and base, the theoretical power density is 11 kWh/m3. Energy densities of vanadium redox flow batteries are in the order of 25 kWh/m3 [7].

### **2. Experimental**

In our stack experiments hydrochloric acid (HCl) and sodium hydroxide solutions (NaOH)of different concentrations have been used as the respective acid and base, with 0.5 mole/L NaCl as the resulting salt solution. As shown in Figure 2, the respective solutions were circulated through the corresponding compartments from external storage tanks. Through the electrode compartments, a 0.25 mole/L sodium sulphate solution (Na2SO4was circulated instead of NaCl, to avoid chlorine formation in the electrode reactions. In the set-up of Figure 1, tested in this contribution, an additional salt compartment (NaCl) was placed between the sequence of repeating cell units and each electrode compartment, separated from the electrode compartments by a cation exchange membrane (CM). This should prevent the penetration of Cl− into the electrode compartments, which would lead to a release of Cl2 in the electrode reactions. Prime advantages of the chosen electrolytes are their low price, ample availability and low ecological concern.

**Figure 2.** Parallel flow concept for acid, base and salt solution through four repeating cell units in the middle of a stack. The black arrows indicate the flow directions of acid (HCl), of base (NaOH) and of salt solution (NaCl), and the colored arrows indicate the main fluxes of specific ions during discharge, also at OCV.

Table 1 shows the expected electrode reactions during charge and discharge. Their standard potentials represent a voltage sink of 1.23 + 0.83 = 2.06 V, both at charge and at discharge. This loss could be reduced by modified electrode reactions with gas-consuming electrodes, but this was not in the scope of the present study.


**Table 1.** Electrode reactions during charge and discharge.

Since the open circuit voltage (OCV) obtained in the single cell experiments with 1 M acid and base was 0.76 V, only a stack with three and more cell units would be able to compensate the electrode losses and deliver a net voltage during discharge. The influence of the electrode losses obviously decreases with an increasing number of cell units in the stack. However, the stack voltage is not fully proportional to the number of cell units. This is partly due to the so-called shunt currents, which flow through the feed and exit lines of the different electrolytes, transforming electrical energy into heat.

### *2.1. Parallel Flow Concept*

In our set-up, a parallel flow concept for the electrolytes has been used, as shown schematically in Figure 2 for four repeating cell units in the middle of a stack. The parallel flow concept has the advantage of reduced pressure drop, compared to a serial flow concept where all of the acid, base and salt solution flows consecutively through adjacent cell units. The potential disadvantage is that the flow through adjacent parallel chambers may differ, if the flow resistances of the chambers are not equal. In addition, a parallel flow concept leads to larger ionic shunt currents through the feed/exit lines, as discussed in the next subsection. The colored arrows in Figure 2 indicate the main ionic fluxes across the membranes at discharge. At a smaller extent, similar fluxes also occur at OCV conditions. They are caused by a breakthrough of co-ions (Cl<sup>−</sup> through CM and H<sup>+</sup> through AM), as discussed more detailed in [8]. This breakthrough is responsible for the fact that the open current voltage (OCV) across a repeating cell unit for 1 mol/L HCl and NaOH at 25 ◦C decreases from the theoretical value of 0.828 V to a measured value of 0.764 V.

Like conventional electrodialysis stacks, the stack is formed from cell frames with pathways for the respective electrolytes, grid spacers in the cell frame windows, with the respective membranes and with rubber gaskets between cell frames and membranes. The cell frames for the stack experiments of 0.5 mm thickness have been obtained from DEUKUM Co, Frickenhausen, Germany (www.deukum.de), see Figure 3. They proved to be well suited for our lab-scale experiments, since they allow feeding of four different electrolytes through the holes in the frame. Their active membrane area is about 100 cm2. As for the single cell tests, the membranes (fumasep® FKB, fumasep® FAB and fumasep® FBM,) have been provided by Fumatech Co., Bietigheim-Bissingen, Germany. Their properties are listed in Table 2.

**Table 2.** Properties of membranes fumasep® FAB, fumasep® FBM with fumasep® PEEK as reinforcement, provided by Fumatech [9].


**Figure 3.** Subsequent cell frames for one repeating cell unit: (**a**) for up-flow of HCl, (**b**) for up-flow of NaOH and (**c**) for cross-flow of the NaCl solution.

Figure 3 shows the flow through subsequent cell frames, (a) for up-flow of HCl, (b) for up-flow of NaOH and (c) for cross-flow of the NaCl solution. Together with the respective membranes, separated by rubber gaskets and grid spacers, the three frames form the repeating cell units shown schematically in Figure 1. The feed and exit flows of the respective electrolytes pass through the holes in the cell frames as indicated by the respective colors in Figure 3. All electrolytes are circulated by separate pumps from reservoirs through the respective compartments as indicated in Figure 2. Further details of the design are given in [8].

### *2.2. Voltage Measurements*

The most direct information during stack operation can be obtained from voltage measurements along the stack. In the single repeating cell unit studied in [2], Haber-Luggin capillaries with calomel electrodes have been used to measure the voltage drop across the respective membranes. To reduce the ample space required for the Haber-Luggin capillaries, in the present stack experiments the voltages between different repeating cell units have been measured by thin Pt wires. The insulated Pt wires were introduced and sealed between two rubber gaskets. To minimize concentration effects, the Pt wires were only positioned in the acid compartments, where ionic conductivities are high and about equal throughout the stack.

In the following, voltage measurements will be presented where Pt wire probes were positioned in the center of the 1st, the 2nd, the 6th, the 11th and the 19th acid cell of the stack. This allows measurement of the voltages across the 1st repeating cell unit, across the first five-cell unit (next to the electrode) and across the second five-cell unit (in the middle of the stack). In order to determine the total voltage of the 20-cell stack, the measured first cell voltage was added to the voltage measured between the first and the 19th acid compartment, assuming that the voltage across the repeating cell units next to both electrodes was equal.

### *2.3. Operating Conditions*

In all cases, Na2SO4 electrolyte of 0.25 M was circulated through the electrode chambers and 0.5 M NaCl electrolyte was circulated through the salt chambers, while the concentration of HCl and NaOH, circulated through their chambers, was varied between 0.5 M and 1 M. If not specified differently, the empty-space flow velocity through the electrolyte chambers was adjusted to a mean value of about 2 to 3 cm/s.

### **3. Results and Discussion**

In this section, the measured voltages in stacks consisting of 5 to 20 repeating cell units are presented and discussed for different acid and base concentrations, for different charge and discharge currents and for different charge/discharge periods.

### *3.1. Open Circuit Voltage (OCV)*

With the Pt wire probes, the open circuit voltage can be measured between different cell units as mentioned above, but OCV for the whole stack can also be measured directly between the electrodes of the stack. Table 3 presents OCV results for increasing stack size between 1 and 20 repeating cell units if the electrolytes flow through the stack continuously (constant ion concentrations). In an ideal stack, OCV should amount to the theoretical single cell voltage (0.828 V for 1 M acid and base at 25 ◦C) times the number of repeating cell units. With the measured single cell voltage of 0.764 V, the theoretical values of the stack are given in the second row. The actually obtained OCV values, measured across the electrodes, are given in the third row. While the measured values are still close to the calculated values for up to 10 repeating cell units, an increasing voltage loss is observed if the cell unit number is raised to 15 and further to 20 cells. Reasons for this trend are discussed in the next section.

**Table 3.** Calculated and measured open circuit voltage (OCV) for 1 M acid and base at 25 ◦C, depending on the number of repeating cell units.


### *3.2. Self-Discharge at OCV*

The experimental results of Table 2 have been obtained with stacks where the required electrolytes were continuously pumped through the stack, compensating the effect of any side reactions or shortcuts, which may consume acid and base at OCV. In the OCV experiments reported in Figure 4, the electrolyte flows have been stopped at the beginning of the experiment. Now the resulting drop of the voltage in stacks with different numbers of repeating cell units can be analyzed over time. This should provide indications as to why OCV across stacks with a larger number of repeating cell units increases less than proportional to the number of units, as shown in Table 2. A self-discharge due to limited membrane permselectivities, as indicated by the ionic flux arrows at OCV in Figure 2 and discussed in detail in [2], could contribute to the effect. However, it should not be responsible for the less than proportional increase of OCV with the number of repeating cell units.

Figure 4a–d show the results of the self-discharging tests of (a) 5-cell, (b) 10-cell, (c) 15-cell and (d) 20-cell stacks with initial 1 M acid and base. In a 5-cell stack (Figure 4a), the self-discharge is very low, and OCV remains close to 4V during the whole measurement period. In a 10-cell stack (Figure 4b), the stack voltage is about doubled and a slow self-discharge is followed by an accelerated discharge after about 30 min down to a level of 4 V. If the cell number is further increased (Figure 4c,d), the accelerated drop from the higher stack voltage down to about 4 V occurs at about 20 min for the 15-cell stack, and at about 12 min for the 20-cell stack. This faster decrease of OCV is consistent with the assumption that ionic shunt current through the feed/exit lines of the electrolytes should be responsible for the OCV decrease. A similar influence of shunt currents through parallel feed/exit lines of the electrolytes is well known from redox flow batteries [10–12].

**Figure 4.** Self-discharge test at 25 ◦C for different stack sizes and different acid and base concentrations. (**a**) 5-cell stack, 1M acid and base, (**b**) 10-cell stack, 1M acid and base, (**c**) 15-cell stack, 1M acid and base, (**d**) 20-cell stack, 1M acid and base, (**e**) 20-cell stack with 0.5 M acid/base concentration. OCV\_Stack (black) has been measured between the electrodes, OCV\_AllSCs (blue) between the first and the last repeating cell unit of the stack, OCV\_5cells\_1st (green solid) and OCV\_5cells\_2nd (green, dotted) have been measured across the first and second 5-cell unit.

In stacks with 10 and more cells, OCV shows an accelerated decline to about 4 V after a certain period of low decline, which becomes shorter with increasing cell number. An explanation can be drawn from Figure 4d. Here, the difference between OCV over the first 5 cells of the stack (OCV\_5cells\_1st, green solid line) and OCV over the second 5-cell package (OCV\_5cells\_2nd, green dotted line) is shown. The voltage over the 5 middle cells was smaller than that over the first 5 cells from the very beginning. It dropped to (almost) zero at the time of the accelerated decline of the total cell voltage, whereas the voltage over the first 5 cells was less affected. The voltage decline indicates that the acid and base have been completely consumed (neutralized) in the middle of the stack, but are still present in the cells next to the electrodes. This can be explained by the ionic shunt currents through the feed/exit lines of the electrolytes. Since H<sup>+</sup> and OH<sup>−</sup> have much higher ionic mobility than the other ions present, only the influence of H<sup>+</sup> and OH<sup>−</sup> shunt currents in the acid and base feed/exit lines will be considered in the following, using the simplified picture in Figure 5.

**Figure 5.** (**a**) Simplified sketch of ionic shunt currents of H<sup>+</sup> (red arrows) and OH<sup>−</sup> (blue arrows) through the feed/exit lines of acid and base, and compensating ionic currents inside of four repeating cell units in the middle of the stack. The different arrow thickness indicates the increased shunt currents in the stack center. (**b**) Corresponding electric potential profiles ϕ across the repeating cell units (black) and in the acid (red) and the base feed/exit lines (blue). The arrows indicate the direction of the ion fluxes between cell compartments and feed/exit lines.

In Figure 5a only four repeating cell units in the middle of the stack are shown, together with their feed/exit lines for acid (red) and base (blue). At OCV, a voltage profile ϕ (black) develops across the stack as shown in Figure 5b. It mainly consists of the voltage steps in each of the bipolar membranes. The voltage gradient across the stack causes a counter-movement of H<sup>+</sup> ions and of OH<sup>−</sup> ions along their feed/exit lines as indicated by the red and blue arrows in Figure 5a. The ions originate from the acid/base compartments at one side of the stack and turn back into the respective compartments at the opposite side. This leads to an ion flux maximum in the stack center. To ensure electroneutrality, the ion fluxes through the feed/exit lines have to be compensated by fluxes of opposite direction through the cell compartments. This leads to the indicated fluxes with neutralization of H<sup>+</sup> and OH<sup>−</sup> inside the bipolar membranes, which is maximal in the middle of the stack. Comparable to a simple mixing of acid and base, this neutralization energy is converted into heat, representing a self-discharge with loss of electrical energy. The fluxes of H<sup>+</sup> and OH<sup>−</sup> into the bipolar membranes have to be compensated by equivalent fluxes of Cl<sup>−</sup> and Na<sup>+</sup> into the adjacent salt compartments.

Arrows in Figure 5b qualitatively depict the driving potentials between the feed/exit lines and the potentials on both sides of the respective bipolar membranes. This explains why the ion fluxes of H<sup>+</sup> and OH- between the stack cells and the feed/exit lines change direction in the middle of the stack, as indicated by the arrows in Figure 5a,b. The fluxes of H<sup>+</sup> and OH<sup>−</sup> through the feed/exit lines accumulate in the stack center as indicated by the arrow thickness. The compensating fluxes inside the repeating cell units are therefore also strongest in the middle of the stack. This leads to the pronounced self-discharge and the resulting drop of OCV across the bipolar membranes in the middle of the stack.

In summary, the ionic shunt currents of H<sup>+</sup> and OH<sup>−</sup> ions in the acid/base feed exit lines are driven by the potential along the acid/base feed/exit lines and result in a neutralization reaction in the bipolar membranes. Estimation shows that in the stack used, shunt currents due to H<sup>+</sup> migration through the feed/exit lines of HCl amount to more than 60% and due to OH− migration through the feed/exit lines, NaOH amounts to more than 30% of the total shunt current [8]. This leads to the consumption of acid and base, primarily in the stack center, which reduces the respective cell voltages. If under OCV the flow of acid and base is stopped, it leads to the gradual neutralization of all acid and base in the stack center as can be concluded from the voltage drop to zero of "OCV\_5cells\_2nd" after about 12 min in Figure 4d.

In Figure 4e, results are presented where the initial acid/base concentration has been reduced from 1 to 0.5 M. This reduction of 50% should result in an accelerated discharge in about half the time. However, since the reduced concentrations also lead to reduced ionic conductivities and hence to a reduced shunt current, only a small difference between Figure 4d,e can be observed. Again, OCV across the second 5-cell unit in the stack center drops to zero, now after around 10 min.

The above details of the observed OCV decline after stopping the acid/base feed help to elucidate the influence of shunt currents. However, it should be mentioned that shunt currents are only of importance for the REDBP battery during operation. The consumption of acid and base by self-discharge at OCV is limited to the amount present in the stack and does not affect the concentration of acid and base in the storage tanks. This is a general advantage of flow batteries.

To determine all the details of ionic flux interactions, a numerical simulation of the whole stack would be required, where ionic species balances for all ionic species as well as electroneutrality have to be enforced in the whole stack. But this is out of scope for the present contribution.

### *3.3. Charge and Discharge Behavior*

Figure 6 shows the charge and discharge behavior of a 20-cell stack with 1 M acid/base concentration at 25 ◦C. Starting with OCV conditions, after 1 min a charge current density of 9 mA/cm2 is imposed for 5 min, followed by zero current for 1 min and a discharge current density of 9 mA/cm2 for another 5 min. The voltage measured across the stack (black line) is higher under charging and lower under discharging than the voltage measured with the Pt wire probes between all repeating cell units, because the electrode losses need to be added to the blue curve under charge and subtracted under discharge. The mean voltages measured along the 5 min charge and discharge periods are constant. Under constant acid and base concentrations, the measured voltages even remain constant during extended charge and discharge periods, indicating a stable operation. Again, the voltage measured across the second 5-cell package (V5cells\_2nd, green dotted line) is lower than V5cells\_1st (green solid line) at OCV as well as under charge/discharge, indicating the increased influence of the shunt currents in the stack center.

The absolute difference between the blue and the black curves should, however, be about equal for charge and discharge and vanish at OCV. This would be the case if the blue lines were shifted upwards by 0.7 V, leading to an absolute difference between both curves at charge and discharge of about 3.4 V. This seems to be a reasonable estimate of the losses at the electrodes and across the electrode compartments during charge/discharge with 9 mA/cm2, since under ideal conditions electrode losses of 2.07 V have been predicted by Table 1. The required voltage shift of 0.7 V points to a measurement error of the Pt wire probes (blue curves). It is well known that measurements with Pt wire probes are less accurate than with Haber-Luggin electrodes as applied in [2]. In the summary of the measured results, presented in Figure 7, this voltage shift has been considered and compensated.

**Figure 6.** Charge and discharge behavior of a 20-cell stack with 1 M acid/base at 25 ◦C and 9 mA/cm2 current density. The black line is the stack voltage measured across the electrodes; the blue line is the stack voltage across the 20-cell units. The voltages measured across the first repeating cell unit next to the electrode (red), the first 5-cell unit (green, solid) and the second 5-cell unit (green, dotted) are also shown.

**Figure 7.** Discharge voltages and net discharge power densities of the 20-cell stack with 1 M acid/base at 25 ◦C, plotted over discharge current density.

In Figure 8, measured voltages are presented for increasing current densities under charge (increasing voltages) and under discharge (decreasing voltages). At each current density, the measurement continued until a stable voltage was established. Figure 8a shows the results of a 20-cell stack and Figure 8b of a 15-cell stack, for 1M acid and base. Both results show the same general trends. The voltages measured across the electrodes of the whole stack are given by black dots and the measurements across all repeating cell units (excluding the electrodes) by blue dots. In addition, measurements across the first 5-cell unit (next to one electrode) are given by green solid

lines and the measurements across the second 5-cell unit (in the middle of the stack) by green broken lines. Measurements for one repeating cell unit next to one electrode are displayed in red. If a shift of +0.7 V is imposed on the blue curves, as discussed before, the electrode losses between charge and discharge are about equal. They increase moderately with increasing current density.

**Figure 8.** Current density measurements during charge and discharge with 1 M acid and base at 25 ◦C and 0.1 MPa at steady state. Current-voltage curves of positive slope represent charging, of negative slope discharging. Results for a 20-cell stack are shown in (**a**) and for a 15-cell stack in (**b**). Displayed are the measured voltages, VSC (across BP (bipolar membranes) of the first single cell, red), V5cell\_2nd (across the second 5-cell unit, green dotted), V5cell\_1st (across the first 5-cell unit, green solid), VAllSCs (across all cell units measured between Pt electrodes, blue) and VStack (stack voltage between electrodes, black).

A comparison of Figure 8a,b shows that the voltage losses increase with the number of unit cells. This is a consequence of the fact that the shunt current accumulates in the stack center, as discussed qualitatively with Figure 5. This shunt current has to be compensated by an increased neutralization reaction of H<sup>+</sup> and OH<sup>−</sup> in the bipolar membrane. The unit cells in the stack center therefore contribute less to the stack voltage if the cell number increases. A decrease of the voltage of the second 5-cell unit, located in the middle of the stack, compared to the first 5-cell voltage is also obvious in Figure 8a, demonstrating the influence of the increasing shunt current in the middle of the stack, both at charge and discharge.

### *3.4. Discharge Power Density and E*ffi*ciency*

The focus of the present contribution is on the discharge behavior of flow batteries. Here, as in [2], power density PD, defined as net electric power, divided by stack volume, is useful for characterizing the discharge. For the stack volume, only the active area of the membranes and the thickness of each repeating cell unit (1.96 mm) have been considered, while the volume of both electrode chambers has not been included. In Figure 7, discharge voltage and net discharge power density of the tested 20-cell stack with 1 M acid/base at 25 ◦C are displayed over discharge current density. As explained in Section 3.3, the voltages measured across all single cells (the blue dotted line) have been moved up by 0.7 V (compared to Figure 8), to compensate for measurement errors of the Pt probes.

Considering an ideal behavior without losses across the stack (zero resistances) and with the theoretical single cell voltage of 0.828 V, a constant voltage of 16.56 V would result, leading to a linear increase of power density over current density, marked by the solid red line. The actually measured voltages across all 20 cells of the stack over current density lead to the power densities as displayed by the red line with measurement bullets. The losses partly result from the Ohmic resistances across the membranes and the electrolyte compartments, as already considered in the single-cell experiments in [2]. Another part results from the shunt currents, discussed in Section 3.2.

If we include the measured voltage drop over the electrodes, the pink line with bullets results in the measured power density of the stack. The difference between the red and the pink line is attributed to additional losses in the electrode compartments. The fact that the total losses are substantially higher than predicted from the single-cell experiments can be attributed to the shunt currents through the feed/exit lines. In addition, a nonuniform flow distribution of the electrolytes through adjacent cells could have contributed. This can be concluded from differences between measured voltages of adjacent cell units.

The highest power density for the 20-cell stack of 15 mW/cm3 has been reached at 10 mA/cm2, corresponding to a stack power output of 6 W with 36% voltage efficiency. Neglecting the electrode chamber losses, the highest stack power goes up to 9 W or 24 mW/cm−<sup>3</sup> with 50% voltage efficiency.

The power density decreases at higher current densities due to increasing resistive losses. As observed in the single-cell experiments [2], discharge current densities well above 20 mA/cm2 resulted in unstable behavior, with voltage breakdown due to water accumulation inside BP. However, this is above the discharge power densities reached in the stack experiments.

### **4. Conclusions**

The experimental results obtained in this study show that the single cell results presented in [2] can be extrapolated to multicell stacks, but additional losses have to be considered. Most importantly, the influence of shunt currents has to be taken into account. As demonstrated with open current experiments under self-discharging conditions in Section 3.2, the influence of shunt currents through the parallel acid/base feed/exit lines increases with increasing number of cell units. It is most pronounced in the middle of the stack, where it can substantially reduce the open cell voltage.

The single cell experiments reported in [2] showed that at higher acid/base concentrations, a breakdown of the discharge voltage occurred, if the current density was raised above a certain value between 20 and 40 mA/cm2. This was attributed to a delamination of the bipolar membrane, caused by the excessive formation of water through the neutralization reaction of H<sup>+</sup> and OH<sup>−</sup> in the reaction layer. However, these current densities were well above the discharge current densities reached in the stack experiments of Figures 6 and 8. Delamination of bipolar membranes during discharge was not observed in the stack results presented.

Figure 7 summarizes the experimental discharge results for the 20-cell stacks and points to possibilities for improvement. A reduction of the influence of shunt currents would directly reduce the difference between the (ideal) red solid line and the measured red dotted line. A straightforward approach would be to change the acid and base flow pattern from parallel to a serial flow, where all of the acid and base flows consecutively upwards in one and downwards in the next acid/base chamber. Then, minor shunt currents could only occur between adjacent cells, driven by the single cell voltages but not by the total stack voltage, as in case of the parallel feed/exit lines. Preliminary experiments at open current under self-discharge (similar to Figure 4) proved that the self-discharge of a 5-cell stack with serial flow of acid and base was considerably slower, compared with a similar stack with parallel flow. However, for the present design the pressure drop across the acid/base feed/exit lines was substantially increased. This would result in an unacceptable pressure drop for a 20-cell stack with serial flow.

Since Figure 4 has shown that a 5-cell stack with parallel flow was only marginally affected by shunt currents, a combination of parallel and serial flow could be a compromise between pressure drop and shunt current. The parallel flow through the acid and base compartments should change direction after, for example, every 5 to 10 cells. Electrode losses (the difference between the red and the pink lines in Figure 7) could be reduced, e.g., by improved, gas-consuming electrodes [13].

With respect to higher power densities, the positive influence of larger concentrations of the electrolytes has to be balanced against reduced permselectivities of the ion exchange membranes used. Here, the challenge is to develop monopolar and bipolar membranes with increased fixed-charge concentrations as already discussed in [2]. In the single-cell experiments [2], water accumulation at the BP interface resulted in a breakdown of the cell voltage at current densities above about 30 mA/cm<sup>2</sup> for 1M acid/base. This was well above the current densities reached in the stack experiments. Nevertheless, an improved water transport through BP will be a further challenge for BP membranes optimized for REDBP.

**Author Contributions:** Conceptualization, G.E., H.S. and U.N.; formal analysis, J.X.; investigation, J.X.; methodology, J.X. and H.S.; resources, J.X. and H.S.; supervision, G.E., H.S. and U.N. All authors have read and agreed to the published version of the manuscript.

**Funding:** Jiabing Xia gratefully acknowledges the Ph.D. scholarship from the GREES program at University of Stuttgart in Germany.

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

### **List of Abbreviations**


### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
