# Coupling morphological evolution and multiphase transport to Decipher degradation in PEM water electrolyzers

Youssef Kharchouf, Loïc Assaud, Mireille Turmine, Vincent Vivier

► To cite this version:

Youssef Kharchouf, Loïc Assaud, Mireille Turmine, Vincent Vivier. Coupling morphological evolution and multiphase transport to Decipher degradation in PEM water electrolyzers. Chemical Engineering Science, 2026, 327, pp.123687. 10.1016/j.ces.2026.123687. hal-05533902

HAL Id: hal-05533902

https://hal.science/hal-05533902v1

Submitted on 3 Mar 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Coupling morphological evolution and multiphase transport to Decipher degradation in PEM water electrolyzers

Youssef Kharchouf<sup>a, ∗</sup>, Loïc Assaud<sup>b, ∗</sup>, Mireille Turmine<sup>a</sup>, Vincent Vivier<sup>a, †</sup>

ARTICLE INFO

Keywords:

PEM Electrolysis

Degradation

Porous Transport Layer

Multiphase flow

Level set method

A B S T R A C T

In this work, we present a multiphase numerical model for the corrosion degradation of proton exchange membrane water electrolyzers. The numerical framework comprises two main components, a geometry tracking component and a performance evaluation component. The geometry tracking is based on the Level Set method that models the evolution of the geometry due to the deposition or erosion of corrosion products in the porous transport layer along the lifecycle of the device given a specified corrosion rate. The performance of the device at each time step is evaluated using a multiphase model with the Navier-Stokes equation for the two-phase mixture, and a transport equation for a dispersed phase of bubbles in the continuous water phase. Catalyst layer degradation is included through appropriately chosen boundary conditions in the geometry. This integrated approach allows the model to directly link structural changes in the porous transport layer to evolving transport phenomena and the overall electrochemical performance of the device. The corrosion process is introduced through an effective corrosion velocity, which provides a phenomenological description of PTL degradation rather than a mechanistic prediction based on local electrochemical conditions. Thus, the framework should be considered a multiphysics tool that propagates assumed degradation scenarios to device-scale performance rather than a first-principles corrosion model.

1. Introduction

Hydrogen is highly regarded as a clean and promising energy carrier, offering both an effective means of energy storage and a solution to the intermittency challenges of renewable energies. With the current growing contribution of renewable sources such as solar, wind and biomass, hydrogen can serve as an efficient carrier to enhance their overall effectiveness. However, the majority of hydrogen is still produced by steam methane reforming (Abdin, 2020), a process that relies on fossil fuels. This paves the way for considerable research to go into direct electrochemical water splitting, which is more straightforward but it is energy intensive and currently constitutes only about 4% of global hydrogen production (Ji and Wang, 2021).

Despite alkaline water electrolysis (AWE) being currently the more
commercially mature technology (Ursua, 2012), proton exchange
membrane water electrolyzers (PEMWE) offer a number of advantages,
such as higher current density, energy efficiency and a wider range of
operating temperatures (Wang, 2022; El-Shafie, 2023). Their quick
response time under dynamic loads also makes them more suitable for

pairing up with intermittent renewable sources (Zhang, 2022) when
compared to AWE. Despite these advantages, wide adoption of PEM
electrolysis is hampered by the cost of its components and electro-
catalysts, as well as their degradation (Wang, 2022; Zhang, 2022; Nor-
azahar, 2024; Prestat, 2023), which makes the investigation of the
mechanisms and effects of the degradation of these devices an important
part of the research literature on PEMWES (Norazahar, 2024; Grigoriev,
2014; Omrani and Shabani, 2021; Siracusano, 2018; Papakonstantinou,
2020). In particular, the porous transport layer (PTL) can represent a
significant portion of the capital expenditure associated with this tech-
nology due to the need for titanium along with expensive noble metal
coatings (Liu, 2018). In fact, the PTL and bipolar plates (BPPs) account
for up to 50% of the overall stack cost (Yuan, 2022). Most of the previous
studies in the literature, however, evaluate degradation at the device
scale solely through changes in cell voltage over the lifetime, providing
limited insight into the relative contributions of individual degradation
mechanisms to overall performance (Suermann, 2019). Unlike fuel cells,
there is also a lack of a common standard for accelerated stress testing
these devices (Babic, 2017), making it difficult to compare experimental

<footnote><sup>a</sup> Sorbonne Université, CNRS, Laboratoire de Réactivité de Surface (UMR 7197), 4 place Jussieu, 75005 Paris, France</footnote>

<footnote><sup>b</sup> IMAP, UMR8004, University Paris Sciences & Letters, ESPCI, CNRS, ENS, 10 rue Vauquelin, 75005 Paris, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: vincent.vivier@sorbonne-universite.fr (V. Vivier).</footnote>
data across different studies.

A large number of works in the literature used mathematical models and numerical simulation to study different aspects of the PEMWE (Arbabi, 2016; Paliwal, 2021; Tijani and Rahim, 2016; Falcao and Pinto, 2020). For instance, Arbabi et al. (Arbabi, 2016) introduced a volume of fluid (VoF) based model to investigate oxygen bubble penetration in the PTL. They performed experimental validation using a Poly-dimethylsiloxane (PDMS) microfluidic platform that emulates the geometry of a PTL, along with ethanol emulating water as it provides the correct contact angle with PDMS as water with titanium (Arbabi, 2014). More recently, Paliwal et al. (Paliwal, 2021) constructed a similar model but using the Lattice Boltzmann Method with the Shan-Chen technique to incorporate multiphase flow. One-dimensional models are also widely used as they enable the inclusion of multiple physical phenomena without the computational constraints typically associated with complex two- or three-dimensional geometries. Among these models is the work by Zeng et al. (Zeng, 2024) whose 1D model includes heat, mass and charge transport, and was used to study the effect of the operating temperature and pressure on the mass transport and electrolytic performance of the cell.

However, few have directly addressed degradation modeling, especially for the degradation of the PTLS. In fact, a large class of models rely on an empirical or semi-empirical relationship between the voltage and current (Harrison, 2005), and using the Nernst equation to compute the open circuit voltage (Ruuskanen, 2017). Notably, Chandesris et al. used a one-dimensional transport model to study the effect of oxygen crossover on the degradation of the membrane degradation, including the effects of temperature and current density (Chandesris, 2015). In the same vein, Norazahar et al. performed a reliability analysis on PEM water electrolyzers, but used a Bayesian model and not a physics based transport or electrochemical model (Norazahar, 2024). Such approaches rely on relatively large experimental datasets in order to infer correlations between observable parameters such as temperature and current density, and potential failure modes. While valuable for operational monitoring, they do not resolve the detailed multiphase transport and electrochemical phenomena within the device, limiting their utility for design optimization or predictive evaluation.

More recently, machine learning approaches have also been explored to investigate and optimize PEMWEs, offering data-driven alternatives to capturing performance trends beyond traditional models. Lagela et al. (Legala, 2022) used an artificial neural network (ANN) along with a Support Vector Machine (SVM) to model a PEM Fuel cell, using cell current, temperature, pressure and other operating conditions as model inputs to predict the cell voltage and membrane resistance. Similarly, Madhavan et al. (Madhavan, 2025) compared Extreme Gradient Boosting (XGBoost) and ANN models in predicting the corrosion resistance of various PTL coatings using experimental datasets of corrosion currents and end-of-life voltages. However, these data driven approaches require sizeable, well-constructed, and balanced datasets in order to maximize the accuracy of the models, which implies extensive and systematic experimentation that can be time-consuming, especially when it comes to studying aging and degradation phenomena. The objective of this work is not to predict corrosion kinetics from local electrochemical variables, but rather, to link prescribed morphological evolution scenarios with the associated changes in transport and electrochemical performance. The PTL is represented by a simplified, two-dimensional (2D), pore-scale cross-section that acts as a representative elementary volume (REV), consistent with previous 2D studies of multiphase transport in PEMWE cells. This deliberate choice enables us to focus on the coupling between geometric evolution and two-phase transport.

In this work, we propose a numerical model that accounts for the degradation of a PEMWE due specifically to corrosion. The model is composed of two main components: a geometry tracking component and the performance evaluation component. The geometry tracking model is based on the Level Set Method, in which the initial porous geometry of

the PTL evolves as a function of the deposition or dissolution of corrosion products, which will affect the transport properties inside the PTL. This component provides snapshots of the state of the geometry at different periods over the lifetime of the device. These snapshots are then used as the simulation domain by the performance evaluation component, which is a multiphase model that solves for the volume fraction occupied by the oxygen bubbles in the liquid water phase within the porous transport layer. In this component, the oxygen evolution reaction (OER) is implemented as a boundary condition at the PTL interface representing the catalytic layer.

While a large number of existing studies quantify degradation solely through empirical or semi-empirical correlations between cell voltage and operating time, such approaches provide little insight into the contributions of individual mechanisms such as catalyst degradation, bubble-induced transport losses, or pore blockage due to corrosion. The novelty of this study lies in the explicit coupling of the geometry evolution with performance evaluation within a multiphysics framework. Geometry evolution is treated with the Level Set method, which allows for dynamic tracking of corrosion-induced deposition and topology changes in the porous medium. This evolving geometry is then directly linked to a two-phase transport model that evaluates the effect of structural changes on bubble dynamics and mass transport at the catalyst layer (CL). This physics-based degradation model therefore bridges the gap between device-scale performance data (e.g., polarization degradation curves) and the local morphological and physicochemical changes occurring in the PTL. This model has been validated against experimental data from the literature, showing close agreement with polarization curves in the pristine (no degradation) regime while also capturing local phenomena such as bubble coverage in the device. Furthermore, it is capable of reproducing observed voltage variation trends under aging conditions, confirming its ability to connect device-scale performance metrics with the underlying multiphase transport and morphological changes that drive degradation. The governing equations were solved using the finite element method with the COMSOL Multiphysics software (version 6.2, COMSOL Multiphysics) while the geometry generation and analysis was performed with in-house python code.

## 2. Description of the PTL degradation model

### III. The Level Set Method.

The Level Set Method (LSM) is a numerical technique for modeling evolving interfaces on a Cartesian grid or a fixed mesh, eliminating the need for explicit interface parameterization, and allowing for arbitrary topological changes (Sethian, 1996).

The core idea behind LSM is to represent a higher-dimensional surface using a function $\psi$, whose level set is the interface itself (Fig. 1). Assuming the initial interface is defined by the curve $\Gamma(t = 0)$, the corresponding higher-dimensional surface can therefore be expressed as follows:

$$ \psi(\mathbf{x}, t) = \pm d \qquad (1) $$

Where $\mathbf{x} \in \mathbb{R}^n$, and n is defined by the dimensionality of the interface as $n = n_{int} + 1$, and d is the distance from $\Gamma(t = 0)$, with the sign indicating one of the two regions that are separated by $\Gamma(t = 0)$. The function $\psi$ can therefore be initialized as:

$$ \Gamma(t = 0) = [\mathbf{x} | \psi(\mathbf{x}, t) = 0] \qquad (2) $$

Once $\psi$ has been initialized, instead of finding an evolution mechanism for the interface $\Gamma$, an equation describing the evolution of $\psi$ is used, with the condition that the interface remains defined by the zero level set of $\psi$. The advantage of such a description is that the level set may change, merge or split etc. (i.e. corresponding to a change in its topology) while $\psi$ itself evolves smoothly in a fixed mesh. It should be noted that the normal vector (n) to the interface can also be obtained from the normalized gradient of $\psi$ on the interface:
<figure><img src="image_3.png" /><figcaption>Fig. 1. Diagram of a two-dimensional surface plot of ψ(x) and a contour of its zero-level set indicating the interface. A smooth evolution of the function ψ(x) can lead to topological changes in the level set. (a) ψ(x) before split. (b) ψ(x) after split.</figcaption></figure>

$$n = \frac{\nabla \psi}{|\nabla \psi|} \qquad (3)$$

In this work, the evolution of the interface will be controlled by an underlying velocity field **u**. Thus, in order to move the interface according this velocity field, a classical convection formula can be used, with added numerical stabilization terms on the right-hand side:

$$\frac{\partial \psi}{\partial t} + u \cdot \nabla \psi = \gamma \nabla \cdot (\epsilon \nabla \psi - \psi(1 - \psi)\mathbf{n}) \qquad (4)$$

Where γ is the reinitialization parameter, which stabilizes the convection of the field ψ by constraining its gradients to the interface zone, the thickness of which is defined by ε. For clarity, the full sequence of operations involved in the level-set update is summarized below.

(i) The PTL domain is initialized as a signed-distance field, ψ(x, t = 0);

(ii) A prescribed corrosion velocity field, u<sub>corr</sub>, normal to the interface, is applied using Eq. (5);

(iii) The convection–reinitialization Eq. (4) is solved using a stabilized finite-volume method;

(iv) The zero-level set of psi defines the updated PTL geometry for the next time step.

The evolution of the geometry results from the deposition or erosion of corrosion products over the entire PTL. However, the deposition rate will depend on a number of factors that are tied to corrosion kinetics and the local operating conditions. We therefore define a corrosion rate variable u<sub>corr</sub> that encapsulates all these factors. The underlying velocity field that we use is obtained by the following formula:

$$u = u_{\text{corr}} \psi(1 - \psi) \cdot n \qquad (5)$$

A positive u<sub>corr</sub> represents the deposition of products on the PTL over time, while negative values would represent their erosion. However, since we use this variable to encapsulate multiple degradation mechanisms, we restrict this work to the case where u<sub>corr</sub> is positive as it leads

to the progressive reduction of mass transport in the PTL. A more comprehensive model which would be able to represent the erosion case would need to also include electronic transport and the evolution of electrical resistance in order to accurately capture the degradation tendencies of the device. As a result, the geometry is initialized by setting the function ψ(x, t = 0) = ±1 in the bulk of the solid and fluid regions, and by equation (1) to interpolate smoothly in the interface. The function ψ then evolves through the advection equation (4) indicating the shape of the geometry at each time step through its zero-level set. The reinitialization parameter γ is set to equal the corrosion rate u<sub>corr</sub> while the interface thickness is dynamically set according to the mesh element size. These solutions are subsequently used by the second component of the model to construct the appropriate simulation domain at each time step in the lifetime of the device.

The level set formulation adopted in this study allows for the continuous tracking of corrosion-driven interface evolution. This includes pore throat constriction, interface merging, and changes in the topology of the porous transport layer. Unlike simplified approaches based on resizing idealized particles directly, this method captures connectivity loss and progressive pathway blockage, which strongly influence multiphase transport. Additionally, future extensions of the framework can easily incorporate spatially and temporally varying corrosion velocities without altering the numerical structure of the model.

### 3. Performance evaluation model

Once the geometry is obtained at a given time step, we proceed to the multiphase model in order to evaluate the performance of the device. It is based on a momentum transport model of a dispersed phase (oxygen bubbles) in a continuous phase (liquid water). The density and viscosity of the mixture are given by

$$\rho = \phi_c \rho_c + \phi_d \rho_d \qquad (6)$$

$$\mu = \phi_d \mu_d + \phi_c \mu_c \qquad (7)$$

Where φ<sub>c</sub> and φ<sub>d</sub> are the volume fraction of the continuous and dispersed phase, respectively, while ρ<sub>c</sub> and ρ<sub>d</sub> are their corresponding densities. The transport equation of the dispersed phase is given as:

$$\frac{\partial \phi_d}{\partial t} + u_d \cdot \nabla \phi_d = 0 \qquad (8)$$

And the momentum conservation is given by:

$$\rho \frac{\partial \mathbf{u}}{\partial t} + \rho(\mathbf{u} \cdot \nabla)\mathbf{u} + (\rho_d - \rho_c)(\mathbf{u}_{\text{slip}} \cdot \nabla)\mathbf{u} = -\nabla p + \nabla \cdot \mathbf{K} + \rho g \qquad (9)$$

In this equation, K represents the viscous stress tensor, and g the gravitational acceleration, which will significantly affect the flow in this case as the continuous phase density is about three orders of magnitude higher than the dispersed phase density. The slip velocity between the two phases is determined by the drag force acting on the bubbles, which is quantified using the widely adopted Schiller-Naumann correlation (Guan, 2020; Karunarathne and Tokheim, SIMS(2017)) which is defined by the following formula for the drag coefficient:

$$\left\{ \begin{array}{l} C_D = \frac{24}{\text{Re}_b} (1 + 0.15\text{Re}_b^{0.687}), \text{Re}_b \le 1000 \\ C_D = 0.44, \text{Re}_b > 1000 \end{array} \right. \qquad (10)$$

Here, C<sub>D</sub> is the drag coefficient, which quantifies the resistance experienced by a dispersed bubble in the continuous liquid phase. The term Re<sub>b</sub> = ρ<sub>c</sub>d<sub>b</sub>|u<sub>c</sub> - u<sub>d</sub>|/μ<sub>c</sub> is the bubble Reynolds number, where ρ<sub>c</sub> and μ<sub>c</sub> are the density and dynamic viscosity of the continuous phase, d<sub>b</sub> is the bubble diameter, and |u<sub>c</sub> - u<sub>d</sub>| the magnitude of the slip velocity between the
continuous and disperse phases. The drag coefficient is then used to evaluate the interphase momentum transfer, which determines the actual slip velocity of the bubbles within the liquid phase. The multi-phase formulation used in this work follows a dispersed-phase mixture approach, which allows for efficient simulation of two-phase transport across the entire porous transport layer and flow channel domain. Although capillary pressure effects dominate pore-scale liquid-gas in-vision processes, resolving liquid-gas interfaces explicitly through separated-interface methods, such as VOF, LBM, or phase-field models, requires a spatial resolution that is too fine for device-scale simulations. Thus, the present reduced-order formulation prioritizes capturing large-scale transport trends and their coupling to evolving porous geometry rather than detailed pore-scale capillary dynamics.

The boundary conditions that are applied are no-slip on the walls $u = 0$, and two inlet boundary conditions for water inlet and the oxygen evolution reaction at the catalytic layer. In the case of the water inlet, the flow rate is guaranteed by the imposed velocity, as well as a dispersed phase volume fraction equal to zero, which corresponds to the fact that only liquid water is entering the simulation domain from this boundary.

$$ \mathbf{u}_{\text{inlet}} = \frac{Q}{S_{\text{inlet}}}, \phi_d = 0 \quad (11) $$

Where $Q$ is the desired flowrate and $S_{\text{inlet}}$ is the surface area of the inlet. However, for the catalytic layer, the velocity field is composed of both the dispersed and continuous components. Conceptually, the continuous phase is being consumed while the disperse phase is generated at the catalytic layer, according to the reaction:

$$ 2\text{H}_2\text{O}_{(l)} \rightarrow \text{O}_{2(g)} + 4\text{H}^+ + 4\text{e}^- \quad (12) $$

$$ \mathbf{u}_{\text{CL}} = \frac{\mathbf{j}_{\text{O}_2} + \mathbf{j}_{\text{H}_2\text{O}}}{\rho_d + \rho_c} \quad (13) $$

Where $\mathbf{j}_{\text{O}_2}$ and $\mathbf{j}_{\text{H}_2\text{O}}$ are the mass fluxes of the disperse and continuous phases, respectively. The mass fluxes of the two phases are not independent, and in order to ensure mass conservation during the reaction the following formula is used:

$$ \mathbf{j}_{\text{O}_2} = -\mathbf{j}_{\text{H}_2\text{O}} \cdot \frac{2\text{M}_{\text{H}_2\text{O}}}{\text{M}_{\text{O}_2}} \quad (14) $$

Where $\text{M}_{\text{H}_2\text{O}}$ and $\text{M}_{\text{O}_2}$ are the molar masses of the $\text{H}_2\text{O}$ and $\text{O}_2$, respectively.

The cell voltage that is necessary across the device for water-splitting is composed of multiple reversible, activate, ohmic and transport components. Since our model is focused on the anode side, we assume a constant voltage drop across the cathode and model the various over-potentials on the anode side. For instance, the transport overpotential that results from oxygen bubbles partially blocking access to the CL surface is modelled by a Tafel equation:

$$ \eta_{\text{bub}} = A \ln \left( \frac{i}{i_0 (1 - \phi_{\text{CL}})} \right) \quad (15) $$

Where $A$ is the Tafel slope, and $i$ and $i_0$ are the current and exchange current densities, respectively, and $\phi_{\text{CL}}$ is the volume fraction of the oxygen phase at the catalyst layer, representing the surface that is blocked by bubbles.

The total cell voltage is therefore the sum of the contributions of the different loss mechanisms:

$$ V = E_{\text{rev}} + \eta_{\text{bub}} + \eta_{\text{act}} + \eta_{\text{eat}} + \eta_{\text{ohm}} \quad (16) $$

## 4. Simulation domain

In order to achieve a compromise between model accuracy and

computational efficiency, a two-dimensional cross-sectional domain of the anode porous transport layer was employed as the simulation domain. The cross-section spans from the CL-PEM interface to the BPP across the flow channel and was taken to be representative of an in-plane vertical slice through the PTL microstructure, as shown in Fig. 2.

We assume that out-of-plane gradients over the velocity field are negligible so that the essential features of the oxygen-water two-phase transport (nucleation of bubbles at the CL-PEM interface, bubble rise through the porous microstructure, and evacuation towards the flow channels) are captured without resorting to a fully resolved 3D domain. This assumption is supported by the fact that oxygen content showed periodicity at the order of 100 µm, as well as the fact that oxygen tends to take preferential pathways irrespective of the water flow rate or current density (Satjaritanun, et al., 2020). Furthermore, several modeling works in literature have adopted 2D cross-sections of the device, as continuum geometric domains or pore-network structures, in order to study two-phase transport or electrochemical performance, demonstrating that the dominant physics of the flow can be successfully captured with a 2D model (Altaf, 2020; Chen, 2020). In this study, the 2D pore geometry is understood as a reduced-order surrogate that captures essential pore-scale features relevant to bubble transport. Although this approach cannot reproduce the full 3D tortuosity or connectivity of commercial PTLs, it allows for the coupling of geometry evolution at an affordable computational cost. Additionally, it provides a foundation that can be replaced in future work with µCT-derived 3D structures without altering the model framework. Similarly, experimental two-dimensional microfluidic platforms were used to visualize two-phase transport that is representative of a PEM water electrolyzer, with the porous characteristics such as pore size, porosity or tortuosity being obtained from the actual PTLs using micro-computed tomography (Arbabi, 2016; Arbabi, 2014). In this case, the model can handle different studies on porosity or tortuosity as long as the geometry generation method is compatible. While this work focuses mainly on porosity as the main study parameter, more sophisticated micro-porous structure generation methods include the effect of tortuosity, and the resulting generated structures can directly be used with the current model without any modifications to the model itself. Furthermore, the in-plane extent is chosen as a representative elementary volume, large enough to include multiple pores and throats, while still remaining at a spatial scale on the order of $10^2$ µm that was experimentally observed with Operando X-ray computed tomography (Satjaritanun, et al., 2020).

The computational domain was discretized using a hybrid unstructured mesh which is composed of triangular elements in the bulk and progressively refined quadrilateral elements in the boundary layers. As the geometry evolves with degradation, the mesh size increases from approximately 130,000 elements to 150,000 in order to capture the progressive narrowing of pore throats as corrosion products constrict them (Fig. 3). The mesh convergence analysis was carried out on an aged porous transport layer (PTL), in which constricted pore throats are much more prevalent. Such geometrical features require a finer spatial resolution to be accurately captured, leading to a higher number of mesh elements compared to the pristine configuration. As shown in the convergence study (Fig. 3), increasing the mesh density beyond approximately 150 000 elements does not induce any noticeable variation in the computed overpotential, indicating that the numerical solution is mesh-independent at this resolution.

The velocity, pressure, and gas volume fraction equations were discretized using linear Lagrange finite elements. Furthermore, numerical stabilization was enabled using COMSOL's inconsistent stabilization with isotropic diffusion to ensure a stable incompressible-flow solution and to prevent spurious oscillations in the transported volume fraction field. The resulting linear systems for velocity, pressure, and volume fraction were solved using the direct PARDISO solver. The computations were carried out on a standard scientific workstation. The boundary conditions are no-slip conditions on the external walls of the domain and the fibers everywhere except at the CL, inlet and outlet. At the CL,
<figure><img src="image_4.png" /><figcaption>Fig. 2. Example of the 2D simulation domain used. The porous transport layer (PTL) fibers are represented as red disks, while the interconnected fluid region (pore space) is shown in blue. The catalyst layer (CL) is located at the bottom boundary, with the fluid inlet on the right and outlet on the left.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Mesh convergence graph using the cell overpotential as the monitored quantity, plotted against the total number of mesh elements. The stabilization of the overpotential indicates mesh-independent results.</figcaption></figure>

oxygen generation is prescribed by a flux boundary condition derived from the imposed current density (from equation (13), as well as water availability through the water volume fraction. It is important to note that, because the fluxes of oxygen $j_{O_2}$ and water $j_{H_2O}$ have opposite signs, the catalyst layer boundary condition corresponds to an inflow of oxygen and an outflow of water. The inlet velocity is prescribed based on the imposed flow rate and inlet length, and applies only to the liquid phase. A pressure boundary condition is applied at the outlet for both phases.

## 5. Results and discussion

### VII. Pristine regime.

Before addressing the degradation aspects of the model, we first conduct a validation step using the polarization curves generated by the model under degradation-free conditions. For this purpose, we used experimental measurements reported by Liu et al. (Liu, 2018), who

investigated the influence of iridium oxide coatings on the resistance of the PTL-CL interface and overall device performance. This dataset is particularly suitable because it provides IR-Free polarization curves, which are preferable for comparison with our model. The necessity of IR-Free curves arises from the fact that the model in this work does not explicitly solve for electric charge transport within the solid, electronically conducting phases of the device, namely the titanium phase in the PTL and BPPs that act as electrical contacts. In other words, the model developed in this work does not account for any electrical, contact or electronic resistances within the PTL, but does account for its hydraulic resistance.

In the experiments of Liu et al. (Liu, 2018) titanium PTLs with a thickness of 350 µm and a porosity of approximately 68% were used with fiber diameters ranging between 15 and 50 µm. The cells were operated at a temperature of 80°C and a pressure of 1 bar. Fig. 4 compares the experimental polarization curves under these conditions with
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Comparison of experimental polarization curves of PEM electrolyzer cells with titanium PTLs and model predictions (Liu, 2018).</figcaption></figure>

the predictions of our model. The primary comparison is made with IR-Free data. For completeness, an IR-included curve was also generated by introducing a single fitted resistance parameter. This fitted value, approximately 218 mΩ·cm², yields a good match to the experimental IR-included curve and is in close agreement with the internal resistance of 220 mΩ·cm² measured by electrochemical impedance spectroscopy (EIS) in Liu et al. (Liu, 2018), serving as a useful consistency check.

Fig. 5a illustrates the predicted bubble coverage within different regions of the PEMWE geometry. In the PTL region, bubble coverage increases linearly with current density, which reflects independence from water flow in the BPP. In contrast, the flow channels are strongly influenced by convective flow, which evacuates the oxygen bubbles even at high current densities, explaining the deviation of the flow channel curves from linearity. It should be noted that the approximately linear increase in bubble coverage with current density predicted here reflects the simplified geometry of the model and absence of explicit coalescence/detachment kinetics. Therefore, the linear trend in Fig. 4a should be interpreted as a qualitative response to the simplified pore morphology rather than a quantitative prediction of PTL saturation under industrial conditions.

To further validate the bubble coverage predictions, we compared the model output with experimental observations reported by Sinapan et al. (Sinapan, 2024), who used the convolutional neural network YOLOV7 on images of the flow channels on the anode side of the PEMWE to detect, classify and analyze the generated oxygen bubbles and their coverage. Their setup measured bubble coverage of the flow channels at current densities between 0.3 and 1.4 A/cm², given a water flow rate of 0.4 L/min at 60°C and ambient pressure. The PTL was a 180 µm titanium grid. In our model the flow rate was normalized by the experimentally reported active surface area of 40.71 cm².

Fig. 5b presents a direct comparison between the experimental bubble coverage data of Sinapan et al. (Sinapan, 2024) and the predictions of our model under equivalent operating conditions. The agreement is excellent, with deviations remaining around 2% throughout most of the tested current densities. In fact, the error rises to 4% at higher current densities where the model's accuracy decreases due to the increasing prevalence of bubble coalescence, which is not taken into account. This close match confirms that the multiphase

formulation captures the essential physics of bubble accumulation and evacuation in the flow channels. The distribution of the volume fraction over the PTL region that is predicted by the model is shown in Fig. 6.

The close agreement between simulation and experiment supports the validity of the model in the pristine regime, thereby enabling us to proceed to the next phase of this work, in which degradation effects are explicitly considered. This strong agreement lends confidence that the framework accurately captures the baseline physics, ensuring that any subsequent deviations can be unambiguously attributed to degradation phenomena.

## 6. Degradation mechanisms

In this section we describe how degradation mechanisms are included in the model. The first aspect of degradation is introduced through the geometry evolution, which impacts mass transport and leads to an increase in mass transport overpotential at the CL, due to the coupling of the active surface area to the volume fraction of the dispersed phase in the multiphase model. The evolution of the corrosion-driven geometry is governed by an effective corrosion velocity, which is the net result of multiple degradation processes, such as titanium dissolution, oxide deposition, and local chemical effects. It is important to emphasize that the corrosion rate used in level set evolution is an effective phenomenological parameter. It does not arise from local electrochemical kinetics, such as titanium (Ti) dissolution rates, potential distribution, pH, or species transport. These factors would require a fully coupled electrochemical corrosion model. Instead, the corrosion rate aggregates the net effect of these processes into an assumed geometry-evolution rate. In the present study, this rate is treated as constant, corresponding to an accelerated aging scenario. It does not arise from a mechanistic electrochemical corrosion law based on local potential or species concentrations. This approach aligns with previous degradation and reliability studies that use empirical or semi-empirical descriptions of PTL or membrane aging. Therefore, the purpose of the present framework is to propagate prescribed degradation scenarios to the device scale, rather than to mechanistically predict corrosion initiation or propagation. The catalyst degradation is introduced independently through a progressive reduction of the exchange current density
<figure><img src="image_7.png" /><figcaption>Fig. 5. (a) Bubble coverage in various components of the electrolyzer. The coverage is linear with current density at the PTL area but not at the outlet or BPP channel where incoming water evacuates the bubbles; (b) Comparison of predicted bubble coverage in the electrode flow channels with experimental measurements by Sinapan et al. (Sinapan, 2024) at 60°C, ambient pressure, and a water flow rate of 0.4 L/min.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Comparison of the distribution of the volume fraction of oxygen in the PTL at different current densities. (a) 1.4 A/cm²; (b) 1 A/cm²; (c) 0.3 A/cm².</figcaption></figure>

at the catalyst layer.

This model also reflects the fact that when gas bubbles are generated, they momentarily block access to the CL and part of the reaction surface, until the bubbles are evacuated towards the flow channel and outlet. Given that the corrosion velocity is considered to be constant, the change in the effective porosity of the porous transport layer is almost linear, with a slight deviation caused by random regions of the porous

space being isolated from the rest of the geometry as the pathways leading to them are blocked by the growth of the solid region, as shown in Fig. 7.

Another degradation factor that the model takes into consideration is the loss of catalyst performance. In real PEMWE systems, PTL corrosion, Ir dissolution and redeposition, loss of contact resistance, and catalyst inactivation are tightly coupled processes. The present model treats PTL geometry evolution and catalyst degradation as independent inputs to isolate their respective impacts on performance. A coupled treatment of Ti dissolution, Ir redeposition, CL/PTL interfacial resistance, and electronic contact degradation would require additional reaction-transport models, which are beyond the scope of this initial implementation but would be a logical next step.

Since the reaction is modelled as a boundary condition, we model the catalyst's evolution as a reduction in the available active surface area for the oxygen evolution reaction.

The effect on the overpotential is shown in Fig. 8a, highlighting that most of the performance loss occurs during the final steps once degradation exceeds 50%. It should be noted however, that most of the reduction in performance that occurs in real electrolyzer devices is mainly due to transport losses and ohmic losses, especially at higher current densities (Suermann, 2015).

## 7. Degradation regime

In this section, we validate the model by comparing its performance
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** (a) Example porous domain under initial conditions, where the black regions represent the solid phase, and the white regions the fluid domain; (b) Simulation domain after 150 h; (c) Evolution of the effective porosity, which is mainly linear due to the constant corrosion rate, with slight deviations arising as certain regions of the domain become inaccessible due to progressively blocked pathways.</figcaption></figure>

output with experimental measurements of the voltage evolution under constant current operation. Despite the fact that this criterion conflates kinetic, ohmic and transport losses, it remains the most common method for tracking device aging in literature, as it is straightforward to implement experimentally and provides an overall measure of the level of degradation of the device.

A key challenge in validating degradation models for PEMWEs is the scarcity of literature datasets that report, in a consistent manner, the combination of long-term degradation curves, detailed operating conditions, and PTL structural characteristics. Many studies present only polarization curves or voltage-time trajectories without sufficient information on the cell architecture, making direct quantitative comparison difficult. To address this limitation, we use the study of Suermann et al. (Suermann, 2019), which uniquely provides both

overpotential-time degradation data, with a detailed description of the PTL properties and operating conditions, as well as other electrochemical characterization techniques. This comprehensive data makes it a particularly suitable reference for benchmarking our model predictions against experimentally observed aging behavior. Indeed, the degradation rate is defined as the rate of voltage variation during aging experiments which is given in µV/h.

In their work, Suermann et al. reported using a normalized recirculation flow of 20 mL/(min.cm²), at ambient pressure and a temperature of 60°C. The corresponding operating conditions and model parameters adopted in the present study are summarized in Table 1. In this case, the PTL consists of 1 mm thick titanium with a fiber diameter of 20 µm at a porosity of 56%. Fig. 8b shows a direct comparison between the experimentally measured IR-Free degradation curve and the model prediction given identical conditions. The quasi-linear behavior of the degradation rate shown in Fig. 7b stems from the constant effective corrosion rate imposed by the geometry-evolution module. As the PTL porosity decreases linearly over time, the associated mass transport overpotential increases at a nearly constant rate under galvanostatic conditions. Minor deviations from perfect linearity result from pore isolation events and the local constriction of transport pathways, which are visible in Fig. 8. The magnitude (about 10 µV/h) is in good agreement with the IR-free values reported by Suermann et al. (Suermann, 2019) given a corresponding corrosion rate of 0.085 µm/day, which reinforces the idea that the imposed corrosion velocity corresponds to a physically reasonable accelerated-aging scenario. Fig. 9 illustrates the spatial distribution of the oxygen volume fraction within the PTL and flow channels at the start (left panel) and end (right panel) of the aging test. Initially, oxygen is relatively evenly distributed throughout the PTL and flow channels. By the end of the test, regions with high oxygen volume fraction (>0.9) become more prevalent, particularly near blocked pores and constricted pathways, underscoring a reduction in

<caption>**Table 1**</caption>

Summary of operating conditions, physical properties, and electrochemical parameters used in the simulations.

<table><tr><td>Temperature</td><td>60°C</td></tr><tr><td>Flowrate</td><td>20 mL/(min.cm²),</td></tr><tr><td>Bubble diameter</td><td>10 µm</td></tr><tr><td>Continuous phase density</td><td>1000 kg/m³</td></tr><tr><td>Disperse phase density</td><td>1.17 kg/m³</td></tr><tr><td>Exchange current density</td><td>1.9 × 10⁻⁸ A/cm²</td></tr><tr><td>Thermodynamic cell voltage</td><td>1.1948 V</td></tr><tr><td>Corrosion rate</td><td>0.085 µm/day</td></tr></table>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** (a) Effect of the exchange current density on the overpotential, as a function of time, and given a linear degradation of the catalyst performance, reaching 100% at 1000 h; (b) Comparison of the experimentally measured IR-Free degradation rate from Suermann et al. (Suermann, 2019) at end-of-test conditions, and the model prediction given the same current density. The degradation rate is defined as the rate of change of the cell voltage in µV/h.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Comparison of the oxygen volume fraction distribution within the PTL and flow channels at the start (left) and end (right) of the aging test. At the end of the test, regions with very high oxygen volume fraction (>0.9) become more prevalent, reflecting reduced transport efficiency and a more constricted porous network due to geometric evolution.</figcaption></figure>

mass-transport efficiency driven by corrosion induced geometric evolution. This visualization provides mechanistic insight into how local structural changes within the PTL progressively translate to the overall increase in cell overpotential.

## 8. Limitations of the model

The present model is intentionally reduced-order. First, corrosion is represented by a constant effective corrosion velocity instead of a mechanistic kinetic law. Second, the PTL is represented by a two-dimensional (2D) cross-sectional geometry that acts as a representative elementary volume. This geometry cannot reproduce full three-dimensional (3D) pore connectivity, tortuosity, or mechanical support. Third, the model is isothermal and does not incorporate membrane aging, electronic transport in the solid, or thermo-mechanical stresses, all of which play important roles in long-term degradation. It also does not yet incorporate electronic conduction through the PTL or bipolar plates through a charge transport or current conservation model, nor does it capture the impact of dynamic operating conditions such as load cycling. Despite these simplifications, the close agreement with both polarization data and long-term degradation curves highlights the robustness of the approach.

Several avenues exist for further refinement of the current model. One avenue is incorporating non-uniform and time-varying corrosion kinetics to capture more localized degradation phenomena, potentially coupled to the potential distribution over the entire PTL domain. Another is to extend to three-dimensional pore structures reconstructed from micro-computed tomography imaging of real PTLs instead of synthetically constructed porous geometries. Import Alternatively, more generalized pore structure geometries could be generated without relying on randomly positioned disks, taking full advantage of the level set method's capabilities to better represent the PTL and to allow explicit control over permeability and tortuosity. While the size of individual PTL particles can in principle be prescribed directly as a function of aging time, such an approach remains valid only in highly idealized configurations involving non-overlapping, perfectly spherical features. The use of the Level Set method therefore provides a more general and rigorous framework for geometry evolution, as it naturally handles particle overlap, non-spherical morphologies, and complex porous microstructures. In addition, it offers the flexibility to incorporate spatially and temporally varying corrosion rates in future developments.

It should also be noted that coupling electronic transport would enable direct prediction of full cell voltages without the need for IR correction, integrating the model into accelerated stress testing protocols and design studies for corrosion-resistant PTL materials and coatings. However, it is not straightforward to implement in 2D geometries as the porous domains are not connected, which makes this approach more viable for fully resolved 3D models. Furthermore, the Schiller-Naumann drag model should be interpreted as a closure approximation consistent with the dispersed-phase mixture formulation rather than as a detailed representation of individual bubble dynamics. According to the continuum-scale assumptions of the model, the influence of solid particles, pore confinement, and bubble deformation is intentionally neglected. Consequently, the current framework replaces explicit capillary and pore-scale effects with a reduced-order formulation that can capture essential transport trends across the entire PEM electrolyzer geometry.

By combining physics-based degradation mechanisms with numerical efficiency, the presented model provides a valuable tool for predicting device lifetime, guiding materials development, and ultimately supporting the durability-driven scale-up of PEMWE technology. However, the modular structure of this framework makes it directly compatible with more advanced sub-models, such as those involving mechanistic corrosion kinetics, membrane degradation, heat transport, and 3D µCT-resolved PTL geometries.

This work focuses on performance degradation; however, PEMWE degradation is linked to safety–critical phenomena, including hydrogen crossover, membrane thinning, pinhole formation, and local hot spots caused by increased interfacial resistance. Moreover, it is well known that hydrogen crossover approaching 4 vol% in oxygen reaches flamability limits. Additionally, uneven PTL degradation can modify local pressure fields and accelerate membrane mechanical failure. While the present model does not compute crossover or temperature fields, the predicted evolution of PTL morphology and oxygen saturation can be fed directly into established crossover and thermo-mechanical models.

## 9. Conclusions

In this work, we developed and implemented a two-component multiphase numerical framework to investigate the corrosion-driven degradation of proton exchange membrane water electrolyzers. The approach couples a geometry tracking module, based on the Level Set method, with a performance evaluation module, based on a mixture model for oxygen-water two-phase flow in the porous transport layer. This combination enabled us to dynamically link morphological changes in the PTL to transport and kinetic losses at the catalyst layer, thereby providing a mechanistic description of how corrosion progressively reduces cell performance over time.

The main outcomes of this study can be summarized as follows:

*   Validation in the pristine regime: The model successfully reproduced experimental IR-free polarization curves reported in the literature. After accounting for an effective internal resistance, the predicted polarization curves showed close agreement with the full experimental data, supporting the physical fidelity of the transport and reaction formulations in the absence of degradation.

*   **Geometry evolution under corrosion:** The Level Set method is effective in tracking porosity loss and topological changes in the PTL due to deposition of corrosion products. The results indicate that, under a constant corrosion velocity, the effective porosity decreases almost linearly, but with localized deviations arising from the isolation of pore regions as transport pathways become blocked. While in real electrolyzers it is difficult to assign and quantify the different contributions of degradation to different mechanisms, the reduction in transport that is implemented in this model effectively encapsulates the overall degradation and voltage increase. Future integration with extensive electrochemical characterizations techniques could enable the precise quantification of each mechanism's contribution.

*   **Coupled impact of catalyst degradation and bubble-induced transport losses:** The model demonstrated that catalyst degradation, represented as a progressive reduction of exchange current density at the catalyst layer, mainly affects the late stages of operation, where performance collapses rapidly once more than ~ 90% of the catalyst
activity is lost. By contrast, transport losses due to oxygen bubble accumulation and pore blockage grow more gradually and contribute significantly throughout the device lifetime. This supports experimental observations that mass transport and ohmic losses dominate at high current densities. It also shows that catalyst degradation in electrolyzers is not as limiting as transport and resistive losses.

Reproduction of long-term degradation curves: Beyond polarization validation, the model was applied to simulate constant-current operation over extended periods. The predicted voltage rise as a function of time matched well with degradation curves reported in experimental studies by tuning only physically meaningful parameters such as corrosion rate and catalyst decay constant. This demonstrates that the framework can serve as an effective tool for interpreting and deconvoluting experimental aging data.

Overall, this work establishes a reduced-order multiphysics framework that propagates prescribed morphology evolution scenarios to device-scale voltage degradation. Although it does not mechanistically predict corrosion kinetics, membrane failure, or thermo-mechanical stresses, the framework offers a physically grounded approach to isolating and quantifying the sensitivity of PEMWE performance to evolving PTL geometry. This lays the groundwork for coupling more detailed mechanistic or safety-oriented models.

## CRediT authorship contribution statement

**Youssef Kharchouf:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Data curation, Conceptualization. **Loïc Assaud:** Writing – review & editing, Validation, Supervision, Resources, Project administration, Funding acquisition, Data curation, Conceptualization. **Mireille Turmine:** Writing – review & editing, Visualization, Validation, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Vincent Vivier:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

The authors gratefully acknowledge the financial support of the French National Research Agency (agreement ANR-22-PEHY-0009 under the France 2030 program and agreement ANR-22-CE93-0008-01).

## Data availability

Data will be made available on request.

## References

Abdin, Z., et al., 2020. Hydrogen as an energy vector. Renew. Sustain. Energy Rev. 120, 109620. https://doi.org/10.1016/j.rser.2019.109620

Ji, M., Wang, J., 2021. Review and comparison of various hydrogen production methods based on costs and life cycle impact assessment indicators. Int. J. Hydrogen Energy 46, 38612–38635. https://doi.org/10.1016/j.ijhydene.2021.09.142

Ursua, A., et al., 2012. Hydrogen production from water electrolysis: current status and future trends. Proc. IEEE 100, 410–426. https://doi.org/10.1109/JPROC.2011.2156750

Wang, T., et al., 2022. PEM water electrolysis for hydrogen production: fundamentals, advances, and prospects. Carbon Neutrality 1, 21. https://doi.org/10.1007/s43979-022-00022-8

El-Shafie, M., 2023. Hydrogen production by water electrolysis technologies: a review. Results Eng. 20, 101426. https://doi.org/10.1016/j.rineng.2023.101426

Zhang, K., et al., 2022. Status and perspectives of key materials for PEM electrolyzer. Nano Res. Energy 1, e9120032. https://doi.org/10.26599/NRE.2022.9120032

Norazahar, N., et al., 2024. Degradation modelling and reliability analysis of PEM electrolyzer. Int. J. Hydrogen Energy 50, 842–856. https://doi.org/10.1016/j.ijhydene.2023.07.153

Prestat, M., 2023. Corrosion of structural components of proton exchange membrane water electrolyzer anodes: a review. J. Power Sources 556, 232469. https://doi.org/10.1016/j.jpowsour.2022.232469

Grigoriev, S.A., et al., 2014. Failure of PEM water electrolysis cells: case study involving anode dissolution and membrane thinning. Int. J. Hydrogen Energy 39, 20440–20446. https://doi.org/10.1016/j.ijhydene.2014.05.043

Omrani, R., Shabani, B., 2021. Hydrogen crossover in proton exchange membrane electrolysers. Electrochim. Acta 377, 138085. https://doi.org/10.1016/j.electacta.2021.138085

Siracusano, S., et al., 2018. Degradation issues of PEM electrolysis MEAs. Renew. Energy 123, 52–57. https://doi.org/10.1016/j.renene.2018.02.024

Papakonstantinou, G., et al., 2020. Degradation study of a PEM water electrolyzer under dynamic operation. Appl. Energy 280, 115911. https://doi.org/10.1016/j.apenergy.2020.115911

Liu, C., et al., 2018. Performance enhancement of PEM electrolyzers through iridium-coated titanium porous transport layers. Electrochem. Commun. 97, 96–99. https://doi.org/10.1016/j.elecom.2018.10.021

Yuan, X.-Z., et al., 2022. The porous transport layer in PEM water electrolysis, Sustainable. Energy Fuels 6, 1824–1853. https://doi.org/10.1039/D2SE00260D

Suermann, M., et al., 2019. Degradation of PEM water electrolysis cells: looking beyond the cell voltage increase. J. Electrochem. Soc. 166, F645. https://doi.org/10.1149/2.1451910jes

Babic, U., et al., 2017. Critical review—identifying critical gaps for polymer electrolyte water electrolysis development. J. Electrochem. Soc. 164, F387–F399. https://doi.org/10.1149/2.1441704jes

Arbabi, F., et al., 2016. Three-dimensional CFD modelling of oxygen bubble transport in porous transport layers. J. Electrochem. Soc. 163, F3062–F3069. https://doi.org/10.1149/2.0091611jes

Paliwal, S., et al., 2021. Lattice Boltzmann study of water–oxygen distributions in porous transport layers. Int. J. Hydrogen Energy 46, 22747–22762. https://doi.org/10.1016/j.ijhydene.2021.04.112

Tijani, A.S., Rahim, A.H.A., 2016. Numerical modeling of the effect of operating variables on Faraday efficiency in PEM electrolyzer. Procedia Technol. 26, 419–427. https://doi.org/10.1016/j.protcy.2016.08.054

Falcão, D.S., Pinto, A.M.F.R., 2020. A review on PEM electrolyzer modelling: guidelines for beginners. J. Clean. Prod. 261, 121184. https://doi.org/10.1016/j.jclepro.2020.121184

Arbabi, F., et al., 2014. Feasibility study of microfluidic platforms for visualizing bubble flows. J. Power Sources 258, 142–149. https://doi.org/10.1016/j.jpowsour.2014.02.042

Zeng, Q., et al., 2024. Numerical simulation of comprehensive performance of PEM electrolytic cell with microporous layer. Int. J. Hydrogen Energy 77, 495–510. https://doi.org/10.1016/j.ijhydene.2024.06.219

Harrison, K.W., et al., 2005. Semiempirical model for determining PEM electrolyzer stack characteristics. J. Fuel Cell Sci. Technol. 3, 220–223. https://doi.org/10.1115/1.2174072

Ruuskanen, V., et al., 2017. PEM water electrolyzer model for power-hardware-in-loop simulator. Int. J. Hydrogen Energy 42, 10775–10784. https://doi.org/10.1016/j.ijhydene.2017.03.046

Chandesris, M., et al., 2015. Membrane degradation in PEM water electrolyzer. Int. J. Hydrogen Energy 40, 1353–1366. https://doi.org/10.1016/j.ijhydene.2014.11.111

Legala, A., et al., 2022. Machine learning modelling for PEM fuel cell performance. Energy AI 10, 100183. https://doi.org/10.1016/j.egyai.2022.100183

Madhavan, P.V., et al., 2025. Data-driven modelling of corrosion behaviour in coated porous transport layers. Artif. Intell. Chem. 3, 100086. https://doi.org/10.1016/j.aichem.2025.100086

Sethian, J.A., 1996. Level set methods: evolving interfaces in geometry, fluid mechanics, computer vision, and materials science. Cambridge Univ. Press.

Guan, X., et al., 2020. CFD simulation of gas–liquid flow in stirred tanks. Chem. Eng. J. 386, 121554. https://doi.org/10.1016/j.cej.2019.04.134

Karunarathne, S.S., Tokheim, L.-A., 2017. Comparison of the influence of drag models in CFD simulation. Proc. SIMS 151–156. https://doi.org/10.3384/ecp17138151

P. Satjaritanun et al., Observation of preferential pathways for oxygen removal through porous transport layers, *iScience* 23 (2020) 101783. https://doi.org/10.1016/j.isci.2020.101783

Altaf, H., et al., 2020. Steady-state water drainage by oxygen in anodic porous transport layers. Processes 8, 362. https://doi.org/10.3390/pr8030362

Chen, Q., et al., 2020. Two-dimensional multiphysics modeling of porous transport layer in PEM electrolyzer. Int. J. Hydrogen Energy 45, 32984–32994. https://doi.org/10.1016/j.ijhydene.2020.09.148

Sinapan, I., et al., 2024. Revealing anodic multi-class bubble dynamics using deep learning. Fuel 364, 131112. https://doi.org/10.1016/j.fuel.2024.131112

Suermann, M., et al., 2015. Investigation of mass transport losses in polymer electrolyte electrolysis cells. ECS Trans. 69, 1141–1148. https://doi.org/10.1149/06917.1141ecst