CFD Modeling and Experimental Validation of an Alkaline Water Electrolysis Cell for Hydrogen Production

: Although alkaline water electrolysis (AWE) is the most widespread technology for hydrogen production by electrolysis, its electrochemical and ﬂuid dynamic optimization has rarely been addressed simultaneously using Computational Fluid Dynamics (CFD) simulation. In this regard, a two-dimensional (2D) CFD model of an AWE cell has been developed using COMSOL ® software and then experimentally validated. The model involves transport equations for both liquid and gas phases as well as equations for the electric current conservation. This multiphysics approach allows the model to simultaneously analyze the ﬂuid dynamic and electrochemical phenomena involved in an electrolysis cell. The electrical response was evaluated in terms of polarization curve (voltage vs. current density) at di ﬀ erent operating conditions: temperature, electrolyte conductivity, and electrode-diaphragm distance. For all cases, the model ﬁts very well with the experimental data with an error of less than 1% for the polarization curves. Moreover, the model successfully simulates the changes on gas proﬁles along the cell, according to current density, electrolyte ﬂow rate, and electrode-diaphragm distance. The combination of electrochemical and ﬂuid dynamics studies provides comprehensive information and makes the model a promising tool for electrolysis cell design.


Introduction
The current fossil fuel-based energy system has to face the most important challenge of our time: climate change. The energy transition to a zero-carbon model requires an ambitious growth plan for renewable energies [1]. However, renewable energy sources (mainly wind and solar power) involve a critical issue that needs to be addressed [2]; their dependence on weather conditions makes them very intermittent, causing a mismatch between supply and demand. Furthermore, in the medium-long term, it is expected that these discrepancies can produce failures in the electricity grid. So, to prevent these problems, adequate energy storage systems are required.
Regarding long-term storage of renewable energies, only chemical energy carriers can be considered nowadays as real solution. Among them, hydrogen has been identified as the best candidate because it can be more efficiently produced from renewable energy surplus in combination with water electrolysis [3]. Furthermore, hydrogen can be transported by its own hydrogen distribution grid as admixture in the natural gas grid, or it can be also stored in appropriate facilities. When required, hydrogen can be efficiently transformed into electricity by fuel cells and returned to the electricity grid or used as fuel for the transport sector or even for new industrial processes [4,5]. In this way, electrolysis enables large scale production of green hydrogen and electrical energy storage.
With respect to water electrolysis, different technologies are commercially available: alkaline, proton exchange membrane (PEM) and solid oxide (SO) electrolysis cells, according to the electrolyte used. Among them, alkaline water electrolysis (AWE) stands out as being a highly developed technology in the industry and the main way to obtain "green hydrogen". Compared to other technologies, the key advantages of alkaline electrolysis are its proven durability, maturity, and low specific costs. Nowadays, electrolyzers are commercially available in the MW range, offering high hydrogen production capacities. The cost of alkaline electrolysis systems is currently in the range of 600-1000 €/kW at the MW level. Regarding the disadvantages, these are mainly the low current densities and the problems derived from working in dynamic conditions. In practice, hydrogen production is limited to an operating range of 20% to 100% of the nominal power in order to prevent the formation of flammable mixtures due to gas diffusion (crossover) through the diaphragm at very low current densities (<0.1 A/cm 2 ) [2, [6][7][8].
For these reasons, the research and development efforts in this technology are mainly focused on increasing the current density and solving the problems derived from working in dynamic conditions. The fluctuations in power supplies to the electrolyzer could cause problems such as the generation of explosive mixtures, corrosion of materials, lower efficiency, pressure drops, and temperature changes, etc. Therefore, the design of alkaline electrolyzers powered by renewable energy is a critical issue [9]. In this sense, the optimization of the cell design and the space between electrodes is a fundamental step to reduce ohmic overpotentials and increase the nominal current density of the electrolyzer [6,10,11].
Therefore, in the last few decades, a significant effort has been made to model, characterize, and analyze the operation of alkaline electrolyzers [2]. In this sense, Olivier et al. [12] conducted an exhaustive and comprehensive analysis of the existing modelling works concerning low temperature electrolysis systems. Therefore, numerous empirical and semi-empirical models have been proposed. One of the most widely used models was proposed by Ulleberg [13], which provides a mathematical description of the polarization curve with parameters fitted empirically. However, these models usually do not include geometric aspects, nor are they capable of modelling those effects that depend on the complete design of the cell.
In this context, Computational Fluid Dynamics (CFD) simulation is a critical and powerful design tool that can be used to improve flow distribution and minimize energy consumption, allowing efficient electrolyte inlet and gas removal from the cell as soon as possible [9]. CFD models allow for the description of thermal, electrochemical, and fluid dynamic phenomena that occur at the same time, from a multiphysics perspective, according to the considered cell geometry. In this way, the voltage required to carry out electrolysis, the generation of H 2 /O 2 bubbles at the electrodes, the flow distribution within each electrolysis cell, and the heat generation according to current density are some of the processes that happen simultaneously in an electrolyzer during its operation.
However, few analytical models based on CFD simulation that describe the behavior of alkaline electrolysis cells have been reported and validated. With this in mind, Aldas [14] applied a two-phase mathematical model to the numerical investigation of gas evolution in a vertical electrochemical cell using the PHOENICS computer code. In other work, Mat et al. [15] investigated the hydrogen evolution, flow field, and current density distribution in an electrochemical cell using a two-phase flow model. The predicted results satisfactorily agreed with data available in the literature. Later, Hawkes at al. [16] used a three-dimensional (3D) CFD electrochemical model to study a high-temperature electrolysis stack using the commercial code FLUENT. Alternatively, Hreiz et al. [17] simulated bubble hydrodynamics on vertical plane electrodes using an Euler-Lagrange CFD approach, and the results were compared with gas velocity fields obtained by a Particles Image Velocimetry (PIV) algorithm, which were in good agreement. Later, El-Askary et al. [18] modeled the hydrodynamics characteristics of the hydrogen evolution process to predict the flow characteristics, gas release rate, and void fraction distribution in electrolysis cells; subsequently, an experimental study was also carried out to verify the calculations made. Zarghami et al. [19] used an Euler-Euler model in FLUENT to simulate the multiphase flow in an alkaline water electrolyzer and compared the results to existing experimental data. Recently, Le Bideau et al. [20] developed a two-phase hydrodynamics model and validated the results with the experimental velocity profiles measured using the Laser Doppler Velocimetry (LDV) method.
Considering the increasing industrial development of water electrolysis that is expected in the coming years [21], it is crucial to consider the effect of renewable energies on the performance of electrolyzers. For these reasons, in the present paper, a two-dimensional (2D) CFD model of an alkaline electrolysis cell implemented in COMSOL is reported based on a work previously published by the authors [9]. The aim of the work is to show the capacity of the model to predict the behavior of the electrolysis cell under different operational conditions and also show how it can influence its performance. Thus, a novel CFD model that simultaneously simulates the electrochemical and fluid dynamic behavior is described, taking a real alkaline water electrolysis cell as a reference. In this way, the response of an AWE cell is simulated under different temperatures and electrolyte concentrations. Polarization curve and gas-liquid fraction can be analyzed at the same time using the presented model. Furthermore, the model shows a wide versatility regarding geometrical parameters, being able to consider the effect of modifying the electrode-diaphragm distance. In order to evaluate the accuracy of the model, simulated results were validated against the real response of the AWE cell at different operating conditions, reporting a very good agreement with a mean relative absolute difference error lower than 1%. In addition, fluid dynamics phenomena such as gas-liquid distribution, turbulence, or formed gas profile were predicted. Therefore, the capacity to simulate both phenomena (electrochemistry and fluid dynamics) in the same model make it a complete design tool.

Alkaline Water Electrolysis
Water electrolysis deals with decomposition of water into hydrogen and oxygen by passing an electric current (DC) between two electrodes separated by an aqueous electrolyte with high ionic conductivity [22]. The overall reaction for water splitting is: (1) Figure 1 shows a scheme of an alkaline electrolysis cell. The two electrodes (anode and cathode) are immersed in a KOH aqueous solution (electrolyte), and they are separated by a porous diaphragm that allows the ionic transport (OH − ) but is impermeable to gases. In a typical operation, electrolyte enters into the anodic and cathodic compartments by the bottom inlets. At the anode, oxygen bubbles are produced and, at the cathode, hydrogen is generated. These bubbles grow until they detach from the electrode surface. The mixture of generated gases and electrolyte then leaves each compartment through the upper outlets [9]. In this way, in an alkaline electrolysis cell, when the required potential between the electrodes is applied, the following semi-reactions take place simultaneously in the cathode (hydrogen evolution reaction (HER)) and anode (oxygen evolution reaction (OER)) [  In this way, in an alkaline electrolysis cell, when the required potential between the electrodes is applied, the following semi-reactions take place simultaneously in the cathode (hydrogen evolution reaction (HER)) and anode (oxygen evolution reaction (OER)) [6]: Thus, for alkaline water electrolysis to occur, a minimum voltage is required (reversible voltage, U rev ), which, according to the semi-reactions potential (E 0 ) from Equations (2) and (3), is equal to 1.23 V (E 0,OER -E 0,HER ) at standard conditions (1 bar and 25 • C). This value, in agreement with the second law of thermodynamics, corresponds to ∆g • = 273.2 kJ/mol.
However, the real cell voltage (U) is always higher than the latter because of irreversibilities or overpotentials. Therefore, the real cell voltage can be defined as the sum of reversible voltage (U rev ) and the overpotentials (η), as shown in Equation (4): The term η is the sum of activation, ohmic, and concentration overpotentials. These overpotentials are defined as follows [23]: (1) Activation overpotentials: related to activation energies of hydrogen and oxygen formation reactions on the surface of electrodes; (2) Ohmic overpotentials: sum of the electrical resistance of several components such as electrodes, current collectors, etc., and the transport resistance related to gas bubbles, ionic transfer in the electrolyte, and resistivity of the diaphragm; (3) Concentration overpotentials: due to mass-transport limitations occurring on the surface of the electrodes at high currents.
The total contribution of these overpotentials to the cell voltage (U) can be analyzed through the polarization curve of an electrolysis cell, as represented in Figure 2.

Model Geometry and Mesh
As mentioned previously in Figure 1, a conventional alkaline water electrolysis cell is divided into two compartments separated by a diaphragm. In the cathodic chamber, the reduction of water to produce H2 takes places, and the oxygen evolution reaction occurs at the anode. The function of the diaphragm is to avoid the mixture of the two gases while maintaining a low resistivity.

Model Geometry and Mesh
As mentioned previously in Figure 1, a conventional alkaline water electrolysis cell is divided into two compartments separated by a diaphragm. In the cathodic chamber, the reduction of water to produce H 2 takes places, and the oxygen evolution reaction occurs at the anode. The function of the diaphragm is to avoid the mixture of the two gases while maintaining a low resistivity.
Taking as reference a rectangular laboratory electrolysis cell (Figure 3a), the geometry of the model was built following fluid dynamic requirements. In this way, simplifications were made in order to reduce the model complexity ( Figure 3b). As result, a good approximation can be made just by 2D geometry, which allows an optimal study of the main involved variables (Figure 3c). The electrode surfaces (33 mm height) were considered by two simple boundary conditions (1 and 5 in Figure 3c) on both sides of the cell. These domains work as current collectors where potential was applied. Regarding the separator (3 in Figure 3c), it was introduced in the geometry as a thin rectangle (0.5 mm) that separates both compartments. In this case, only ionic conductivity is allowed in this component. Concerning the anodic and cathodic compartments (2 and 4 in Figure 3c), they were defined as two thick rectangles (1.5, 4, or 10 mm depending on the model considered) between the separator and the corresponding electrode with the electrolyte inlet at the bottom and the biphasic mixture (electrolyte and gases) outlet at the top.

Model Geometry and Mesh
As mentioned previously in Figure 1, a conventional alkaline water electrolysis cell is divided into two compartments separated by a diaphragm. In the cathodic chamber, the reduction of water to produce H2 takes places, and the oxygen evolution reaction occurs at the anode. The function of the diaphragm is to avoid the mixture of the two gases while maintaining a low resistivity.
Taking as reference a rectangular laboratory electrolysis cell (Figure 3a), the geometry of the model was built following fluid dynamic requirements. In this way, simplifications were made in order to reduce the model complexity ( Figure 3b). As result, a good approximation can be made just by 2D geometry, which allows an optimal study of the main involved variables (Figure 3c). The electrode surfaces (33 mm height) were considered by two simple boundary conditions (1 and 5 in Figure 3c) on both sides of the cell. These domains work as current collectors where potential was applied. Regarding the separator (3 in Figure 3c), it was introduced in the geometry as a thin rectangle (0.5 mm) that separates both compartments. In this case, only ionic conductivity is allowed in this component. Concerning the anodic and cathodic compartments (2 and 4 in Figure 3c), they were defined as two thick rectangles (1.5, 4, or 10 mm depending on the model considered) between the separator and the corresponding electrode with the electrolyte inlet at the bottom and the biphasic mixture (electrolyte and gases) outlet at the top.  Regarding the meshing of the control domain, a rectangular mesh was generated, taking advantage of the regular geometry of the cell, with a total of 7200 elements (when the electrode-diaphragm distance was 10 mm). Previously, a mesh independence study was carried out in order to determine the number of nodes and elements that should be used. In the anode and cathode chamber, a maximum element size of 0.8 mm was established, reducing this value to 0.1 mm in the diaphragm domain (these values in the mesh were scaled when the electrode-diaphragm distance decreased from 10 mm to 1.5 mm). Furthermore, at the boundary of the electrodes and on the surfaces between the diaphragm and the electrolytic chambers, inflation layers were incorporated to correctly model the effects that take place in these contours ( Figure 3d).

Mathematical Procedure and Governing Equations
As previously mentioned, the proposed model is based on previous work carried out by the authors [9] using the CFD simulation program COMSOL Multiphysics v4.3a. In this way, the model has been built using the "electric currents" module to calculate both the current distribution and the polarization curve, and the "two-phase turbulent bubbly flow" module was used to study the generation and distribution of gas, one for each electrolytic chamber. An implicit backward differentiation formula (BDF) method was used for the numerical integration of the differential equations. The tolerance was fixed at 0.001 in the solver to control the absolute error.
The coupling between the electrochemical and fluid dynamic results was linked by the current density. The intensity supplied to the electrolysis cell allows the electric field to be resolved and at the same time determines the generation of hydrogen and oxygen. The gas formation in the cell is studied by means of the conservation equations for the biphasic mixture (hydrogen/electrolyte in cathode and oxygen/electrolyte in anode) that solves, among other variables, the fraction of gas (Φ g ) in each electrolytic chamber. This fraction of gas is used in each instance to determine how the conductivity of the electrolyte varies and again to resolve the electric field inside the cell.

Electrical Current Conservation
The electric charge transport was studied using the "electric currents" module of COMSOL, which solves the current conservation problem (J in A/m 2 ) for the scalar electric potential (U in V). For time dependent studies, the equation takes the following form [24]: where σ (S/m) is the equivalent electrical conductivity, Q (A/m 3 in a 3D approximation) is the current source and J e (A/m 2 ) is the externally generated current density. Regarding the electric displacement it is a constitutive relation that describes the macroscopic properties of the medium. Otherwise, in planar 2D, it is assumed that the model has symmetry and the electric potential varies only in the "x" and "y" directions and is constant in the "z" direction. This implies that the electric field is tangential to the "xy" plane. To solve it, the thickness in the "z" direction (electrode depth) is introduced in the equation by the parameter d (m) [24], as shown in Figure 3a.
On the other hand, to define the equivalent electrical conductivity (σ), the following three equations were incorporated to define the conductivity of the medium (electrolyte) and the diaphragm (see parameters in Table 1): Table 1. Constants, initial, and boundary conditions for the two-dimensional (2D) model [25][26][27][28]. 1. An empirical relationship for specific conductivity of electrolyte (σ 0 in S/m) with respect to temperature (T) and KOH concentration (w) was used [29]:

Symbol
The conductivity of the diaphragm (σ m in S/m) was defined as a function of the conductivity of electrolyte (σ 0 ) and geometric parameters such as porosity (ε m ) and tortuosity (τ m ), according to [28]: 3.
The Bruggeman equation (σ e in S/m) relates the variation of conductivity of electrolyte (σ 0 ) with the volume fraction of gas (Φ g ) inside the cell [30]. The gas fraction for both electrolytic chambers at each current density value was calculated according to the equations described in the "liquid-gas flow distribution" section: Finally, according to Equation (5), the electric transport charges on the electrodes were not considered (only ionic transport was studied inside the electrolysis cell). In order to incorporate the activation overpotentials and the reversible voltage at different temperatures in the polarization curve, the following two auxiliary equations, whose parameters are defined in Table 1, were introduced into the model (as boundary conditions in the electrode): 1.
The reversible potential (U rev in V) was defined according to the LeRoy et al. [31] equation as a function of temperature (T): 2. Activation overpotentials (η act in V) were defined for the cathode and the anode by the Butler-Volmer equation (Tafel equation form) for each current density (i), according to: Liquid-Gas Flow Distribution The gas and liquid flow distribution was modeled, applying the "two-phase turbulent bubbly flow" module of COMSOL. This mathematical approach describes the two-phase flow using an Euler-Euler model. The module solves the volume fraction occupied by each of the two phases, without defining each bubble in detail. This mathematical approach is suitable for modeling the macroscopic behavior of many gas bubbles rising through a liquid. It treats the two phases as interpenetrating media, tracking the averaged concentration of the phases. One velocity field is associated with each phase, and the dynamics of each of the phases are described by a momentum balance equation and a continuity equation [9,32].
Based on these assumptions, the sum of the momentum equations for the two phases gives a momentum equation for the liquid velocity, a continuity equation, and a transport equation for the volume fraction of the gas phase. The momentum equation is [32]: where u is velocity (m/s), p is the pressure (Pa), Φ is the phase volume fraction, ρ is the density (kg/m 3 ), g is the gravity vector (m/s 2 ), F is any additional volume force (N/m 3 ), ψ l is the dynamic viscosity of the liquid (Pa·s), and ψ T is the turbulent viscosity (Pa·s). The subscripts "l" and "g" denote quantities related to the liquid phase and the gas phase, respectively [32]. The continuity equation is: And the transport of the volume fraction of gas is given by Equation (13): where . m gl is the mass transfer rate from gas to liquid (g·m −3 ·s −1 ). In the simplified 2D geometry considered in this study, the gas flows generated (g·m −2 ·s −1 ) on the active surfaces of the electrodes (H 2 and O 2 ) were defined by the Faraday equation for each current density (i): . Table 1 summarizes the parameters corresponding to Equation (14) for hydrogen production (boundary condition at cathode) and to Equation (15) for oxygen production (boundary condition at the anode).
Regarding the turbulence, a k − ε model was used. This model solves two extra transport equations for two additional variables: the turbulent kinetic energy, k (m 2 /s 2 ), and the dissipation rate of turbulent energy, ε (m/s 3 ). The turbulent viscosity (ψ T ), the transport equations for the turbulent kinetic energy (k), and the evolution of the turbulent energy's dissipation (ε) are defined as [32]: where C ε , C µ , Φ k , and Φ ε are constants of the turbulent model. The term S k is related to the bubble-induced turbulence. Finally, the gas velocity ( → u g in m/s) is calculated according to: where → u slip is the relative velocity between the phases and → u dri f t is an additional contribution from the turbulence model used [32].

Initial and Boundary Conditions
To solve the different equations described above, the following simplifications were considered in order to improve convergence and reduce the computational cost of the model [5,9,32]: (1) The gas density is negligible compared to the liquid density; (2) The movement of gas bubbles in relation to the liquid phase is determined by a balance between viscous drag and pressure forces; (3) The two phases share the same pressure field; (4) Hydrogen and oxygen crossover through the diaphragm is negligible; (5) The electrolysis cell works with a high enough flow rate of electrolyte to avoid the accumulation of gas bubbles in the cell; (6) The bubbles have a diameter of less than 1 mm, so the Hadamard-Rybczynski drag law for spherical gas bubbles in liquid is used for the gas velocity; (7) The electrical resistance of the electrodes (Ni) is negligible with respect to the rest of the elements of the electrolysis cell. Figure 4 shows the boundary conditions implemented in this study. Regarding the different materials used in the model, the COMSOL material library was used to define hydrogen, oxygen, and water. In the latter case, the density and dynamic viscosity were also corrected to take into account the KOH dissolution (electrolyte). Table 1 summarizes, the different parameters and constants considered in the model.
where is the relative velocity between the phases and is an additional contribution from the turbulence model used [32].

Initial and Boundary Conditions
To solve the different equations described above, the following simplifications were considered in order to improve convergence and reduce the computational cost of the model [5,9,32]: (1) The gas density is negligible compared to the liquid density; (2) The movement of gas bubbles in relation to the liquid phase is determined by a balance between viscous drag and pressure forces; (3) The two phases share the same pressure field; (4) Hydrogen and oxygen crossover through the diaphragm is negligible; (5) The electrolysis cell works with a high enough flow rate of electrolyte to avoid the accumulation of gas bubbles in the cell; (6) The bubbles have a diameter of less than 1 mm, so the Hadamard-Rybczynski drag law for spherical gas bubbles in liquid is used for the gas velocity; (7) The electrical resistance of the electrodes (Ni) is negligible with respect to the rest of the elements of the electrolysis cell. Figure 4 shows the boundary conditions implemented in this study. Regarding the different materials used in the model, the COMSOL material library was used to define hydrogen, oxygen, and water. In the latter case, the density and dynamic viscosity were also corrected to take into account the KOH dissolution (electrolyte). Table 1 summarizes, the different parameters and constants considered in the model. Regarding the charge transfer coefficient ( ), it is worth noting that 0.5 is the value usually reported in most of the models in which electrochemical kinetics are considered. According to the literature [33], this value is usually used in the absence of actual measurements. However, , in most Regarding the charge transfer coefficient (α), it is worth noting that 0.5 is the value usually reported in most of the models in which electrochemical kinetics are considered. According to the literature [33], this value is usually used in the absence of actual measurements. However, α, in most systems, turns out to be between 0.3 and 0.7 [33]. In this work, as can be seen from Table 1, 0.60 and 0.77 were used for anodic and cathodic reactions, respectively. These values were chosen because they reported the best fit with the experimental values. In the case of the anode, α a is in the common experimental range. On the other hand, the cathodic charge transfer coefficient (α c = 0.77) is slightly higher than the upper limit of the usually experimental range. However, this value well agrees with experimental values reported by authors working on alkaline water electrolysis [34]. In fact, some authors [35] working with real electrolysis stacks have noted that α cannot be considered as a constant, because it varies with temperature, electrode materials, etc., and, thus, it is expected to change even with the experimental system. Furthermore, differences are expected when α is determined in an electrochemical laboratory cell than when it is done in a real water electrolysis cell or even a commercial stack. Figure 5 shows the alkaline water electrolysis test bench used in this study, located in the Centro Nacional del Hidrógeno (www.cnh2.es). The system is constituted by all the elements of a typical electrolyzer, but on a laboratory scale, which allows for the study of a wide range of parameters. Therefore, the system is suitable for analyzing electrochemical and fluid dynamic processes. 0.77 were used for anodic and cathodic reactions, respectively. These values were chosen because they reported the best fit with the experimental values. In the case of the anode, is in the common experimental range. On the other hand, the cathodic charge transfer coefficient ( = 0.77) is slightly higher than the upper limit of the usually experimental range. However, this value well agrees with experimental values reported by authors working on alkaline water electrolysis [34]. In fact, some authors [35] working with real electrolysis stacks have noted that cannot be considered as a constant, because it varies with temperature, electrode materials, etc., and, thus, it is expected to change even with the experimental system. Furthermore, differences are expected when is determined in an electrochemical laboratory cell than when it is done in a real water electrolysis cell or even a commercial stack. Figure 5 shows the alkaline water electrolysis test bench used in this study, located in the Centro Nacional del Hidrógeno (www.cnh2.es). The system is constituted by all the elements of a typical electrolyzer, but on a laboratory scale, which allows for the study of a wide range of parameters. Therefore, the system is suitable for analyzing electrochemical and fluid dynamic processes. As can be seen in Figure 5, the test bench has three electrolysis cells (ELECTROCELL, Micro Flow Cell) that are electrically connected in series, a DC power supply (Elektro-Automatik, EA-PSI 6000), a centrifugal magnetic drive pump (IWAKI, MD-30RVM-220N), two gas-liquid separators to separate the oxygen and hydrogen produced from the electrolyte and a heating system (Hillesheim, H300 DN12). Additionally, the facility has different measurement and control devices that are monitored by a Supervisory Control and Data Acquisition (SCADA) system [10].

Experimental Protocol
The polarization curves were obtained according to a previously established experimental procedure [5,10]. These curves describe the electrochemical behavior of an electrolysis cell and allow determination of the values of voltage and current in which it works. For this purpose, several temperatures (from 30 to 70 °C) were tested. The electrolyte used in the test bench was a concentrated KOH solution in water (from 22 to 42 wt% KOH), and the overall flow rate was fixed at 1.4 l/min per cell. Regarding the electrodes, Ni > 99% was used for the anodes and cathodes. In all cases, the active As can be seen in Figure 5, the test bench has three electrolysis cells (ELECTROCELL, Micro Flow Cell) that are electrically connected in series, a DC power supply (Elektro-Automatik, EA-PSI 6000), a centrifugal magnetic drive pump (IWAKI, MD-30RVM-220N), two gas-liquid separators to separate the oxygen and hydrogen produced from the electrolyte and a heating system (Hillesheim, H300 DN12). Additionally, the facility has different measurement and control devices that are monitored by a Supervisory Control and Data Acquisition (SCADA) system [10].

Experimental Protocol
The polarization curves were obtained according to a previously established experimental procedure [5,10]. These curves describe the electrochemical behavior of an electrolysis cell and allow determination of the values of voltage and current in which it works. For this purpose, several temperatures (from 30 to 70 • C) were tested. The electrolyte used in the test bench was a concentrated KOH solution in water (from 22 to 42 wt% KOH), and the overall flow rate was fixed at 1.4 l/min per cell. Regarding the electrodes, Ni > 99% was used for the anodes and cathodes. In all cases, the active area was 10 cm 2 . The diaphragms used were Zirfon ® Perl 500 UTP (AGFA) with a pore size of 0.15 µm and a thickness of 500 µm [10]. The electrode-diaphragm distance was fixed from 1.5 to 10 mm according to the case considered (see Figure 3b).
For each of the tests, the following protocol was used: Once the cell was assembled, both the electrolyte pump and the heating system were turned on. Once the flow rate and temperature were stable, the power supply was switched on and a gradual current sweep from 50 to 400 mA/cm 2 was performed. For each current intensity value, different operating parameters, such as the cell voltage to obtain the different polarization curves, were monitored and recorded [5,10].

Polarization Curve
The electrochemical response of the model was validated by comparing simulated polarization curves (U model cell ) against those experimentally obtained for the electrolysis cell (U test cell ). The accuracy of the CFD model was also studied by varying some of the most critical operating parameters in an alkaline water electrolysis system. Additionally, the mean relative absolute difference error (MRAD error ) was calculated according to Equation (20) in order to ensure the accuracy and validity of the proposed model [5]: The error calculated was lower than 0.51%, which indicates an excellent correlation between the experimental and modeled results. Figure 6 shows a parity chart between the model and the experimental voltage results. From this figure, it is confirmed that the main discrepancy between the model and real data occurs at low current densities (corresponding to the lowest voltages in the polarization curve). In this zone, activation overpotentials represent the major contribution to the energy required for the electrolysis process. In other words, at these potentials the cathodic and anodic reactions are the limiting steps and, thus, it can be considered that the differences are related to electro-kinetics. In the model, the main electrochemical kinetics parameters involved are i 0 and α. It is well known [35] that both concepts strongly depend on operating conditions (temperature, pressure, reaction, electrode material, etc.). While the influence of i 0 has been widely studied on electrolysis [25,26,35], the charge transfer coefficient has not been treated in much detail. In fact, as described above, in the literature, typically a constant value of 0.5 is given to α for both reactions. In the present work, different values of i 0 , depending on temperature and reaction, have been used (see Table 1). In the case of α, although in the model different values were used for oxidation and reduction, these chosen values were constant for all the operating conditions. The use of these approximations, instead of the real value for each temperature, can be the cause of the discussed deviations at low current densities. Nevertheless, it is worth noting that, despite these differences, the error reported value is very low and the model works satisfactory, as shown in the following sections. Polarization curves simulated by the model were compared with experimental ones in Figure 7, at different temperatures. As expected, the cell potential at a certain current density value decreases as temperature increases within the analyzed range (Figure 7a). This trend is due to the fact that an

Influence of Temperature
Polarization curves simulated by the model were compared with experimental ones in Figure 7, at different temperatures. As expected, the cell potential at a certain current density value decreases as temperature increases within the analyzed range (Figure 7a). This trend is due to the fact that an increase in temperature favors the reaction kinetics, decreasing the reversible voltage and therefore the required energy. This behavior was previously reported by other authors for AWE systems [2,5,10,13,22,23]. From Figure 7b, it is evidenced that the best correlation between experimental and calculated results can be seen at 30 • C, 50 • C, and 70 • C. Figure 6. Parity chart between the model and experimental voltage results for all the conditions of temperature, electrolyte conductivity, and electrode-diaphragm distance analyzed.

Influence of Temperature
Polarization curves simulated by the model were compared with experimental ones in Figure 7, at different temperatures. As expected, the cell potential at a certain current density value decreases as temperature increases within the analyzed range (Figure 7a). This trend is due to the fact that an increase in temperature favors the reaction kinetics, decreasing the reversible voltage and therefore the required energy. This behavior was previously reported by other authors for AWE systems [2,5,10,13,22,23]. From Figure 7b, it is evidenced that the best correlation between experimental and calculated results can be seen at 30 °C, 50 °C, and 70 °C.

Influence of Electrolyte Conductivity
Electrolyte conductivity is a parameter that has a significant effect on the ohmic overpotential. As described from Equation (6) to Equation (8), it is a critical aspect to determine not only the electrolyte ohmic losses, but also those related to diaphragm and generated gas bubbles in both cathodic and anodic compartments. In the case of the diaphragm, it must be taken into account that it is just a physical barrier, without ionic or electronic conductivity, so its conductivity is determined by the conductivity of the electrolyte inside the diaphragm channels. Regarding the gas bubbles, when H2 and O2 are generated, a biphasic mixture gas-electrolyte is formed, in which the gas fraction

Influence of Electrolyte Conductivity
Electrolyte conductivity is a parameter that has a significant effect on the ohmic overpotential. As described from Equation (6) to Equation (8), it is a critical aspect to determine not only the electrolyte ohmic losses, but also those related to diaphragm and generated gas bubbles in both cathodic and anodic compartments. In the case of the diaphragm, it must be taken into account that it is just a physical barrier, without ionic or electronic conductivity, so its conductivity is determined by the conductivity of the electrolyte inside the diaphragm channels. Regarding the gas bubbles, when H 2 and O 2 are generated, a biphasic mixture gas-electrolyte is formed, in which the gas fraction is a non-conductor. The greater the number of gas bubbles in this biphasic mixture, the greater the fraction of gas (Φ g ). As a result, the electrolyte conductivity will be lower according to Equation (8), especially in the close electrode region.
In Figure 8, the polarization curves for different electrolyte concentration values calculated with the model are compared with those experimentally obtained. From this figure, the following points are discussed: (1) As previously reported [10], increasing the electrolyte conductivity implies a better electrolysis performance due to a reduction of ohmic losses, which turns out in a lower required energy. (2) An optimum value of electrolyte concentration is identified at 32 wt% KOH, which corresponds to 94.54 S/m at 50 • C [29]. Above this value, mass transfer limitations can occur, but below it, the ohmic losses are too high because the electrolyte conductivity is very low [29].
(3) The model shows an accurate agreement with the experimental polarization curves for the different electrolyte conductivity values considered (Figure 8b).
performance due to a reduction of ohmic losses, which turns out in a lower required energy. (2) An optimum value of electrolyte concentration is identified at 32%wt KOH, which corresponds to 94.54 S/m at 50 °C [29]. Above this value, mass transfer limitations can occur, but below it, the ohmic losses are too high because the electrolyte conductivity is very low [29].
(3) The model shows an accurate agreement with the experimental polarization curves for the different electrolyte conductivity values considered (Figure 8b).

Influence of Electrode-Diaphragm Distance
In AWE cells, the ohmic overpotential is strongly affected by the proximity of the electrode to the diaphragm. The reason for this is that a reduction of this distance means a lower quantity of electrolyte between both components, and thus the ohmic contribution of this compartment decreases. However, very narrow cells can favor high gas fraction during the electrolysis, which can result in a significant increasing of ohmic overpotential. Hence, it seems reasonable that there is an optimal electrode-diaphragm distance to achieve the best AWE cell performance. In fact, several works have calculated this optimum distance [11,36]. Figure 9 shows the simulated and experimental polarization curves when different electrodediaphragm distances were used, for an electrolyte flow rate of 1.4 l/min. The tendency of the modeled curves well coincide with the real ones; the decrease of the distance between electrode and diaphragm leads to a lower electrolysis potential.

Influence of Electrode-Diaphragm Distance
In AWE cells, the ohmic overpotential is strongly affected by the proximity of the electrode to the diaphragm. The reason for this is that a reduction of this distance means a lower quantity of electrolyte between both components, and thus the ohmic contribution of this compartment decreases. However, very narrow cells can favor high gas fraction during the electrolysis, which can result in a significant increasing of ohmic overpotential. Hence, it seems reasonable that there is an optimal electrode-diaphragm distance to achieve the best AWE cell performance. In fact, several works have calculated this optimum distance [11,36]. Figure 9 shows the simulated and experimental polarization curves when different electrode-diaphragm distances were used, for an electrolyte flow rate of 1.4 l/min. The tendency of the modeled curves well coincide with the real ones; the decrease of the distance between electrode and diaphragm leads to a lower electrolysis potential.

Gas Generation Profile
Evolution of gas profile with current density is shown in Figure 10. Hydrogen and oxygen void fractions (space occupied by gas bubbles in each electrolytic chamber) draw a typical profile next to the electrode surface [17,18,37,38]: void fraction progressively increases from the bottom to the top of the electrode due to the accumulation of generated gas and the effect of the flow rate [15].

Gas Generation Profile
Evolution of gas profile with current density is shown in Figure 10. Hydrogen and oxygen void fractions (space occupied by gas bubbles in each electrolytic chamber) draw a typical profile next to the electrode surface [17,18,37,38]: void fraction progressively increases from the bottom to the top of the electrode due to the accumulation of generated gas and the effect of the flow rate [15].

Gas Generation Profile
Evolution of gas profile with current density is shown in Figure 10. Hydrogen and oxygen void fractions (space occupied by gas bubbles in each electrolytic chamber) draw a typical profile next to the electrode surface [17,18,37,38]: void fraction progressively increases from the bottom to the top of the electrode due to the accumulation of generated gas and the effect of the flow rate [15]. Hydrodynamic characteristics of the biphasic flow (generated gas/electrolyte) in the AWE cell strongly influence the performance of the system. In these systems the generated gas bubbles have a multiple and critical role. On one hand, when they are produced, bubbles form a curtain of increasing thickness along the surface of the electrode in a vertical direction, causing a turbulence phenomenon that contributes to the mixing and distribution of present species and so too to mass transfer. On the other hand, when bubbles are produced, they can adhere to some sections of the electrode surface, Hydrodynamic characteristics of the biphasic flow (generated gas/electrolyte) in the AWE cell strongly influence the performance of the system. In these systems the generated gas bubbles have a multiple and critical role. On one hand, when they are produced, bubbles form a curtain of increasing thickness along the surface of the electrode in a vertical direction, causing a turbulence phenomenon that contributes to the mixing and distribution of present species and so too to mass transfer. On the other hand, when bubbles are produced, they can adhere to some sections of the electrode surface, making them inactive for electrochemical reaction. This causes an increase in the ohmic overpotential, a higher required energy and, as a result, a worse cell performance [19].
According to the previous discussion, the model results show that gas fraction reaches the highest value next to the electrode surface (x = 0 cm) and it decreases towards the diaphragm (x = 0.15 cm) for all studied current densities ( Figure 11). In addition to the expected increase in the fraction of gas with the current, from this graph it can also be deduced that the higher the current density, the greater the penetration of gas in the lateral direction, due to the increase in the lateral gas velocity [14]. making them inactive for electrochemical reaction. This causes an increase in the ohmic overpotential, a higher required energy and, as a result, a worse cell performance [19]. According to the previous discussion, the model results show that gas fraction reaches the highest value next to the electrode surface (x = 0 cm) and it decreases towards the diaphragm (x = 0.15 cm) for all studied current densities ( Figure 11). In addition to the expected increase in the fraction of gas with the current, from this graph it can also be deduced that the higher the current density, the greater the penetration of gas in the lateral direction, due to the increase in the lateral gas velocity [14].

Conclusions and Future Work
In the present work, the development and experimental validation of an alkaline water electrolysis cell 2D CFD model was reported. The strong points of the model are its versatility to simulate different operating conditions as well as the simultaneous analysis of fluid dynamic and electrochemical phenomena. These aspects make it a powerful and cheap design tool, particularly

Conclusions and Future Work
In the present work, the development and experimental validation of an alkaline water electrolysis cell 2D CFD model was reported. The strong points of the model are its versatility to simulate different operating conditions as well as the simultaneous analysis of fluid dynamic and electrochemical phenomena. These aspects make it a powerful and cheap design tool, particularly suitable for the study of new cell configurations.
Regarding the variation of electrolysis parameters, the response of the model was evaluated in terms of polarization curves at different values of temperature, electrode-diaphragm distance, and electrolyte concentration. Performance of the electrolysis cell was improved for high temperatures, small electrode-diaphragm distances, and high electrolyte conductivity, as observed from both experimental and simulated data. In this regard, the mean relative absolute difference error was calculated in order to ensure the accuracy and validity of the proposed model. The error calculated was lower than 0.51% for the polarization curve, which indicates an excellent correlation between the experimental and modeled results. Therefore, the proposed model is suitable for predicting the behavior of the AWE cell under a wide range of operating conditions.
In addition to electrochemical studies, the proposed model also offers a fluid dynamic analysis. From simulation results, it was concluded that produced gases generate a "curtain profile" in the electrode, increasing the void fraction in a vertical direction due to the accumulation of gas. Furthermore, the model simulates the gas profile along the x-axis of the electrolytic chamber (from the electrode to the diaphragm), proving that the higher the current, the greater the gas distribution in the electrolytic chamber space. These results provide very useful information for cell design purposes.
On the other hand, although the results confirm that the model well fits the experimental results, various points still need to be studied: (1) Charge transfer coefficient (α): the influence of this term over the performance of real electrolysis systems has hardly been studied. Even though it is possible to find different values in the literature, some authors have claimed significant differences when working with a laboratory electrochemical cell and with an electrolysis stack. Experimental values up to 1.035 [35] have been reported, which contrast with that which is considered to be the common range (0.3-0.7). In fact, it greatly depends on operating conditions: electrode materials, electrolyte, etc. For these reasons, it would be very interesting to experimentally determine this parameter in the test bench shown in Figure 5 in order to improve the accuracy of the model, especially in the activation overpotentials zone where the model fits worse with real data. This is extensible to exchange current density, which is another critical kinetic parameter. (2) Zero-gap cell design: this is one of the most extended alkaline water electrolysis cell designs in commercial systems. This configuration allows us to "ideally" eliminate bubbles between the electrodes because cathode and anode are placed directly over the diaphragm. For this to happen, the electrodes have to be porous. Therefore, in order to improve the model applicability to real systems, this geometry should be implemented in the model, making it a more useful design tool.

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, and in the decision to publish the results.