# Understanding the Operation of a Gas Diffusion Electrode Setup for the Oxygen Reduction Reaction: Experiment versus 3D Multiphysics Modeling

Erwan Tardy, Raphaël Riasse, Florent Vandenberghe, Florence Druart, Marian Chatenet, Antoine Bonnefont

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

## To cite this version:

Erwan Tardy, Raphaël Riasse, Florent Vandenberghe, Florence Druart, Marian Chatenet, et al. Understanding the Operation of a Gas Diffusion Electrode Setup for the Oxygen Reduction Reaction: Experiment versus 3D Multiphysics Modeling. ChemElectroChem, 2025, 12 (17), 10.1002/celc.202500172. [https://hal.science/hal-05190352v1](https://hal.science/hal-05190352v1)

HAL Id: hal-05190352

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

Submitted on 28 Jul 2025

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_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Understanding the Operation of a Gas Diffusion Electrode Setup for the Oxygen Reduction Reaction: Experiment versus 3D Multiphysics Modeling

Erwan Tardy, Raphaël Riasse, Florent Vandenberghe, Florence Druart, Marian Chatenet,* and Antoine Bonnefont*

Proton exchange membrane fuel cells (PEMFC) require highly efficient oxygen reduction reaction (ORR) electrocatalysts. The intrinsic ORR performance of advanced ORR catalysts (measured in rotating disk electrode, RDE) is often not obtained in membrane electrode assembly (MEA), which denotes for RDE inability to forecast intrinsic activity at large current density/low potential, owing to severe mass-transport limitation. The gas diffusion electrode (GDE) is a relevant tool to assess the intrinsic ORR activity of catalysts at high current density/low potential, so it enables to better forecast their performance in PEMFC MEA. Herein, ORR kinetics is studied in a GDE via cyclic voltammetry and electrochemical

impedance spectroscopy; the polarization curve is modeled using multiphysics and multicomponent 3D simulations. The model allows to investigate the interplay between the electrochemical kinetics and the mass-transport of reactant and products in the flow channels of the monopolar plate, gas diffusion layer, porous electrode, and solution. The simulations highlight the significant impact of partial water flooding in the catalyst layer—even at a minimal thickness of 0.6 µm—on the shape of the GDE polarization curves and suggest that mass-transport limitation inside the catalyst layer may be limiting in a GDE setup, specifically at high current density.

## 1. Introduction

Although proton exchange membrane fuel cells (PEMFC) are already commercialized in various systems, the development of highly efficient oxygen reduction reaction (ORR) electrocatalysts still remains crucial to enhance their performance, reduce capital cost of the systems, and address the availability of the platinum group metals that are used in their electrodes.[1] While the scientific community has been extremely successful in developing new and sometimes more active electrocatalysts for the ORR, it still struggles to transfer their intrinsic performance in operating PEMFC membrane electrode assemblies (MEAs): very often, their practical performance remains far from those expected from model experiments, that are still to date usually performed using the rotating disk electrode setup, RDE (the setup of choice to assess a material's ORR performance[2-7] in most labs on the planet[8]).

There are multiple reasons to account for the RDE inability to predict MEA performance. First, and most importantly, in RDE configuration, the $O_2$ reactant must be dissolved in the electrolyte, where it has much smaller diffusivity ($10^{-5}$ cm$^2$ s$^{-1}$ range) and concentration/activity ($10^{-6}$ mol cm$^{-3}$ range for the solubility) than in the gas phase, which implies much slower mass-transfer rate and related smaller limiting current density than in real PEMFC conditions, typically by ≈3 orders of magnitude. As a result of this and second, the ORR activity data can only be extracted at high potential values in RDE ($E > 0.85$ V vs reversible hydrogen electrode, RHE), although PEMFC optimal performance is reached at $0.6 < E < 0.8$ V versus RHE, the cell voltage essentially resulting from the cathode potential. Third, in a PEMFC MEA, proton conduction in the thin layer of ionomers covering the electrocatalyst nanoparticles can be limiting, which can affect their operation at the high current densities at stake in MEA configuration.[9,10] Because of this, one cannot be sure that the electrocatalysts are properly implemented in MEA, where emphasized mass-transport of reactants, $H^+$ and $O_2$, to the catalytic sites might not be granted for present ionomers and MEA structure.[11,12] To go beyond these inherent difficulties of the RDE, recent studies focused on advanced electrochemical characterization setups like the floating electrode (FE),[13] the gas-diffusion electrode (GDE),[14,15] or differential cell (DC),[16,17] that enable to reach much higher ORR current densities at the lab scale and therefore better match MEA data, a prerequisite for relevant benchmarking. As a matter of fact, these setups enable reliable comparison of ORR predictions from model measurements with unit PEMFC on a broad voltage/current range, which is only possible on a narrow range with the RDE.[18]

E. Tardy, R. Riasse, F. Vandenberghe, F. Druart, M. Chatenet, A. Bonnefont
CNRS, Grenoble INP, LEPMI
Université Grenoble Alpes, Université Savoie Mont Blanc
F-38000 Grenoble, France
E-mail: marian.chatenet@grenoble-inp.fr
pierre-antoine.bonnefont@grenoble-inp.fr

R. Riasse
CEA, Liten
Université Grenoble Alpes
38000 Grenoble, France

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

Supporting information for this article is available on the WWW under https://doi.org/10.1002/celc.202500172

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

© 2025 The Author(s). ChemElectroChem published by Wiley-VCH GmbH. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
<figure><img src="image_6.png" /></figure>

The common goal of these techniques is to assess the intrinsic kinetic activity of the studied ORR catalyst, which requires that the only kinetics limitation originates from the charge transfer reaction (the ORR) at the Pt catalyst | electrolyte/ionomer interface. This naturally requires that the high-frequency resistance is measured and corrected, and that any kind of mass-transfer process is not limiting. In the particular case of the GDE, even using low Pt loading (so to be catalyst-limited) to evaluate the intrinsic kinetics parameters (the strategy employed in[16]), the ORR curves might not be totally free of any mass-transfer limitation. To be more specific, one can a priori not exclude that proton mass-transport in the electrolyte solution, electronic and ionic conductivity within the catalyst layer thickness or water flooding in the catalyst layer or in the gas diffusion layer (GDL) do affect the ORR data, which implicitly means that the operating conditions e.g., gas humidification and electrolyte temperature do probably influence the ORR data, which of course renders their interpretation complex.

In this article, we combined galvanostatic electrochemical impedance spectroscopy (GEIS) and 3D-GDE modeling to understand ORR data on GDE for a bimetallic state-of-the-art catalyst Pt₃Co/C. It shows that, even for low-loaded catalyst layers in differential conditions with high flux of reactant, mass-transfer to/within the catalyst layer cannot be ignored, stressing that the device, despite its huge fundamental interest, must be used wisely to predict appropriately the intrinsic catalytic activity of a given catalyst material.

## 2. Experimental Section

### 2.1. GDE Cell, Electrode, Electrolyte, and Gases

The full GDE setup, displayed in Figure 1a, has been first detailed and published in ref. [16]. Inspired by the design of Wiberg et al.,[19] it was further developed in collaboration with the CEA LITEN, so that the working electrode current collector is removable and can be cleaned. It is a graphite monopolar plate (MP) that separates the stainless-steel lower body of the cell from

the upper body: a polytetrafluoroethylene (PTFE) beaker filled with the electrolyte. The porous working electrode (3 mm-diameter – 0.071 cm² geometric surface area), a catalyst layer (CL) deposited onto a hydrophobic GDL with microporous layer (Freudenberg H23C8), is placed onto the graphite monopolar plate at the center of a “donut” of the same catalyst-free GDL, as depicted in Figure 1b. This donut has external and internal diameters of 6 and 3 mm respectively. The GDL-supported active layer (with no membrane), directly separates the gas flow field from the liquid electrolyte contained in the PTFE beaker, the interface hence being a gas | liquid one, schematized in Figure 1c. This is a strongly different situation than the one encountered in typical fully-flooded thin-film RDE interface and is closer to real single-cell gas | hydrated-polymer interface. To ensure a good control of the electrochemical parameter, a three-electrode configuration is used, with the counter (platinum mesh) and reference (RHE, Gaskatel) electrodes located in the liquid upper part, the reactant being fed to the GDE from underneath via the gas channel of the graphite monopolar plate. The latter bears 11 parallel channels of ≈250 µm-width (same for ribs) and 400 µm-depth, which bring the gas (oxygen, carbon monoxide, or argon) to the active layer. This leads to a strongly enhanced mass-transport rate of the gaseous reactant to the catalyst sites, the great advantage of this setup compared to the classical flooded thin-film RDE.

The electrolyte preparation is made with extremely pure products, here 70% perchloric acid (HClO₄, Rotipuran Ultra, Roth) and ultrapure water (18.2 MΩ cm, total organic carbon (TOC) < 3 ppb, Elix + Milli-Q, Millipore, Merck); a concentration of 1 mol L⁻¹ of HClO₄ is used in these experiments to ensure that enough protons are available at the interface during the ORR measurements and to prevent any limitation by the transport of protons.

Pure gases (from Messer) are used for the electrochemical measurements: Ar (99.999%), CO (99.997%), O₂ (99.999%). The mass flow controllers offer regulated and reproducible gas flows for most of them, CO being added to the Ar flow. All gases are humidified with a bubbler at room temperature prior to entering the cell.

<figure><img src="image_7.png" /><figcaption>**Figure 1.** Schematics of a) the gas diffusion electrode setup, b) the donut-like working electrode (the center part being the catalyzed (active) electrode, the external ring being the electron-conductive electrode “holder”), and of c) the GDE gas | liquid interface.</figcaption></figure>
## 2.2. Catalyst Materials

The active part of the GDE (catalyst layer) was prepared from platinum-cobalt alloyed nanoparticles supported on carbon, purchased from Tanaka Kikinzoku Kogyo: Pt₃Co nanoparticles on Vulcan (TEC36V52, referred as PtCo/VC). The received powder is used without any additional treatment for the realization of the catalytic inks. The preparation of the active layers, the chemical analyses to qualify the catalyst loading and the cleaning of the experimental setup are more detailed in ref. [16]. A very low catalyst loading of 20 µgPt cm⁻² was used, so to be catalyst-limited in the GDE experiments, as detailed in [16] this corresponds to a very thin CL (see hereafter and Figure 2).

## 2.3. Electrochemical Characterizations

For the electrochemical measurements with this three-electrode setup, the protocol used is the same as the one realized in ref. [16]. The cell main body is thermally regulated at 25 °C using a thermostatic bath (F12, Julabo), and the electrodes are controlled with a Bio-Logic SP-150 potentiostat equipped with a frequency analyzer. At the beginning of each (except impedance) electrochemical technique, the high-frequency resistance is determined at the open-circuit potential at a frequency of 100 kHz with a rather small amplitude (10 mV). This resistance is then dynamically corrected at 85% and the 15% remaining are post-corrected to access ohmic-drop-corrected (IR-free) experimental curves.

After a well-known break-in protocol consisting in potential cycling that has been used in many studies as it is fast and enables to reach the highest specific activity possible for a given catalyst,[20,21] a sequence of cyclic voltammetry (CV, 3 cycles) and a carbon monoxide (CO) stripping protocol are performed to measure the electrochemical active surface area (ECSA) used in the following to determine the active specific surface area $a_v$. An example of CO stripping and its subsequent cycle used to correct

the baseline is shown in Figure S1 (Supporting Information). Once again, additional information on how to perform these techniques is given in details in ref. [16].

Polarization curves in GDE are done using slow scan rate CVs (5 mV s⁻¹). The potential of the working electrode is held for one minute at 0.2 V versus RHE (the lower potential boundary) and cycled with an upper potential boundary of 1.05 V versus RHE, in a high gas flow of oxygen (60 mL min⁻¹). Seven cycles are nevertheless necessary to reach a steady state: the performance keeps on increasing as the cycles go. The experimental curves displayed are the increasing potential scan of the 7th cycle realized for each catalyst.

After these polarization curves, the maximum current density reached is set as a limit to not exceed for GEIS. A succession of impedance spectra is done at various currents: 0, 1, 10, 25, and 50 mA, which correspond to the following current densities: 0, 14, 140, 350, and 700 mA cm⁻². The frequency varies from 100 kHz to 1 Hz with 8 points per decade and an amplitude equal to 10% of the steady state point ($I=0.1$ mA being the amplitude for the impedances made at 0 mA).

The GEIS spectra are fitted with a transmission line equation considering the ORR and the proton conduction within the catalyst layer, while disregarding gas transport in the channel, GDL, and catalyst layer, as well as electron conduction. For this simple case, the impedance can be expressed by:[22,23]

$$Z_p = \sqrt{\frac{Z_{\text{loc}}}{a_v \sigma_i^{\text{eff}}}} \coth \left( \sqrt{\frac{a_v}{\sigma_i^{\text{eff}} Z_{\text{loc}}}} \delta_L \right) \quad (1)$$

with

$$Z_{\text{loc}} = \frac{R_{\text{ct}}}{j\omega C_{\text{dl}} R_{\text{ct}} + 1} \quad (2)$$

where $a_v$ is the ratio of the platinum surface area to the volume of the porous catalyst layer, $\delta_{CL}$, the thickness of the catalyst layer, $\omega$,

<figure><img src="image_8.png" /><figcaption>Figure 2. a) Modeled geometry of the MP through which the 3D gas diffusion electrode is fed by reactants, b) 3D model with domain discretization of the liquid electrolyte and catalyst layer (GDE), and c) scanning electron microscopy micrograph of the cross-section of the GDL-supported active layer; the average active layer thickness is around 0.6 µm for the considered catalyst loading (21 µgPt cm⁻²).</figcaption></figure>
<figure><img src="image_9.png" /></figure>

<table><caption>**Table 1.** Model configuration.</caption><caption>Dimensions of the modeling domains from experimental setup.</caption><thead><tr><th>Symbol</th><th>Parameter</th><th>Values</th><th>Unit</th></tr></thead><tbody><tr><td><strong>Greek</strong></td><td></td><td></td><td></td></tr><tr><td>δ<sub>CL</sub></td><td>CL thickness</td><td>0.6</td><td>μm</td></tr><tr><td>δ<sub>GDL</sub></td><td>GDL thickness<sup>a)</sup></td><td>200</td><td>μm</td></tr><tr><td>δ<sub>le</sub></td><td>liquid electrolyte thickness</td><td>Variable</td><td>μm</td></tr><tr><td><strong>Latin</strong></td><td></td><td></td><td></td></tr><tr><td>h<sub>ch</sub></td><td>channel height</td><td>0.4</td><td>mm</td></tr><tr><td>l<sub>cell</sub></td><td>cell length</td><td>5.5</td><td>mm</td></tr><tr><td>l<sub>ch</sub></td><td>channel length</td><td>5.5</td><td>mm</td></tr><tr><td>S<sup>geo</sup><sub>CL</sub></td><td>CL geometric surface area<sup>b),c)</sup></td><td>0.071</td><td>cm<sup>2</sup></td></tr><tr><td>w<sub>cell</sub></td><td>cell width</td><td>5.25</td><td>mm</td></tr><tr><td>w<sub>ch</sub></td><td>channel width</td><td>0.25</td><td>mm</td></tr><tr><td>w<sub>fib</sub></td><td>rib width</td><td>0.25</td><td>mm</td></tr></tbody></table>

<table><caption>**Operating and boundary conditions.**</caption><thead><tr><th>Symbol</th><th>Parameter</th><th>Values</th><th>Unit</th></tr></thead><tbody><tr><td><strong>Greek</strong></td><td></td><td></td><td></td></tr><tr><td>Φ<sup>GDL,MP</sup><sub>e</sub></td><td>electronic potential applied at the GDL | MP interface (solid red line in Figure 4)</td><td>1–0.3</td><td>V</td></tr><tr><td>Φ<sup>CL,electrolyte</sup><sub>l</sub></td><td>ionic potential applied at the CL | liquid electrolyte interface (solid black line in Figure 4)</td><td>0</td><td>V</td></tr><tr><td><strong>Latin</strong></td><td></td><td></td><td></td></tr><tr><td>p<sup>in</sup><sub>g</sub></td><td>inlet gas pressure</td><td>1</td><td>bar</td></tr><tr><td>RH</td><td>relative humidity</td><td>100%</td><td>–</td></tr><tr><td>T</td><td>temperature</td><td>25</td><td>°C</td></tr><tr><td>x<sup>in</sup><sub>O<sub>2</sub></sub></td><td>inlet molar fraction of the oxygen</td><td>1 – x<sup>in</sup><sub>vap</sub></td><td>–</td></tr><tr><td>x<sup>in</sup><sub>vap</sub></td><td>inlet molar fraction of the water vapor</td><td>p<sup>in</sup><sub>g</sub>(T)/p<sup>in</sup><sub>g</sub>RH</td><td>–</td></tr><tr><td>–</td><td>inlet gas</td><td>pure O<sub>2</sub></td><td>–</td></tr><tr><td>–</td><td>inlet volume flow rate</td><td>60</td><td>mL min<sup>−1</sup></td></tr></tbody></table>

<table><caption>**Material and electrochemical properties used in the model.**</caption><thead><tr><th>Symbol</th><th>Parameter</th><th>Values</th><th>Unit</th></tr></thead><tbody><tr><td><strong>Greek</strong></td><td></td><td></td><td></td></tr><tr><td>α</td><td>charge transfer coefficient for reduction</td><td>0.5</td><td>–</td></tr><tr><td>β</td><td>charge transfer coefficient for oxidation</td><td>0.5</td><td>–</td></tr><tr><td>ε</td><td>Porosity</td><td>CL: 0.37<br>GC: 1<br>GDL: 0.8</td><td>–</td></tr><tr><td>ε<sub>e</sub></td><td>volume fraction of the electronically conductive phase</td><td>CL: 0.37<br>GDL: 0.2</td><td>–</td></tr><tr><td>ε<sub>l</sub></td><td>volume fraction of the ionically conductive phase</td><td>CL: 0.26</td><td>–</td></tr><tr><td>K</td><td>absolute permeability</td><td>CL: 8 × 10<sup>−12</sup><br>GDL: 6 × 10<sup>−8</sup></td><td>cm<sup>2</sup></td></tr></tbody></table>

<table><caption>**Table 1.** Continued.</caption><caption>Material and electrochemical properties used in the model.</caption><thead><tr><th>Symbol</th><th>Parameter</th><th>Values</th><th>Unit</th></tr></thead><tbody><tr><td>μ<sub>g</sub></td><td>gas viscosity</td><td>O<sub>2</sub>: 2.026 × 10<sup>−5</sup><br>H<sub>2</sub>: 8.789 × 10<sup>−6</sup></td><td>Pa s</td></tr><tr><td>σ<sup>bulk</sup><sub>e</sub></td><td>bulk electronic conductivity<sup>26)</sup></td><td>CL: 120<br>GDL: 120</td><td>S cm<sup>−1</sup></td></tr><tr><td>σ<sup>eff</sup><sub>l</sub></td><td>effective ionic conductivity<sup>d)</sup></td><td>CL: 0.0538</td><td>S m<sup>−1</sup></td></tr><tr><td><strong>Latin</strong></td><td></td><td></td><td></td></tr><tr><td>a<sub>v</sub></td><td>active specific surface area</td><td>l<sub>pt</sub>=ECSA / δ<sub>CL</sub></td><td>m<sup>2</sup> m<sup>3</sup> geo</td></tr><tr><td>c<sub>H<sub>2</sub>,ref</sub></td><td>reference proton concentration</td><td>1</td><td>M</td></tr><tr><td>D<sup>bulk</sup><sub>O<sub>2</sub>/vap</sub></td><td>bulk binary diffusion coefficient for O<sub>2</sub>/vap gaseous species</td><td>1 × 10<sup>−5</sup></td><td>m<sup>2</sup> s<sup>−1</sup></td></tr><tr><td>D<sup>bulk</sup><sub>O<sub>2</sub>/liquid</sub></td><td>bulk diffusion coefficient of gaseous oxygen in liquid water</td><td>1 × 10<sup>−9</sup></td><td>m<sup>2</sup> s<sup>−1</sup></td></tr><tr><td>ECSA</td><td>electrochemically active surface area<sup>e)</sup></td><td>30</td><td>m<sup>2</sup> g<sup>−1</sup> pt</td></tr><tr><td>E<sub>eq</sub></td><td>equilibrium potential for ORR</td><td>1.23</td><td>V</td></tr><tr><td>H<sub>O<sub>2</sub>/liquid</sub></td><td>Henry's constant for gaseous oxygen in liquid water at 25 °C<sup>27)</sup></td><td>1.3 × 10<sup>−3</sup></td><td>mol (L atm)<sup>−1</sup></td></tr><tr><td>i<sup>ORR</sup><sub>0</sub></td><td>exchange current density for ORR</td><td>3 × 10<sup>−6</sup></td><td>A cm<sup>−2</sup></td></tr><tr><td>L<sub>ptCo</sub></td><td>Pt<sub>3</sub>Co loading<sup>e)</sup></td><td>21.3</td><td>μg<sub>pt</sub> cm<sup>−2</sup> geo</td></tr><tr><td>M<sub>H<sub>2</sub></sub></td><td>molar mass of dihydrogen</td><td>2</td><td>g mol<sup>−1</sup></td></tr><tr><td>M<sub>H<sub>2</sub>O</sub></td><td>molar mass of water</td><td>18</td><td>g mol<sup>−1</sup></td></tr><tr><td>M<sub>O<sub>2</sub></sub></td><td>molar mass of dioxygen</td><td>32</td><td>g mol<sup>−1</sup></td></tr><tr><td>p<sup>in</sup><sub>O<sub>2</sub></sub></td><td>initial partial pressure of the oxygen</td><td>96 858</td><td>Pa</td></tr><tr><td>p<sub>ref</sub></td><td>reference pressure</td><td>1</td><td>atm</td></tr><tr><td>R</td><td>universal gas constant</td><td>8.3145</td><td>J (mol K)<sup>−1</sup></td></tr></tbody></table>

a) From technical sheet for a H23C8 Freudenberg GDL @ 1 MPa.

b) Corresponding to a 3 mm-diameter ink deposit. c) HCIO<sub>4</sub> solution. d) From GEIS fit (see *Electrochemical characterizations* section).

e) From experimental tests.

the angular frequency, j, the imaginary number, σ<sup>eff</sup><sub>j</sub>, the effective ionic conductivity across the thickness of the catalyst layer, C<sub>dl</sub>, the double layer capacitance, and R<sub>ct</sub>, the charge transfer resistance (Table 1). By fitting the experimental impedance spectra obtained at 1 mA, where the spectra are assumed to be unaffected by the oxygen and proton mass-transport limitations, the value of C<sub>dl</sub>, R<sub>ct</sub>, and σ<sup>eff</sup><sub>j</sub> can be determined (Table S1, Supporting Information).

## 2.4. Mathematical Model

The mathematical model is developed to represent a gas diffusion electrode as detailed in the previous section (Figure 1a). The model used to study oxygen reduction reaction is
<figure><img src="image_10.png" /></figure>

three-dimensional, multiphysics, and multicomponents. Figure 2 shows the schematic of the modeling domain, which reproduces the geometry of the GDE used to obtain the experimental data. It consists of a part of the setup with 11 straight gas channels (GCs), one gas diffusion layer (GDL), one CL, and one domain to represent the liquid electrolyte. All physics are discussed below, where mass, momentum, species, electronic, and ionic charge are conserved within the model. The main assumptions of the model are as follows: 1) the gas mixture is considered to be ideal; it flows in a laminar regime. 2) A full hydration of the catalyst layer is assumed (water activity equal to the unity), which is reasonable owing to the fact that fully-humidified gases were fed to the working electrode; the working electrode is in direct contact with the liquid aqueous electrolyte and water production proceeds during acidic ORR. 3) The porous media exhibit isotropic characteristics, and the transport phenomena within them are described using effective properties. 4) The pressure in the free fluid and the pressure in the pores is continuously over the interface between the free flow domain and the porous domain. 5) A single-phase flow is considered. 6) Since the experimental GDE is thermostated, the temperature is assumed to be constant in the simulations, corresponding to an isothermal model.

## 2.5. Gas Flow in Channels

Gas flow in channels is described by the conservation laws of mass and momentum:

$$ \nabla \cdot (\rho_g u_g) = 0 \qquad (3) $$

$$ \rho_g(u_g \cdot \nabla)u_g = \nabla \cdot \left[ -p_g l + \mu_g (\nabla u_g + (\nabla u_g)^T) - \frac{2}{3} \mu_g (\nabla \cdot u_g) l \right] \quad (4) $$

where $\rho_g$, $u_g$, $p_g$, and $\mu_g$ are the density, velocity, pressure, and viscosity of the gas mixture, respectively, while $l$ represents the identity matrix. The gas mixture density is computed assuming ideal gases:

$$ \rho_g = \frac{p_g}{RT} \left( \sum_{i=1}^{n} x_i M_i \right) \qquad (5) $$

in which $R$ and $T$ are the universal gas constant and temperature, respectively. In addition, $x_i$ and $M_i$ represent the mole fraction and molar mass of species $i$ ($i$: O$_2$, H$_2$O). The gas mixture viscosity is calculated by the Wilke formula:$^{[24]}$

$$ \mu_g = \sum_{i=1}^{n} \frac{x_i \mu_i}{\sum_{j=1}^{n} x_j \psi_{ij}} \qquad (6) $$

$$ \psi_{ij} = \frac{\left[ 1 + \left( \frac{M_i}{M_j} \right)^{\frac{1}{n}} \left( \frac{\mu_i}{\mu_j} \right)^{\frac{1}{n}} \right]^2}{\sqrt{8 \left( 1 + \frac{M_i}{M_j} \right)}} \qquad (7) $$

where $\mu_i$ and $\psi_{ij}$ are the viscosity of species $i$ and mixing function, respectively.

## 2.6. Gas Flow in Porous Media

Gas phase in porous media is subject to a mass-conservation equation (Equation (8)) and a momentum-conservation equation described by Brinkman's law (Equation (9)):

$$ \nabla \cdot (\rho_g u_g) = R_g \qquad (8) $$

$$ \frac{1}{\epsilon^2} \rho_g (u_g \cdot \nabla) u_g = \nabla \cdot \left[ -p_g l + \mu_g \frac{1}{\epsilon} (\nabla u_g + (\nabla u_g)^T) \right. \\ \left. - \frac{2}{3} \mu_g \frac{1}{\epsilon} (\nabla \cdot u_g) l \right] - \left( \mu_g \kappa^{-1} + \frac{R_g}{\epsilon^2} \right) u_g \qquad (9) $$

where $R_g$, $\epsilon$, and $\kappa$ are the mass source term, porosity, and absolute permeability, respectively. The mass source term represents the oxygen consumption and water production during ORR in the catalyst layer:

$$ R_g = R_{H_2O} + R_{O_2} \quad \text{in the CL} \qquad (10) $$

$$ R_{O_2} = -\frac{|i_v|}{4F} M_{O_2} \quad \text{in the CL} \qquad (11) $$

$$ R_{H_2O} = +\frac{|i_v|}{2F} M_{H_2O} \quad \text{in the CL} \qquad (12) $$

where $i_v$ and $F$ are the volumetric reaction current and Faraday's constant, respectively.

## 2.7. Species Transport in Gas

The mass flux of species in gas consists of convection and diffusion according to the Stefan-Maxwell derivation:

$$ \rho_g(u_g \cdot \nabla)\omega_i = \nabla \cdot \left[ \rho_g \omega_i \sum_k D_{ik}^{\text{eff}} \left( \nabla x_k + (x_k - \omega_k) \frac{\nabla p_g}{p_g} \right) \right] + R_i \qquad (13) $$

where $\omega_i$ and $R_i$ are the mass fraction and mass source term of species $i$, respectively. The mass source term is obtained by Equation (11) for oxygen (consumption) and Equation (12) for water (production). The effective binary diffusion coefficient between species $i$ and $k$ is denoted by $D_{ik}^{\text{eff}}$ and calculated by the Bruggeman correlation to consider the tortuosity of the porous media:

$$ D_{ik}^{\text{eff}} = D_{ik}^{\text{bulk}} \epsilon^{1.5} \qquad (14) $$

in which $D_{ik}^{\text{bulk}}$ is the bulk binary diffusion coefficient between species $i$ and $k$.

## 2.8. Electron and Ion Transport

The electrical current $i_e$ passing through the electronically conductive phase and the ionic current $i_i$ through the ionomer phase are both determined by Ohm's law:

$$ \nabla \cdot i_i = \nabla \cdot (-\sigma_i^{\text{eff}} \nabla \Phi_i) = i_v \quad \text{CL} \qquad (15) $$
<figure><img src="image_11.png" /></figure>

$$ \nabla \cdot i_e = \nabla \cdot (-\sigma_e^{\text{eff}} \nabla \Phi_e) = -i_v \qquad \text{CL, GDL} \qquad (16) $$

where $\Phi$ and $\sigma^{\text{eff}}$ represent the potential and effective conductivity, respectively. The subscripts i and e denote ionic and electronic, respectively. $i_v$ is the volumetric reaction current obtained during the oxygen reduction reaction. The effective ionic conductivity is obtained by GEIS fit (see *Electrochemical Characterizations* section and Table 1) while the effective electronic conductivity is calculated by the Bruggeman correlation:

$$ \sigma_e^{\text{eff}} = \sigma_e^{\text{bulk}} \varepsilon_e^{1.5} \qquad (17) $$

in which $\sigma_e^{\text{bulk}}$ is the bulk electronic conductivity and $\varepsilon_e$, the volume fraction of the electronically conductive phase (Table 1).

## 2.9. Oxygen Reaction Kinetics

The volumetric reaction current related to the oxygen reduction reaction within the catalyst layer is calculated using the Butler-Volmer equation:

$$ i_v = a_v i_0^{\text{ORR}} \left[ \exp\left(\frac{\beta F}{RT}\eta\right) - a_{H^+} a_{O_2} \exp\left(-\frac{\alpha F}{RT}\eta\right) \right] \qquad (18) $$

in which $\alpha$ and $\beta$ are the charge transfer coefficients for reduction and oxidation, respectively, both set to 0.5. $i_0^{\text{ORR}}$ refers to the exchange current density for the ORR and is derived from the GDE measurements (Table 1). In this equation, $a_{H^+}$ and $a_{O_2}$ represent the proton and oxygen activities, respectively. The overpotential $\eta$ is calculated from the electronic and ionic potentials:

$$ \eta = \Phi_e - \Phi_i - E_{\text{eq}} \qquad (19) $$

where $E_{\text{eq}}$ is the equilibrium potential, set to 1.23 V with respect to the counter electrode, which potential is set to 0 V.

## 2.10. Numerical Procedure

The governing equations are solved using COMSOL Multiphysics version 6.3. The computational mesh consists of ≈163 000 tetrahedral elements for the first and third cases, 298 000 for the second case, and 198 000 for the fourth case. A mesh refinement study was performed on the first case to ensure grid independence (see Figure S2, Supporting Information). Simulations are carried out on an Intel Core i7-10875H CPU running at 2.30 GHz. A relative tolerance of $10^{-3}$ is imposed for stationary solver steps. The segregated solution method is employed for the computation. All parameters, correlations, and model definitions are provided in Table 1, which also includes the applied boundary conditions.

## 3. Results and Discussion

**Figure 3a,b** presents the GEIS data obtained for a 0.6 µm-thick PtCo/VC catalyst layer at various applied currents. The GEIS spectra were verified by performing the Kramers-Kronig analysis, as shown in Figure S4 (Supporting Information). Across all currents, the high-frequency resistance remains unchanged, indicating no detectable variation in local conductivity within the $HClO_4$ solution due to proton consumption. Therefore, the high-frequency resistance value can reliably be used for ohmic drop correction along the polarization curve. At low frequencies, no mass-transport loop is observed, confirming that proton diffusion limitations within the $HClO_4$ solution are not present in these conditions. Additionally, the measured GEIS data do not exhibit a clear 45° line in the high-frequency range, as would be expected if ion transport was a strongly limiting phenomenon within the porous electrodes (see example Figure S1-9 in the Supporting Information of$^{[25]}$). This absence might be attributed to the small electrode thickness ($\delta_{\text{CL}}$ = 0.6 µm), and this was totally

<figure><img src="image_12.png" /><figcaption>**Figure 3.** a) GEIS spectra obtained during the ORR on a PtCo/VC GDE in 1 mol L$^{-1}$ at 1, 10, 25, and 50 mA and corresponding fits; b) zoom for the corresponding high-current measurements (10, 25, and 50 mA) and c) variation of the extracted charge transfer resistance ($R_{ct}$) values as a function of the current density for 3 representative replicates (using 3 different PtCo/VC GDEs).</figcaption></figure>
<figure><img src="image_13.png" /></figure>

intentional, hence the chosen ultra-low catalyst loading (20 µgPt cm⁻²). Interestingly, the low-frequency resistance appears to stabilize as the applied current exceeds 10 mA (140 mA cm⁻²) (see Figure 3c). This behavior suggests that O₂ mass-transport limitation within the catalyst layer initiates for currents above 10 mA. 3D-multiphysics simulations are used to gain a better understanding of these phenomena.

This study follows an experimental-to-model approach incorporating the geometric parameters of the experimental GDE into the simulations, such as channel dimensions, catalyst and gas diffusion layers thickness and porosity, PtCo loading. Some model parameters (given in Table 1) were directly derived from the experimental data. Specifically, the ECSA was estimated using CO stripping experiments performed in the GDE cell, while the effective ionic conductivity in the catalyst layer was determined from the GEIS data at 1 mA. The numerical results were then validated through direct comparison with experimental polarization curves. To gain deeper insights into the factors influencing the GDE operation, various simplified cases were explored in the simulations: 1) 1st case: dry catalyst layer, not flooded by liquid water and where the overlying electrolyte aqueous solution is not considered; 2) 2nd case: proton consumption in the active layer and diffusion in solution are the limiting factors; 3) 3rd case: fully flooded catalyst layer, assuming that its hydrophobicity is insufficient to counteract the presence of water produced by the ORR; 4) 4th case (mix between the second and third one): half of the active layer is flooded by water, while proton limitation occurs, due to consumption in the active layer and diffusion in solution.

All these four cases are summarized in Table 2 and illustrated in Figure 4. The hypotheses for each case will be thoroughly presented and described hereafter.

### 3.1. 1st Case: Dry Catalyst Layer

The first case represents the simplest scenario without considering the presence of the aqueous solution. As a result, the protonic activity is equal to unity ($a_{H^+} = 1$) in the Butler–Volmer law (Equation (18)). Additionally, the catalyst layer is assumed to be dry, meaning that oxygen exclusively diffuses in the gas phase, without any contribution from dissolved oxygen in liquid electrolyte. Consequently, the oxygen activity is equal to the following pressure ratio:

<table><caption>Table 2. Different scenarios explored in the simulations.</caption><thead><tr><th></th><th>Catalyst layer flooding</th><th>Proton transport limitation in solution</th><th>Comment</th></tr></thead><tbody><tr><td>1st case</td><td>No</td><td>No</td><td>Dry catalyst layer</td></tr><tr><td>2nd case</td><td>No</td><td>Yes</td><td>Proton limitation</td></tr><tr><td>3rd case</td><td>Yes</td><td>No</td><td>Fully flooded catalyst layer</td></tr><tr><td>4th case</td><td>Yes</td><td>Yes</td><td>1/2 dry | 1/2 flooded + proton limitation</td></tr></tbody></table>

$$a_{O_2} = \frac{p_{O_2}}{p_{\text{ref}}} \qquad (20)$$

where $p_{O_2}$ is the partial pressure of oxygen and $p_{\text{ref}}$, the reference pressure, set to 1 atm. Figure 5a compares the ohmic-drop corrected polarization curves obtained experimentally (dotted line) and numerically (solid line). The model accurately reproduces the activation kinetics behavior (>0.8 V) but clearly underestimates the voltage losses at higher currents (<0.8 V). This leads to a huge divergence between the experimental and the simulated curves at high current densities; the model predicts current densities as large as 4 A cm⁻² at 0.5 V when roughly 0.8 A cm⁻² are reached experimentally. Obviously, the experimental data do not support this simplistic 1st case model.

### 3.2. 2nd Case: Proton Limitation

The second case builds upon the previous one by accounting for proton consumption during the ORR. Consequently, a 55 µm-thick liquid electrolyte layer has been assumed; note that this proton diffusion layer thickness has been chosen to achieve the best fit to the measurements, and does not have any particular physical significance (Figure 4b). The thickness of the proton diffusion layer has been adjusted to improve the agreement between the experimental and modeled polarization curves. The proton transport in the liquid electrolyte and CL is described by a convection-diffusion equation:

$$\nabla \cdot (-D_{\text{HClO}_4}^{\text{eff}} \nabla c_{\text{H}^+}) + u_{\text{H}^+} \cdot \nabla c_{\text{H}^+} = R_{\text{H}^+} \qquad (21)$$

where $D_{\text{HClO}_4}^{\text{eff}}$ is the effective diffusion coefficient of perchloric acid, $c_{\text{H}^+}$, the proton concentration and $u_{\text{H}^+}$, the liquid electrolyte velocity. $R_{\text{H}^+}$ represents the proton consumption during the ORR:

$$R_{\text{H}^+} = - \frac{|i_v|}{F} \quad \text{in the CL} \qquad (22)$$

The proton activity thus depends on the proton concentration:

$$a_{\text{H}^+} = \frac{c_{\text{H}^+}}{c_{\text{H}^+, \text{ref}}} \qquad (23)$$

in which $c_{\text{H}^+, \text{ref}}$ is the reference proton concentration (Table 1). Figure 5b shows the results for this second case. The model again predicts well the activation kinetics as well as part of the ohmic drop (>0.65 V). However, this 2nd case leads to a limiting current caused by significant proton consumption at higher currents (<0.65 V), yielding substantial decrease in proton activity, as can be seen in Figure 6. This limiting current has not been observed experimentally, making this 2nd case rather irrelevant to account for the experimental data.

### 3.3. 3rd Case: Fully Flooded Catalyst Layer

The third case considers a catalyst layer completely flooded by the liquid electrolyte. This assumption arises from the absence of any membrane between the catalyst layer and the liquid electrolyte, leaving the relatively thin (0.6 µm) and rather hydrophilic catalyst
<figure><img src="image_14.png" /></figure>

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

All dimensions in µm.

<caption>**Figure 4.** Schematic of a 2D xz-plane slice from the three-dimensional modeling domain used in this work (Figure 2). The four cases described in Table 2 are shown: a) 1st case, b) 2nd case, c) 3rd case, and d) 4th case. An orange catalyst layer indicates a dry layer, while a blue catalyst layer signifies a fully flooded (3rd case) or partially flooded layer (4th case). Solid lines indicate subdomain boundaries: the red line denotes the GDL | MP interface, and the black line represents the CL | liquid electrolyte interface. For clarity, dimensions in the figure (given in µm) are not shown to scale.</caption>

layer (VC is not a very graphitic carbon black, hence is not very hydrophobic) in direct contact with the liquid electrolyte reservoir. Therefore, a convection-diffusion equation is used to calculate the dissolved oxygen concentration in the liquid electrolyte:

$$ \nabla \cdot \left( -D_{O_2/liquid}^{\text{eff}} \nabla c_{O_2} \right) + u_{O_2} \cdot \nabla c_{O_2} = R_{O_2} \quad (24) $$

where $D_{O_2/liquid}^{\text{eff}}$ is the effective diffusion coefficient of gaseous oxygen in water liquid being ≈25 000 times lower than that

in the gas phase; $u_{O_2}$ is the oxygen velocity in the liquid electrolyte. The oxygen activity is then described by the following concentration ratio:

$$ a_{O_2} = \frac{c_{O_2}}{p_{O_2}^{\text{in}} H_{O_2/liquid}} \quad (25) $$

where $p_{O_2}^{\text{in}}$ is the initial partial pressure of the oxygen and $H_{O_2/liquid}$, the Henry's constant (Table 1). Figure 5c presents
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 5.** Numerical (solid lines) and experimental (dotted line) polarization curves for the four cases studied. The cases considered are a) a fully dry active layer, 1st case, b) a limitation by the diffusion/consumption of the protons in the active layer, assuming a 55 µm-thick proton diffusion layer in the HClO₄ solution, 2nd case, c) a completely flooded active layer, 3rd case, and d) a mix between a half-flooded active layer with limitation by the diffusion/consumption of the protons, 4th case.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 6.** Calculated values of the local proton concentration in the 55 µm-thick liquid electrolyte layer for various cell voltage, obtained for the 2nd case. The abscissa 0 (liquid electrolyte) represents the interface with 1 mol L⁻¹ perchloric acid, while the abscissa 0 (active layer) corresponds to the electrolyte | active layer interface. The abscissa 0.6 µm (active layer) corresponds to the catalyst layer | gas diffusion layer interface. Note that the scale of thickness on the abscissa is not identical for the liquid electrolyte and the active layer (for which the scale is magnified, for better readability).</figcaption></figure>

the results for the third case. The O₂ diffusion in the water-flooded pores manifests itself by a significant increase of the ORR overpotential. As it can be seen, there is, this time, a great overestimation of the voltage losses in the simulated polarization curve compared to the experimental one, even at low current density (range between 50 and 200 mA cm⁻²). Whereas the

activation region matched rather well for the two previous cases, it is not true anymore for this fully-flooded assumption, except at very low current density (<50 mA cm⁻²). As the potential of the simulated polarization curve is lower than that of the experimental curve, it can be concluded that, under the current experimental conditions, the catalyst layer is never fully flooded.
This may be linked to the strong hydrophobicity of the GDL used
and the very small thickness of the catalyst layer (0.6 µm).

### 3.4. 4th Case: Partially Flooded Catalyst Layer and Proton Limitation

The last (and fourth) case considers a catalyst layer that is partially flooded by the liquid electrolyte. The first half of the layer (in contact with the GDL) is considered dry, meaning it behaves as in the first case, where oxygen diffuses in the gas phase. The second half of the layer (in contact with the liquid electrolyte) is considered flooded, meaning it behaves as in the third case, where oxygen is dissolved in the liquid electrolyte. A proton limitation similar to the second case has been added in a more restricted volume (10 µm-thick). Figure 5d presents the results for the fourth case. The model accurately reproduces the activation kinetics behavior (>0.75 V) and follows relatively well the voltage losses at higher currents (<0.75 V). A mass-transport limited current density appears at 900 mA cm⁻² even though no clear signs of proton transport limitation are visible on the experimental curve. This case does not represent a satisfying description of the experimental curves either, but highlights the possibility that a significant portion of the polarization curves can be accurately fitted by adjusting the proton mass-transport parameters in the solution.

### 3.5. Parametric Study: Flooding Levels within the Catalyst Layer

In summary, none of the fourth cases considered to date do satisfactorily fit the experimental GDE polarization plots. The previous results however suggest that, under the experimental conditions used, the catalyst layer is partially flooded. A parametric study based on the "flooded thickness" was conducted to assess whether this parameter alone is sufficient to replicate the experimental curve. Simulations were performed for incremental fractions of 1/6th of the flooded active layer, ranging from a fully dry state (0 × 1/6th) to a completely flooded active layer (6 × 1/6th). All these scenarios are presented in Figure 7.

The various simulations provide valuable insights into the potential effects of partial flooding in the active layer on the shape of the polarization curve. The graph clearly illustrates that the experimental curve cannot be reproduced under the assumption of static catalyst layer flooding level. As flooding gradually increases, the maximum overpotential increases, and the achievable current density decreases. However, the flooding process simultaneously exerts a negative impact in the activation region. Among the simulated cases, the best compromise is found when the fraction of the flooded active layer is 0.5, supporting the validity of the fourth scenario analyzed in the previous section. It is also highly likely that catalyst layer flooding dynamically increases with current density during the experiment, as the acidic ORR performed here is producing water.

Assuming that the shape of the experimental current-voltage curve is influenced by the progressive flooding of the catalyst layer, it is possible to estimate with the help of the simulations how the flooded thickness evolves with current density. Figure 8

<figure><img src="image_19.png" /><figcaption>Figure 7. Parametric study according to the flooded thickness from a fully dry catalyst layer to a fully flooded catalyst layer: numerical (solid lines) and experimental (dotted line) polarization curves.</figcaption></figure>

illustrates the evolution of the catalyst layer flooding as a function of current density to achieve the best agreement between the simulated and experimental polarization curve. At low current densities (below 100 mA cm⁻²), the catalyst layer remains essentially dry, corresponding to the initial simulation case with negligible impact of liquid water accumulation in the catalyst layer. As the current density increases, water production intensifies, leading to liquid water accumulation. Around 300 mA cm⁻², approximately half of the catalyst layer becomes saturated with liquid water, and this increases to two-thirds flooding by 400 mA cm⁻². Flooding continues to progress markedly with rising current density, reaching a level of about five-sixths of the layer near 700 mA cm⁻². Beyond this point, the rate of additional flooding slows significantly as the water reaches the interface with the hydrophobic gas diffusion layer, which is located at ≈0.6 µm in Figure 8. These results highlight the nonlinear flooding behavior with the current density. It likely results from several factors, including the microstructural properties of the catalyst layer—such as porosity, tortuosity, permeability, hydrophilicity, and thickness—as well as external operating conditions like temperature, pressure, gas flow, and humidification.

Finally, we would like to emphasize that a reasonable fit of the experimental polarization curves could also be achieved by incorporating a resistive contribution to the simulated curves. This scenario has been investigated and is shown in Figure S5 (Supporting Information); one can see that the fit between the experiment and modeled data is very good in that case too; however, it must be pointed out that it is not experimentally validated, as the GEIS data do not provide evidence of such an additional resistance. This shows that the best fit is not always the best result, and that the physical meaning of the model must stay a major driver in numerical simulations.
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 8.** Estimated catalyst layer flooding as a function of applied current density derived from the agreement between numerical simulations and the experimental polarization curve.</figcaption></figure>

## 4. Conclusion

The GDE is a powerful tool to assess the intrinsic ORR activity of catalysts at high current density/low potential, and so enables to better forecast their performance in PEMFC MEA. In this contribution, the ORR kinetics of a state-of-the-art PtCo/VC commercial catalyst was studied in a GDE in 1 mol L⁻¹ HClO₄ electrolyte at room temperature via cyclic voltammetry and galvanostatic electrochemical impedance spectroscopy. The experimentally obtained polarization curve is modeled using multiphysics and multicomponent 3D simulations, the model being partially fed by fitting of the EIS graphs. The model enables to investigate the interplay between the electrochemical kinetics and the reactants (O₂ and H⁺) and products (H₂O) mass-transports, that proceed in the flow channels of the MP, the GDL, the porous electrode, and the liquid electrolyte solution. Accurate measurement of the high-frequency resistance and ionic conductivity within the catalyst layer is of crucial importance to prevent misinterpretation of the processes involved and to ensure reliable simulation of the polarization curve. The simulations show that the active layer can neither be considered fully dry, nor fully flooded by the liquid electrolyte. Proton consumption in the active layer and diffusion in solution are also not the most limiting factors. Finally, the model highlights the non-negligible impact of partial flooding of the catalyst layer on the shape of the GDE polarization curves. Comparison of numerical simulations with experimental data enables to estimate the flooding thickness variation with the applied current density. This flooding is possibly caused by the accumulation of water dynamically produced in the course of the ORR; this is noteworthy, as it happens even for a very thin active layer as that used herein (0.6 µm). This shows that the oxygen mass-transport in the partially flooded CL plays a significant role even for small thicknesses of CL; it likely depends on the nature of the carbon support used (here Vulcan XC72), on its wettability and on the applied current density.

## Acknowledgements

Some of this work has been performed within the framework of the Centre of Excellence of Multifunctional Architectured Materials (CEMAM), Grenoble France n° ANR-10-LABX-44-01. The authors sincerely thank Y. Bultel for his precious support during scientific discussions.

## Conflict of Interest

The authors declare no conflicts of interest.

## Author Contributions

Erwan Tardy, Antoine Bonnefont, Florent Vandenberghe: software; validation; visualization; Erwan Tardy, Raphaël Riasse: writing—original draft preparation; Erwan Tardy, Raphaël Riasse, Florent Vandenberghe: investigation; Erwan Tardy, Raphaël Riasse, Florent Vandenberghe, Florence Druart, Marian Chatenet, Antoine Bonnefont: writing—review & editing. Erwan Tardy and Raphaël Riasse contributed equally.

## Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Keywords:** 3D multiphysics modeling · electrochemical impedance spectroscopy · gas diffusion electrode · numerical simulations · oxygen reduction reaction

[1] P. C. K. Vesborg, T. F. Jaramillo, RSC Adv. **2012**, 2, 7933.

[2] O. Le Bacq, A. Pasturel, R. Chattot, B. Previdello, J. Nelayah, T. Asset, L. Dubau, F. Maillard, *ChemCatChem* **2017**, *9*, 2324.
<figure><img src="image_22.png" /></figure>

[3] R. Chattot, T. Asset, J. Drnec, P. Bordet, J. Nelayah, L. Dubau, F. Maillard, *Nano Lett.* **2017**, *17*, 2447.

[4] L. Dubau, J. Nelayah, S. Moldovan, O. Ersen, P. Bordet, J. Drnec, T. Asset, R. Chattot, F. Maillard, *ACS Catal.* **2016**, *6*, 4673.

[5] C. Chen, Y. Kang, Z. Huo, Z. Zhu, W. Huang, H. L. Xin, J. D. Snyder, D. Li, J. A. Herron, M. Mavrikakis, M. Chi, K. L. More, Y. Li, N. M. Markovic, G. A. Somorjai, P. Yang, V. R. Stamenkovic, *Science* **2014**, *343*, 1339.

[6] L. Gan, M. Heggen, R. O'Malley, B. Theobald, P. Strasser, *Nano Lett.* **2013**, *13*, 1131.

[7] C. Cui, L. Gan, M. Heggen, S. Rudi, P. Strasser, *Nat. Mater.* **2013**, *12*, 765.

[8] Y. Garsany, J. Ge, J. St-Pierre, R. Rocheleau, K. E. Swider-Lyons, *J. Electrochem. Soc.* **2014**, *161*, F628.

[9] K. Kodama, T. Suzuki, K. Shinozaki, R. Jinnouchi, *J. Phys. Energy* **2023**, *5*, 014018.

[10] K. Kodama, T. Nagai, A. Kuwaki, R. Jinnouchi, Y. Morimoto, *Nat. Nanotechnol.* **2021**, *16*, 140.

[11] Y. Fukuyama, T. Shiomi, T. Kotaka, Y. Tabuchi, *Electrochim. Acta* **2014**, *117*, 367.

[12] T. A. Greszler, D. Caulk, P. Sinha, *J. Electrochem. Soc.* **2012**, *159*, F831.

[13] C. M. Zalitis, D. Kramer, A. R. Kucernak, *Phys. Chem. Chem. Phys.* **2013**, *15*, 4329.

[14] M. Inaba, A. W. Jensen, G. W. Sievers, M. Escudero-Escribano, A. Zana, M. Arenz, *Energy Environ. Sci.* **2018**, *11*, 988.

[15] B. A. Pinaud, A. Bonakdarpour, L. Daniel, J. Sharman, D. P. Wilkinson, *J. Electrochem. Soc.* **2017**, *164*, F321.

[16] R. Riasse, C. Lafforgue, F. Vandenberghe, F. Micoud, A. Morin, M. Arenz, J. Durst, M. Chatenet, *J. Power Sources* **2023**, *556*, 232491.

[17] F. Vandenberghe, F. Micoud, P. Schott, A. Morin, C. Lafforgue, M. Chatenet, *Electrochim. Acta* **2025**, *511*, 145364.

[18] T. Lazaridis, B. M. Stühmeier, H. A. Gasteiger, H. A. El-Sayed, *Nat. Catal.* **2022**, *5*, 363.

[19] G. K. H. Wiberg, M. Fleige, M. Arenz, *Rev. Sci. Instrum.* **2015**, *86*, 024102.

[20] R. Chattot, T. Asset, P. Bordet, J. Drnec, L. Dubau, F. Maillard, *ACS Catal.* **2017**, *7*, 398.

[21] K. Shinozaki, J. W. Zack, R. M. Richards, B. S. Pivovar, S. S. Kocha, *J. Electrochem. Soc.* **2015**, *162*, F1144.

[22] J. Huang, Y. Gao, J. Luo, S. Wang, C. Li, S. Chen, J. Zhang, *J. Electrochem. Soc.* **2020**, *167*, 166503.

[23] M. Eikerling, A. A. Kornyshev, *J. Electroanal. Chem.* **1999**, *475*, 107.

[24] C. R. Wilke, *J. Chem. Phys.* **1950**, *18*, 517.

[25] C. Lafforgue, P. Toudret, F. Micoud, M. Heitzmann, J.-F. Blachot, M. Chatenet, *J. Power Sources* **2025**, *627*, 235851.

[26] M. R. Gerhardt, L. M. Pant, A. Z. Weber, *J. Electrochem. Soc.* **2019**, *166*, F3180.

[27] R. Sander, *Atmos. Chem. Phys.* **2023**, *23*, 10901.

Manuscript received: April 28, 2025

Revised manuscript received: July 3, 2025

Version of record online: