# Characterizing PEM fuel cell catalyst layer properties from high resolution three-dimensional digital images, Part II: Oxygen effective diffusion and proton effective conductivity tensors

Mohamed Ahmed-Maloum, Joël Pauchet, Michel Quintard, Marc Prat

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

## To cite this version:

Mohamed Ahmed-Maloum, Joël Pauchet, Michel Quintard, Marc Prat. Characterizing PEM fuel cell catalyst layer properties from high resolution three-dimensional digital images, Part II: Oxygen effective diffusion and proton effective conductivity tensors. International Journal of Hydrogen Energy, 2024, 79, pp.394-410. 10.1016/j.ijhydene.2024.07.010. hal-04955037

HAL Id: hal-04955037

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

Submitted on 23 Oct 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.

Public Domain
# Characterizing PEM fuel cell catalyst layer properties from high resolution three-dimensional digital images, Part II: oxygen effective diffusion and proton effective conductivity tensors

Mohamed Ahmed-Maloumª, Joël Pauchetᵇ, Michel Quintardª, Marc Pratª

ªInstitut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse, CNRS –
Toulouse, France

ᵇCEA, LITEN, Université Grenoble Alpes, 38054 Grenoble, France

**Abstract:**

The oxygen effective diffusion and proton effective conductivity tensors of a cathode catalyst layer (CCL) are computed using a two-scale approach on CCL microstructure three-dimensional images obtained by focused ion beam-scanning electron microscopy (FIB-SEM) after reconstruction of the ionomer distribution. The results show that the obtained tensors are not isotropic. It is also found that some off-diagonal components are not negligible in the case of the proton effective conductivity tensor. The simulations confirm that Knudsen diffusion has a significant impact on oxygen diffusion. The results are comparable to those obtained in previous works on synthetic images as regards the oxygen effective diffusion through-plane component but differ as regards the proton conductivity through-plane component. Contrary to previous works using synthetic images which tend to underestimate the proton conductivity, they are globally in better agreement with experimental data. The results also show that the effective proton conductivity is heterogeneous at the scale of the computational domains typically considered in several previous works due to the variability in the ionomer distribution at this scale. The impact of liquid water in the primary pores on the proton transport is explored and found to have a significant impact.

**Keywords:** Ionomer, FIB-SEM, Numerical simulations, Effective diffusion tensor, Effective proton conductivity tensor, CCL, PEMFC
# 1. Introduction

Polymer Electrolyte Membrane Fuel Cells (PEMFC) are considered as a key technology in conjunction with the expected extensive use of clean and renewable energy. However, improvements in performance and durability are needed [1]. With the increase in the current densities in the operation of PEMFCs, a better understanding and characterization of the transport phenomena in the various porous layers forming the electrode become more important than ever in order to develop better designs [2, 3]. In complement to experimental tests, this can be achieved via modelling and numerical simulations [2, 6]. Modelling and simulations are also developed to reduce development cycles via systematic studies of the design and optimization of fuel cell systems [4-6]. As discussed in [4, 7, 8], the numerical simulations can be performed at various scales of interest. The direct numerical simulations (DNS) based on the Lattice Boltzmann Method [9, 10], the Volume of Fluid (VOF) method [11, 12], or finite element [13, 14] methods, or image-based methods, such as the full-morphology (FM) method [14, 15, 16], are performed over three-dimensional digital images of component microstructures. Pore network models have also been frequently used, e.g. [14, 16, 17, 18, 19] and references therein, to perform simulations at the scale of the microstructure. They are less accurate than DNS models but much less demanding in terms of computational cost. This allows to perform simulations over larger computational domains. Due to computational costs, however, PNM and DNS simulations are typically used for simulations over sub-domains of a single component, i.e. the catalyst layer or the gas diffusion layer for instance, and therefore no well-suited to perform simulations at the scale of the membrane electrode assembly (MEA) and beyond. For this reason, as exemplified for instance in [8], the detailed simulations of the many transport phenomena coupled to electro-chemistry phenomena occurring in a PEMFC are generally performed using macro-homogeneous models [4, 5, 20]. This type of models involves effective parameters, such as the permeability tensor or the effective diffusion tensor, e.g. [21], which are input data. A way of obtaining these input data and thus of bridging the scale between the microstructure properties and the macroscopic properties is to rely on DNS or PNM simulations [3, 9, 10, 14, 22-25]. In this work, such a numerical simulation approach to the characterization of effective transport properties is developed in the case of a cathode catalyst layer (CCL), focusing on two properties: the oxygen effective diffusion tensor in the gas phase and the effective proton conductivity tensor. As pointed out in [3], the proton conductivity has become a crucial parameter with the progress towards the operation of PEMFCs at ever higher current densities whereas the oxygen transport is related as well to
the ionomer distribution since the ionomer can partially block the pores and is also a transport pathway for the oxygen.

As discussed in [9, 22, 25, 26], one can distinguish two approaches for obtaining the 3D digital images of the CCL microstructure needed to perform the characterization of transport properties via DNS or PNM simulations. A possibility is to generate synthetic images according to a numerical reconstruction method [27]. Illustrations of this approach for the CCL can be found in [3], [10, 22-24, 28], among others. Instead of relying on a numerical generation of synthetic images, a presumably more attractive approach is to consider 3D images obtained via an imagery technique since it can be expected that the image so obtained is more representative of the real microstructure geometry and topology. Two main techniques have been used in this regard: nano-scale resolution X-ray computed tomography [29] and focused ion beam-scanning electron microscopy (FIB-SEM) [13, 19, 21, 26, 30-33], see [34] for a more general presentation of the technique. In the present paper, the simulations are based on FIB-SEM 3D images of a commercially available catalyst layer (TEC10E50E; Tanaka Kikinzoku Kogyo K.K., Japan) with a Pt loading of 46.2 wt%, an Ionomer/Carbon (I/C) ratio of 0.8 and Ketjen Black high surface area (HSA) carbon. However, as discussed in [26], the ionomer and the platinum loaded carbon phase cannot be distinguished in the segmented solid phase from the FIB-SEM data because of the lack of contrast in the SEM images between the ionomer and the carbon support due to their similar mass densities. Whereas identifying the ionomer in the image is not necessarily mandatory as regards the characterization of the oxygen transport in the segmented pore space, this is of course essential in order to characterize the proton transport across the CCL. For this reason, as detailed in the companion paper [26], a procedure was developed to reconstruct the ionomer distribution in the FIB image. In what follows, the images presented in [26] after reconstruction of the ionomer in the images are used. As mentioned in [26], the approach used here to obtain the CCL microstructure can therefore be considered as hybrid [35], since it combines an imagery technique (FIB-SEM) with a numerical approach to reconstruct in the image the ionomer localization, which cannot be directly seen in the FIB-SEM segmented image.

As mentioned above, once 3D images of the CCL microstructure have been obtained, the next step is to perform numerical simulations of the transport phenomena of interest directly at the pore-scale, that is over the 3D images. In this work, this is performed using a finite volume discretization method within the framework of the open source CFD software OpenFOAM [36] to solve the oxygen transport equation in the gas phase and the proton
transport equation in the ionomer to obtain the respective effective tensors.

The paper is organized as follows. Images from [26] are presented in Section 2 together with a summary of the method developed in the companion paper [26] to reconstruct the various solid components, namely the ionomer, the platinum nanoparticles and the carbon in the FIB-SEM segmented image. The numerical methodology used to compute the oxygen effective diffusion tensor and the proton effective conductivity tensor is presented in Section 3. Results of simulations are presented in Section 4. Section 5 presents the main conclusions of the study.

## 2. Cathode catalyst layer microstructure

As mentioned above and described in [26], a digital image of the CCL microstructure is obtained by FIB-SEM with a voxel size of 5 × 5 × 5 nm³. The obtained segmented image is illustrated in Fig.1. The thickness of the CCL sample shown in Fig.1 is comparable to the actual CCL thickness (~ 6 µm). Details on the image processing used to obtain the segmented image can be found in [26].

<figure><img src="image_2.png" /><figcaption>Fig.1 - CCL segmented microstructure from FIB-SEM image with segmented solid phase in light gray and pore space in black. The voxel size is 5 × 5 × 5 nm³, y is the through-plane direction, z is the FIB-SEM milling direction.</figcaption></figure>

As mentioned above, the ionomer cannot be distinguished from the carbon phase in the segmented image. Also, the platinum nanoparticles cannot be distinguished since their average diameter (3-4 nm) is smaller than the voxel size (5 nm). As a result, both the ionomer and the platinum nanoparticles are reconstructed in the segmented image using
numerical techniques. The method for distributing the Pt nanoparticles in the images is presented in [26]. However, the identification of platinum particles in the images is not needed for the present work. For this reason, the Pt particle reconstruction procedure is not presented again in the present article. The main features of the ionomer reconstruction are now recalled. Again, the details can be found in [26]. Three types of ionomer structures, namely the filament ionomer, the coating ionomer and the inter-carbon grain ionomer, are distinguished. The ionomer filaments are identified by applying a filtering procedure to the distribution of solid inclusions obtained via a maximum sphere inscription algorithm available in the image analysis toolbox Porespy [37]. The ionomer coating is identified via a filtering procedure of the curvature distribution computed over the segmented solid surface considering that the concave regions at the surface are preferentially coated during the drying step of the CCL fabrication procedure. This procedure of coating controlled by the local curvature is applied until a desired value of coverage is reached. The proposed procedure leads to a coating ionomer film of variable thickness in the range [2.5 -10 nm] in agreement with the available experimental data [38-41].

The coating ionomer distribution is illustrated in Fig. 2 over a sub-region of the CCL of 250 x 250 x 250 nm³. Since the voxel size in the original image is of 5 x 5 x 5 nm³, these voxels are subdivided into voxels of size 2.5 x 2.5 x 2.5 nm³ in order to obtain a coating ionomer minimum film thickness of 2.5 nm. As illustrated in Fig.2, the coating by the ionomer is not perfect and some regions of the carbon particle surface are not covered. This is consistent with experimental studies showing that the coverage can vary depending on the I/C ratio, e.g. [38] and references therein. Consistent with the experimental observations, [38-40], the convex zones of the carbon agglomerate surface are less likely to be covered than the concave zones.

As pointed out in [26], the obtained coating ionomer forms a well-connected percolating region in the range of I/C ratios classically considered for CCL. However, as also pointed out in [26], the reconstructed filament ionomer and coating ionomer do not correspond to all the volume (or mass) of ionomer expected from the I/C ratio (0.8 in the case of the considered CCL). The complementary ionomer fraction is referred to as the inter-carbon grain ionomer. The inter-carbon grain ionomer is not explicitly reconstructed. It corresponds to ionomer deeper inside the FIB-SEM segmented solid phase (see [26] for more details) in between the carbon grains.
<figure><img src="image_3.png" /><figcaption>Fig.2 - a) 3D rendering of FIB-SEM image segmented solid phase. b) ionomer filaments (in green) and ionomer coating (in blue). The complementary solid phase is in grey; c) Example of coating ionomer film thickness distribution. All images are for a CCL sub-region of 250 x 250 x 250 nm³.</figcaption></figure>

Thus, after reconstruction of the ionomer coating and the filaments, the FIB-SEM image segmented solid phase can be decomposed into the explicitly reconstructed ionomer (I<sub>R</sub>) consisting of the coating and filament ionomer and the complementary solid phase (CSP). The latter is composed of inter-carbon-grain ionomer, platinum nanoparticles, carbon grains and pores. The Ketjen Black carbon particles in the CCL are indeed nanoporous particles with a porosity of about 0.25 [43] and micropores less than a few nanometers [43]. The carbon grain porosity cannot be detected in the FIB-SEM images due to the voxel resolution of 5 nm. The transport of oxygen in what follows is therefore considered in the pore space
identified from the FIB-SEM image, i.e. over a range of significantly bigger pore sizes. As discussed in [44], it is usual to distinguish the primary pores inside the carbon particles (1-10 nm in radius) and the secondary pores (10-100 nm in radius) between the carbon particles. The pore sizes identified in [26] by using the maximum sphere inscription algorithm available in Porespy [37] vary between 15 and 200 nm with an average value around 60 nm. Thus, the pore space obtained via the FIB-SEM imagery technique corresponds to the secondary pores whereas the pores in the CSP are referred to as the primary pores. The secondary pore space is referred to as the macro-pore space in the present article and the volume fraction of macro-pores in the segmented FIB-SEM image is the macroporosity. The latter is denoted by $ε$. Furthermore, it is often considered that the primary pores are filled with water under usual PEMFC operating conditions [44]. Therefore, as in many previous works, [9] and references therein, the oxygen transport is assumed negligible in the primary pores compared to the oxygen transport in the secondary pores in what follows.

### 3 Computation of oxygen and proton effective transport properties.

Following the above discussion, the phase repartition in the CCL is schematically represented in Fig. 2. In what follows, letter $V$ denotes the pore space in the FIB-SEM segmented image (i.e. the secondary pores, indicated as “void” in Fig.2), referred to as the macro-pore space, $I_R$ is the reconstructed ionomer (coating and filaments), CSP is the complementary solid phase. As discussed above, the CSP is in fact a ternary system made of porous carbon (C), Pt nanoparticles, and inter-carbon grain ionomer, $I_{ICG}$. In this section, we present the various models that are used to estimate the effective diffusion of oxygen in the macro-pores, and the effective conductivity of protons in the ionomer.

#### 3.1 Upscaling Oxygen and proton transport equations.

The complete development of the macro-scale equations for the considered transport processes is not presented again here. One can refer to [45] for a general presentation of effective transport using a volume averaging theory, and to [21] for an application to the oxygen diffusion transport in gas diffusion layer (GDL) with a micro-porous layer (MPL). However, some indications are provided to clarify the computations made in this work.

##### Oxygen Transport

Oxygen transport by diffusion at the pore scale in the pore space is most often described by a Fickian diffusion equation
$$ \frac{\partial c_{O_2}}{\partial t} = \nabla \cdot (D_{O_2}(\mathbf{x})\nabla c_{O_2}) \qquad (1) $$

where $c_{O_2}$ is the oxygen concentration, $\mathbf{x}$ denotes a position vector, and $D_{O_2}$ is the local oxygen diffusion coefficient. Neglecting for simplicity the interfacial resistance [46], the boundary condition at the secondary pore-space and reconstructed ionomer interface, $A_{V-I_R}$, is

$$ c_{IO_2} = c_{O_2} RTH_e \qquad (2) $$

where $H_e$ is the Henry's constant for oxygen dissolved in the ionomer, $T$ is the temperature, $R$ is the gas constant, and $c_{IO_2}$ is the oxygen concentration in the ionomer phase [13]. A zero flux boundary condition can be applied at interfaces $A_{V-CSP}$ if the very small flux potentially coming into the $I_{ICG}$-phase is neglected. In other words, we make the assumption that the electrochemical reaction does not have an impact on the effective diffusion tensor. This problem of transport of a chemical species in a multiphase system in a porous medium with partitioning coefficients has been studied using upscaling techniques, see for instance [47] and references therein.

The upscaled equation associated with Eq. (1) may be written,

$$ \frac{\partial \varepsilon c_{O_2}}{\partial t} = \nabla \cdot (\varepsilon \mathbf{D}_{O_2 \text{eff}}(\mathbf{x}) \cdot \nabla c_{O_2}) - q \qquad (3) $$

where, as indicated before, $\varepsilon$ is the CCL macro-porosity (i.e. the porosity corresponding to the secondary pores), $\mathbf{D}_{O_2 \text{eff}}$ is the effective diffusion tensor, $C_{O_2}$ is the intrinsic average concentration over the averaging volume, and $q$ is the averaged flux exchanged with the $I_R$-phase, see [45] for more details. The theory links $\mathbf{D}_{O_2 \text{eff}}$ and $q$ to the transport problem in the $I_R$-phase [47]. As discussed for instance in [1], [13] or [46], oxygen can dissolve in the ionomer and be transported within the ionomer. This transport process is also frequently represented via a diffusion equation, i.e. according to an equation similar to Eq.(1), e.g. [13], [46]. However as pointed out in [13], the oxygen diffusivity in the ionomer is several orders of magnitude smaller than in the pore space. Also, the porosity of the CCL is 0.55 whereas the reconstructed ionomer volume fraction is 0.124 for the considered I/C ratio of 0.8 [26]. As the result, $O_2$ transport across the CCL via the ionomer pathway can be expected to be negligible compared to the diffusive transport in the pore space, except perhaps when the presence of liquid water in the pore space is sufficiently high for hampering this transport, i.e. when the gas phase is about to cease to percolate. In any case, $O_2$ transport in the ionomer is neglected in the present paper, which therefore focuses on the characterization of $O_2$ transport in the gas phase in the secondary pores. Assuming that the impact of the exchange
with the ionomer on $D_{O_2eff}$ is negligible, it is an acceptable approximation to assume that $D_{O_2eff}$ can be computed like in the case of passive diffusion, i.e., assuming zero flux at interfaces. Under these circumstances, the upscaled equation (Eq.(3)) to be considered for determining the effective diffusion tensor can be reduced to

$$ \nabla \cdot (\varepsilon D_{O_2eff}(\mathbf{x}) \nabla C_{O_2}) = 0 \qquad (4) $$

and the pore-scale problem to be solved over a volume of the CCL to estimate the effective diffusion tensor becomes the following steady-state problem:

$$ \nabla \cdot (D_{O_2}(\mathbf{x}) \nabla C_{O_2}) = 0 \qquad (5) $$

$$ \mathbf{n} \cdot (D_{O_2}(\mathbf{x}) \nabla C_{O_2}) = 0 \quad \text{at } A_{V-SSP} \qquad (6) $$

where $\mathbf{n}$ is the normal at the interface and $A_{V-SSP}$ is the interface between the (macro-)pore space and the segmented solid phase (void – FIB-SEM solid phase interface in Fig. 2a).

Given the pore sizes in the CCL (between 15 nm and 200 nm with $\ell$ of an average value of 60 nm as mentioned in Section 2), the Knudsen number $Kn = \frac{\ell}{d_p}$, where $d_p$ is the pore diameter and $\ell = \frac{RT}{\sqrt{2\pi d_{mol}^2 N_A P}}$ is the mean free path ($R$ is the ideal gas constant, $T$ is the temperature (353 K), $P$ is the pressure (1.5 bar), $d_{mol}$ is $O_2$ molecular diameter and $N_A$ is the Avogadro number) ranges between 0.4 and 5.4 with a value of 1.3 for the average pore diameter of 60 nm. For such high values of the Knudsen number, Knudsen diffusion [48] must be considered in conjunction with Fick's diffusion. As in several previous works, e.g. [9, 49] and references therein, our approach is based on the Bosanquet relation [50] combining the molecule-wall and molecule-molecule interactions according to the formula,

$$ D_{O_2}(\mathbf{x}) = \frac{1}{D_{O_2\text{bulk}}^{-1} + D_{O_2-\text{Kn}}^{-1}} \qquad (7) $$

where $D_{O_2\text{bulk}}$ is the bulk diffusion coefficient and $D_{O_2-\text{Kn}}$ is the Knudsen diffusion coefficient. The latter is expressed as

$$ D_{O_2-\text{Kn}} = \frac{d_p}{3} \sqrt{\frac{8RT}{\pi M}} \qquad (8) $$

where $M$ is the molar mass. As mentioned before, the pore diameter $d_p$ in Eq.(8), is determined using a maximum sphere inscription algorithm [21].

In a fuel cell under operation, oxygen and water vapor both diffuses together in the pore space. A more comprehensive description then implies to consider the gas phase as a more complex multicomponent system [51] than the binary system implicitly considered in our approach. If we consider a mixture of the main components ($O_2$, $N_2$, $H_2O$ vapor), the more general diffusion model is the Maxwell-Stefan description [52], which can be recast into a
generalized Fick's diffusion law involving a matrix of diffusion coefficients multiplied by the vector of concentration gradients. Using this Maxwell-Stefan model, it has been demonstrated in [53] that the macro-scale effective diffusion can be expressed as a macro-scale generalized Fick's law, in which the micro-scale matrix of diffusion coefficients is simply multiplied by the inverse of a tortuosity tensor, which does not depend on the diffusion coefficients. For instance, for an isotropic medium, we can use the same tortuosity coefficient, and it is equal to the one obtained in the case of a binary system. As a consequence, the results presented in the present paper can be equally used in a more general multicomponent description. Thus, the present work can be considered as a first step toward a more comprehensive approach. It can also be noted that the binary system approximation has been considered in several previous works and is actually the most common approach in the fuel cell literature. The consideration of the binary system approximation in the present work allows comparison with some of these works as reported in Section 4.

## Proton Transport

The problem of proton transport in the polymer is not a problem of passive diffusion, i.e., with no interaction with the other phases, in particular because of the reaction with oxygen at the catalytic sites. Like for the oxygen transport, mass exchange at interfaces may impact the mathematical structure of the equations, as well as the effective properties. As illustrated in [54], tortuosity may be affected in the case of strong reaction (i.e., large Damkhöler values). If this effect can be neglected, the tortuosity tensor can be calculated like in the case of passive diffusion. While this certainly requires further investigation, e.g. pore-scale detailed simulations of all coupled transport mechanisms, we make, as in several previous works [3, 9] and references therein, such an assumption.

Hence, the governing equations for the proton transport in the ionomer, necessary to estimate the effective conductivity tensor, is also expressed via a diffusion equation [9], as

$$ \nabla \cdot (\sigma_p(x) \nabla \varphi_I) = 0 \qquad (9) $$

where $\sigma_p(x)$ is the proton conductivity in the ionomer and $\varphi_I$ is the potential in the ionomer [23]. This equation is valid within the reconstructed ionomer $I_R$ and the inter-carbon grain ionomer $I_{ICG}$ distributed within the CSP-phase as well. The problem is completed with a zero flux condition at the interface with the pore space (V) and the carbon phase (C).

As reported for instance in [55], $\sigma_p$ strongly varies with the relative humidity, i.e. the ionomer water uptake $\lambda$. Also, as pointed out in [1], specifying the conductivity in Eq. (9) can be an issue since it has been shown that the conductivity of ionomer thin films differs from
that of the membranes [56, 57]. As in several previous works, e.g. [9] and references therein, a simplified approach is considered, assuming that the ionomer conductivity is spatially constant over the computational domain. In other words, the possible variations of the ionomer conductivity with the ionomer coating local film thickness are ignored. Also, the possible presence of liquid water can make the determination of the proton effective conductivity even more difficult since protons can be transported within the liquid water in regions free of ionomer [58, 59]. This effect is first neglected. Nevertheless, insights on the possible impact of liquid water in the primary pores on the proton transport are presented in Section 4. The upscaled proton transport governing equation, associated with Eq. (9), is expressed as,

$$ \nabla \cdot (\varepsilon_I \boldsymbol{\sigma}_{\text{Peff}} \cdot \nabla \phi_I) = 0 \qquad (10) $$

where $\varepsilon_I$ is the volume fraction of ionomer (filament, coating and inter-carbon grain ionomer) in the computational domain, $\boldsymbol{\sigma}_{\text{Peff}}$ is the proton effective conductivity tensor and $\phi_I$ is the intrinsic average potential.

In what follows, we will first consider the effective conductivity of the reconstructed ionomer, namely Eq. (11),

$$ \nabla \cdot (\varepsilon_{IR} \boldsymbol{\sigma}_{\text{Peff}-IR} \cdot \nabla \phi_I) = 0 \qquad (11) $$

where $\varepsilon_{IR}$ is the volume fraction of reconstructed ionomer (the filament and coating ionomer).

### Alternative notation and tortuosity definition

In the proposed macro-scale equations, the volume fraction is taken outside the effective diffusion tensor. It is to be noticed that another definition, sometimes adopted in some papers, incorporate the volume fraction into the effective diffusion tensor. In this case, Eq. [3] for instance is expressed as

$$ \frac{\partial \varepsilon C_{O_2}}{\partial t} = \nabla \cdot (\mathbf{D}_{O_2 \text{eff}}^* (\mathbf{x}) \cdot \nabla C_{O_2}) - q \qquad (12) $$

with

$$ \mathbf{D}_{O_2 \text{eff}}^* = \varepsilon \mathbf{D}_{O_2 \text{eff}} \qquad (13) $$

and, similarly, the effective proton diffusion tensor is defined by

$$ \boldsymbol{\sigma}_{\text{Peff}}^* = \varepsilon_I \boldsymbol{\sigma}_{\text{Peff}} \qquad (14) $$

and $\boldsymbol{\sigma}_{\text{Peff}-IR}^* = \varepsilon_{IR} \boldsymbol{\sigma}_{\text{Peff}-IR}$ when only the reconstructed ionomer is considered. Also, it is customary to discuss the components of these tensors in dimensionless form using the concept of tortuosity. For instance, the dimensionless through-plane component of oxygen
effective diffusion tensor is expressed as

$$ \frac{D_{O_2effyy}^*}{D_{O_2bulk}} = \frac{\epsilon D_{O_2effyy}}{D_{O_2bulk}} = \frac{\epsilon}{\tau_{O_2yy}} \quad (15) $$

where $D_{O_2bulk}$ is the molecular diffusion coefficient of oxygen for the considered pressure and temperature and $\tau_{O_2yy}$ is the tortuosity for the considered tensor component. Similarly,

$$ \frac{\sigma_{Peffyy}^*}{\sigma_P} = \frac{\epsilon_I \sigma_{Peffyy}}{\sigma_P} = \frac{\epsilon_I}{\tau_{Pyy}} \quad (16) $$

and

$$ \frac{\sigma_{Peffyy-I_R}^*}{\sigma_P} = \frac{\epsilon_I R \sigma_{Peffyy-I_R}}{\sigma_P} = \frac{\epsilon_I R}{\tau_{PI_Ryy}} \quad (17) $$

## 3.2 Numerical methodology for the computation of effective properties (sequential approach)
As exemplified in [21] in the case of a MPL, the method for determining the effective transport properties consists in solving local problems, e.g. the pore-scale boundary-value problems (BVP) expressed by Eqs. (5) and (9) together with appropriate boundary conditions. When working with periodic images, periodicity conditions are logically used to make the BVP problem complete. When the medium image is not periodic, different boundary conditions must be used at the computational domain boundaries. This is discussed for instance in [60].

As detailed below, we use the so-called permeameter conditions. The effective transport coefficients are obtained by (i) solving pore-scale equations over a given block of the image, (ii) solving the macro-scale equation over the same domain, with the same boundary conditions, (iii) identifying the effective diffusion tensor by requiring that the macroscopic fluxes obtained by averaging the pore-scale fluxes from step (i) calculation should be equal to their macroscopic counterparts (this is detailed below).

However, solving Eqs. (5) and (9) directly at the scale of the whole CCL image (shown in Fig.1) is not efficient numerically due to the corresponding quite large number of grid points.

As discussed in details in [21], it is more efficient to perform the resolution of the transport problems at the CCL scale by a sequential grid block technique. This sequential approach is schematically represented Fig. 3. The grid block technique consists in partitioning the CCL matrix initial digital image in adjacent cubic blocks of equal size $L_i$. The pore-scale governing equations, Eqs. (5) and (9), are solved for each grid block with appropriate boundary conditions, and effective transport properties are identified for each block. Then, the macroscopic equations (4) and (10) (or (11)) are solved over the coarser grid corresponding to the grid blocks to obtain the effective transport properties at a larger scale,
denoted by $L_2$, following the same identification procedure.

<figure><img src="image_4.png" /><figcaption>Fig. 3 - Schematic representation of the sequential approach for calculating effective properties</figcaption></figure>

### 3.3 Effective transport properties over a grid block (scale $L_1$).

The method for estimating the effective diffusion tensor for a given grid block is based on the following steps:

1. Perform small-scale simulation on the grid block with given boundary conditions,

2. Perform macro-scale simulations on the same geometry, with the same boundary conditions, but this case with a homogeneous diffusion tensor. The effective diffusion tensor is estimated so as to reproduce the average fluxes through the block as calculated in step 1.

The procedure is detailed below. Dirichlet conditions are used in one direction with Bamberger conditions [61] for lateral conditions. For instance, for imposing a macro-scale gradient in the Ox direction, the used boundary conditions read

$$c_{O_2} = c_{O_2-1} \text{ at cube inlet } (x=0)$$

$$c_{O_2} = c_{O_2-2} \text{ at cube outlet } (x = L_1)$$

$$c_{O_2} = c_{O_2-1} + (c_{O_2-2} - c_{O_2-1}) \frac{x}{L_1} \quad \text{at lateral walls} \quad (18)$$

where the origin for x-axis is taken at the inlet face of the block.

Zero flux conditions are imposed on the internal limiting boundaries via a penalization technique (see again [21]), i.e. a continuous formulation of the transport equation with a very small value of the diffusion coefficient in the inert phases.
For oxygen transport, Eq.(5) is solved together with the above-discussed boundary conditions using the LaplacianFoam solver of the open source CFD software OpenFOAM [36].

Once the pore scale concentration field $c_{O_2}(\mathbf{x})$ is obtained, the average of the diffusive flux over the computational domain is computed as

$$ \langle j_{O_2} \rangle = \frac{1}{\Omega} \int_{\Omega_f} (-D_{O_2}(\mathbf{x}) \nabla c_{O_2}) d\Omega \quad (19) $$

where $\Omega$ is the volume of the computational domain and $\Omega_f$ is the volume occupied by the pore space in $\Omega$.

Then, it is expressed that the average flux given by Eq. (19) must be equal to the average of the macroscopic diffusion flux. The macroscopic flux could be obtained by solving the 3D problem for an unknown tensor $D_{O_2eff}$, this latter being identified through some optimization algorithm. To avoid such computations, it is often admitted, as in [21], that the macro-scale concentration field can be approximated by the solution of the macro-scale problem over a box of infinite lateral extent (this is exact in the case of Bamberger boundary conditions and a homogeneous medium). In this case, the macroscopic flux is constant and the equivalence condition reads

$$ J_{O_2} = -\varepsilon D_{O_2eff} \nabla C_{O_2} = \langle j_{O_2} \rangle \quad (20) $$

The macroscopic gradient is given by $\nabla C_{O_2} = \frac{c_{O_2-2}-c_{O_2-1}}{L_1} e_x$, where $e_x$ is the unit vector in the direction of $x$ and $L_1$ is the dimension of the computational domain in the $x$ direction. Using Eq.(20) leads to the following relationships

$$ D_{O_2effxx} = -\frac{L_1}{\varepsilon(c_{O_2-2}-c_{O_2-1})} \langle j_{O_2x} \rangle \quad (21) $$

$$ D_{O_2effyx} = -\frac{L_1}{\varepsilon(c_{O_2-2}-c_{O_2-1})} \langle j_{O_2y} \rangle \quad (22) $$

$$ D_{O_2effzx} = -\frac{L_1}{\varepsilon(c_{O_2-2}-c_{O_2-1})} \langle j_{O_2z} \rangle \quad (23) $$

Applying the same procedure in the $y$ and $z$ directions allows computing all the components of tensor $D_{O_2eff}$,

$$ \begin{pmatrix} D_{O_2effxx} & D_{O_2effxy} & D_{O_2effxz} \\ D_{O_2effyx} & D_{O_2effyy} & D_{O_2effyz} \\ D_{O_2effzx} & D_{O_2effzy} & D_{O_2effzz} \end{pmatrix}_{L_1} = -\frac{L_1}{\varepsilon_{L_1}(c_{O_2-2}-c_{O_2-1})} \begin{pmatrix} \langle j_{O_2x} \rangle^1 & \langle j_{O_2x} \rangle^2 & \langle j_{O_2x} \rangle^3 \\ \langle j_{O_2y} \rangle^1 & \langle j_{O_2y} \rangle^2 & \langle j_{O_2y} \rangle^3 \\ \langle j_{O_2z} \rangle^1 & \langle j_{O_2z} \rangle^2 & \langle j_{O_2z} \rangle^3 \end{pmatrix} \quad (24) $$

where indices 1, 2 and 3 correspond to the application of the inlet-outlet concentration difference ($c_{O_2-1} - c_{O_2-2}$) in the $x$, $y$ and $z$ directions, respectively. At this point, it must be
emphasized that the resulting tensor is not necessarily symmetric, since symmetry can only be demonstrated exactly in the case of the BVP with periodicity conditions. In addition to the variability in the distribution of effective values, the possible deviation to symmetry characterizes how much the considered grid-blocks are far from a representative elementary volume.

### 3.4 Effective transport properties at the larger scale (scale $L_2$)

Once the effective tensor $D_{O_2eff-L_1}$ has been computed for each grid block of size $L_1$ according to the procedure described in Section 3.3, the diffusion problem at the larger scale $L_2$ is then formulated as (steady-state, no exchange with the surrounding phases)

$$ \nabla \cdot (\varepsilon_{L_1}(\mathbf{x}) D_{O_2eff-L_1}(\mathbf{x}) \cdot \nabla C_{O_2}) = 0. \qquad (25) $$

Such a diffusion problem can also be uspcaled into a larger-scale diffusion problem, and methods similar to the $L_1$ upscaling can be used to estimate the larger-scale effective diffusion tensor $D_{O_2eff-L_2}$.

In this approach, the effective diffusion properties determined for one block and its porosity are assigned to each point in the considered block. Then the procedure to determine the effective diffusion tensor at scale $L_2$ is similar to the one described in Section 3.3. Eq. (25) is solved numerically using the same boundary conditions as presented in Section 3.3 considering successively the directions x, y and z and the effective tensor components at scale $L_2$ are determined using relationships similar to Eq. (24), namely

$$ \begin{pmatrix} D_{O_2effx} & D_{O_2effxy} & D_{O_2effxz} \\ D_{O_2effyx} & D_{O_2effyy} & D_{O_2effyz} \\ D_{O_2effzx} & D_{O_2effzy} & D_{O_2effzz} \end{pmatrix}_{L_2} = - \frac{1}{\varepsilon_{L_2}(c_{O_2-2} - c_{O_2-1})} \begin{pmatrix} H_x \langle j_{O_2x} \rangle^1 & H_y \langle j_{O_2x} \rangle^2 & H_z \langle j_{O_2x} \rangle^3 \\ H_x \langle j_{O_2y} \rangle^1 & H_y \langle j_{O_2y} \rangle^2 & H_z \langle j_{O_2y} \rangle^3 \\ H_x \langle j_{O_2z} \rangle^1 & H_y \langle j_{O_2z} \rangle^2 & H_z \langle j_{O_2z} \rangle^3 \end{pmatrix} \qquad (26) $$

where $H_x = 10.5$ µm, $H_y = 12$ µm, $H_z=4.5$ µm (Table 1) and

$$ \langle j_{O_2} \rangle = \frac{1}{\Omega} \int_{\Omega} (-\varepsilon_{L_1}(\mathbf{x}) D_{O_2eff-L_1}(\mathbf{x}) \cdot \nabla C_{O_2}) d\Omega \qquad (27) $$

As for the first step (Section 3.3), Eq. (25) is solved using the LaplacianFoam solver in OpenFOAM [36].

It can be noted that characterizations of the full tensor are scarce in the literature. In most papers on CCL, only the through plane effective diffusion coefficient is discussed.

The above-described computational procedure, which leads to determine the full tensor, is applied to the characterization of the oxygen effective transport properties considering that the macro-pores are dry and using the 5 nm FIB-SEM voxels as elementary finite volumes
for the simulations at scale $L_1$. This is reasonable since the average pore size (60 nm) is more than 10 times bigger than the voxel size.

### 3.5 *Through-plane proton effective conductivity.*

The situation is different as regards the proton transport since the coating ionomer film thickness can be as low as 2.5 nm (Fig.2). It is considered that a minimum of three computational voxels over the ionomer film thickness is necessary to perform the simulations at scale $L_1$ in the case of the proton transport simulations. This leads to subdivide each FIB-SEM voxel into 9 cubic computational voxels of 0.833 nm of side. This increases the total number of voxels at the scale of each grid block as well as at the whole scale of the CCL FIB-SEM image by a factor of about 10, which increases significantly the computational burden. For this reason, scale $L_2$ is smaller in the case of the proton conductivity characterization than for the $O_2$ transport simulations. The various computational domain sizes considered for each case are summarized in Table 1.

<table><caption>**Table 1** - Computational domain sizes and cubic finite volume sizes considered in the simulations</caption><thead><tr><th></th><th>L<sub>1</sub></th><th>L<sub>2</sub></th><th>Computational<br>finite volume size</th></tr></thead><tbody><tr><td>O<sub>2</sub></td><td>1500 nm ×1500<br>nm ×1500 nm</td><td>7 L<sub>1</sub> × 3 L<sub>1</sub> ×8 L<sub>1</sub><br>(10.5 µm × 4.5 µm × 12 µm)</td><td>5 nm (L<sub>1</sub>)<br>L<sub>1</sub> /10 (L<sub>2</sub>)</td></tr><tr><td>H+</td><td>250 nm ×250<br>nm ×250 nm</td><td>6 L<sub>1</sub> × 6 L<sub>1</sub> ×6 L<sub>1</sub><br>(1.5 µm × 1.5 µm × 1.5 µm)</td><td>0.833 nm (L<sub>1</sub>)<br>L<sub>1</sub> /10 (L<sub>2</sub>)</td></tr></tbody></table>

As can be seen from Table 1, the proton conductivity is computed at scale $L_2$ considering a cubic computational domain of size 1500 nm ×1500 nm ×1500 nm whereas the scale $L_2$ is close to the whole FIB-SEM image (Fig.1) in the case of the oxygen transport computation.

## 4. Oxygen and proton effective transport properties: results of simulations.

### 4.1 *Oxygen effective diffusion tensor*

Simulations are performed considering a total gas pressure of 2 bars and a temperature of 80°C. Fig.4 presents the results obtained at scale $L_1$ ($L_1$ = 1500 nm, i.e. cubic grid blocks of
size 1500 × 1500 × 1500 nm³).

<figure><img src="image_5.png" /><figcaption>**Fig. 4** - Distribution of the dimensionless components of CCL oxygen effective tensor $D_{O_2eff-L_1}$: (a) for grid block size $L_1$ =1500 nm without Knudsen diffusion, (b) for grid block size $L_1$ =1500 nm with Knudsen diffusion, (c) same as in (a) but for components of tensor $D^*_{O_2eff-L_1}$, same as in (b) but for components of tensor $D^*_{O_2eff-L_1}$.</figcaption></figure>

Each tensor component is expressed in dimensionless form in term of tortuosity according to classical relationships of the form $\frac{D_{O_2effxx-L_1}}{D_{O_2bulk}} = \frac{1}{\tau}$ (as introduced in Section 3.1, dropping for simplicity the subscript xx for τ, i.e. writing τ instead of $\tau_{xx}$). We have also plotted in Fig.4 the components of the tensor $D^*_{O_2eff-L_1} = \varepsilon D_{O_2eff-L_1}$ since this is the effective diffusion tensor often considered in the literature as discussed in Section 3.3. In dimensionless form, these components correspond to the plots of $\frac{\varepsilon}{\tau}$ in Fig.4 (in this figure, the notation “τ” thus
refers to a component tortuosity, as clarified in the Figure legend).

As can be seen, the porosity varies from one block to the other in the range [0.49 -0.58]. The diagonal components slightly increase with the porosity. The spreading of the results for a given intermediate porosity is quite moderate, especially when the Knudsen effect is considered. Consistent with the results on the porosity discussed in [26], blocks of the considered size can therefore be considered as representative elementary volume (REV). The effective diffusivity tensor non-diagonal components are smaller than the diagonal components but not completely negligible. The values of the tensor diagonal components for a given porosity depends on the direction for the considered block size, indicating an anisotropy effect.

The through-plane direction corresponds to the y axis. Somewhat surprisingly, the through-plane component $D_{O_2effyy-L_1}$ is in between the two components $D_{O_2effxx-L_1}$ and $D_{O_2effzz-L_1}$ in the in plane-directions. Since the milling direction in the FIB-SEM procedure is the z-direction, the question arises to as whether the difference between the two in-plane directions is due to the impact of the milling. On the other hand, a similar FIB-SEM procedure was applied to a micro-porous layer (MPL) as reported in [21] and no significant difference between the in-plane components of the diffusion tensor was obtained. Nevertheless, it would be interesting to further investigate the possible impact of the milling direction.

<figure><img src="image_6.png" /><figcaption>Fig.5 -Values of the oxygen transport dimensionless effective diffusion tensor diagonal component in through-plane direction for the various grid blocks of size $L_1$.</figcaption></figure>

Fig. 4 illustrates the significant impact of Knudsen diffusion. As pointed out above, the spreading of the results for a given porosity is less. The values of the tensor diagonal
components are significantly smaller compared to the case without consideration of the Knudsen diffusion. This is in line with previous works showing that indeed it is important to consider the Knudsen diffusion in the CCL, e.g. [9, 23, 62]. Fig.5 provides an illustration of the variability of the oxygen effective diffusion tensor over the various blocks of scale $L_1$.

In the fuel cells literature, the tortuosity $\tau$ is often estimated using Bruggeman's relationship, $\tau = \epsilon^{-0.5}$ [63]. As shown in Fig.6, Bruggeman like expressions can be used to represent our simulation results.

<figure><img src="image_7.png" /><figcaption>Fig.6 - Fitting of the evolution of the dimensionless oxygen effective diffusion tensor $D_{O_2eff-L_1}$ diagonal components as a function of porosity with the analytical model of Bruggeman.</figcaption></figure>

However, consistent with [63], it can be seen that Bruggeman relationship with the exponent of -0.5 overestimates the dimensionless effective diffusion coefficients. The fitting leads to significantly smaller exponents, depending on both the considered direction and the consideration or not of Knudsen diffusion. As pointed out in [21], the exponent actually depends on the considered pressure, temperature and gas when Knudsen diffusion is considered.

As mentioned before, the full effective tensor is not considered in most of previous works, which generally focus only on the characterization of the through-plane component, i.e. $D_{O_2effyy}$ or $D^*_{O_2effyy}$ with our notations and considered coordinate system. The comparison
with previous works presented in Fig.7 is therefore performed considering only $D_{O_2effyy}$ or
more precisely $D_{O_2effyy}^*$.

<figure><img src="image_8.png" /><figcaption>Fig.7 - Comparison of through-plane effective diffusivity obtained in the present work with other pore-scale simulation results and experimental data.</figcaption></figure>

The results of our simulations in Fig.7 are for the scale $L_1$ (cubic blocks with 1500 nm sides).
As can be seen, our results when Knudsen diffusion is considered agree with those reported
in He et al. [9] and Lange et al. [23] based on simulations on synthetic images generated
according to the same generation technique, but not with the results reported in Chen et al.
[66], where a different image generation method was used. As discussed in [9], this is an
indication that the details of the image generation are important. Somewhat surprisingly, the
agreement with the results of Chen et al. [66] is much better when Knudsen diffusion is
ignored whereas the agreement with the results reported in Lange et al. [23] become much
less good under this condition. It is likely that this results from differences in the geometry
and morphology between our FIB-SEM images and the synthetic images used by Lange et al.
[23]. Knudsen diffusion is more affected by the smaller pore sizes than Fick diffusion. In
other words, the impact of the various classes of pores on the $O_2$ transport is different
depending on whether the Knudsen diffusion is considered or not. This could explain why
the quality of the agreement with previous results depend on the consideration or not of
Knudsen diffusion.

The comparison with experimental results in Fig.7 is not obvious since the agreement
between two series of experimental measurements by the same group, i.e. Yu et al. [64, 65] is poor. As discussed in [9], this is notably due to the lack of accuracy in the determination of the porosity (the measurement technique used in [64] is different from the one used in [65]). The fact that our results are similar to those reported in He et al. [9] and Lange et al. [23] when Knudsen diffusion is considered can be also seen as an illustration that the protrusion of the ionomer filaments in the pore space has a limited impact on O2 transport since the ionomer filaments are not considered in [9] or [23]. This is consistent with the quite low volume fraction of the ionomer filaments (≈ 1%). As in [9], [23] or [66], our simulations are performed considering a dry pore space. As noted in [9] and consistent with the simulation results presented in [11], the presence of some liquid water in the secondary pores in the experiments might explain in part the overestimation of the effective oxygen through-plane diffusion coefficient compared to the experimental results reported in [6].

Tables 2 and 3 summarizes the result obtained at the scale of the whole CCL FIB-SEM image ($L_2$ scale). Consistent with the results at scale $L_1$, the off-diagonal components can be considered as negligible. An anisotropy effect is still present. As for the results at scale $L_1$, the diagonal component corresponding to the through plane diffusion lies in between the diagonal components in the in-plane directions, $D_{O_2effxx-L_2} < D_{O_2effyy-L_2} < D_{O_2effzz-L_2}$. As for the intermediate scale $L_1$, the impact of Knudsen diffusion is significant with a reduction of the diagonal component by a factor 0.4 compared to the simulations not considering the Knudsen diffusion.

<table><caption>**Table 2** - Dimensionless components of CCL oxygen effective diffusion tensor $D_{O_2eff}$, i.e. $(\frac{1}{\tau})$, at the CCL scale ($L_2$). The through plane coefficient is $D_{O_2effyy-L_2}$. $D_{bulk}$ is a shorthand for $D_{O_2bulk}$ and subscript $O_2$ is omitted for brevity in the tensor components.</caption><thead><tr><th></th><th>D<sub>effxx-L<sub>2</sub></sub>/D<sub>bulk</sub></th><th>D<sub>effxy-L<sub>2</sub></sub>/D<sub>bulk</sub></th><th>D<sub>effxz-L<sub>2</sub></sub>/D<sub>bulk</sub></th><th>D<sub>effyx-L<sub>2</sub></sub>/D<sub>bulk</sub></th><th>D<sub>effyy-L<sub>2</sub></sub>/D<sub>bulk</sub></th><th>D<sub>effyz-L<sub>2</sub></sub>/D<sub>bulk</sub></th><th>D<sub>effzx-L<sub>2</sub></sub>/D<sub>bulk</sub></th><th>D<sub>effzy-L<sub>2</sub></sub>/D<sub>bulk</sub></th><th>D<sub>effzz-L<sub>2</sub></sub>/D<sub>bulk</sub></th></tr></thead><tbody><tr><td>With Knudsen diffusion</td><td>0.14</td><td>0.0001</td><td>0.002</td><td>-0.001</td><td>0.17</td><td>0.05</td><td>0.0001</td><td>0.048</td><td>0.24</td></tr><tr><td>Without Knudsen diffusion</td><td>0.363</td><td>0.004</td><td>0.003</td><td>0.0026</td><td>0.457</td><td>0.11</td><td>0.0017</td><td>0.10</td><td>0.578</td></tr></tbody></table>
<table><caption>**Table 3** - Dimensionless components of CCL oxygen effective diffusion tensor $D_{O_2eff}^*$, i.e. $(\frac{\epsilon}{\tau})$, at the CCL scale ($L_2$). The through plane coefficient is $D_{effyy-L_2}$. $D_{bulk}$ is a shorthand for $D_{O_2bulk}$ and subscript $O_2$ is omitted for brevity in the tensor components.</caption><thead><tr><th></th><th>εD<sub>effxx</sub>-L<sub>2</sub>/D<sub>bulk</sub></th><th>εD<sub>effxy</sub>-L<sub>2</sub>/D<sub>bulk</sub></th><th>εD<sub>effxz</sub>-L<sub>2</sub>/D<sub>bulk</sub></th><th>εD<sub>effyx</sub>-L<sub>2</sub>/D<sub>bulk</sub></th><th>εD<sub>effyy</sub>-L<sub>2</sub>/D<sub>bulk</sub></th><th>εD<sub>effyz</sub>-L<sub>2</sub>/D<sub>bulk</sub></th><th>εD<sub>effzx</sub>-L<sub>2</sub>/D<sub>bulk</sub></th><th>εD<sub>effzy</sub>-L<sub>2</sub>/D<sub>bulk</sub></th><th>εD<sub>effzz</sub>-L<sub>2</sub>/D<sub>bulk</sub></th></tr></thead><tbody><tr><td>With Knudsen diffusion</td><td>0.08</td><td>0.00005</td><td>0.0011</td><td>-0.0005</td><td>0.093</td><td>0.028</td><td>0.00005</td><td>0.026</td><td>0.132</td></tr><tr><td>Without Knudsen diffusion</td><td>0.2</td><td>0.0022</td><td>0.0016</td><td>0.0014</td><td>0.25</td><td>0.06</td><td>0.0009</td><td>0.055</td><td>0.32</td></tr></tbody></table>

## 4.2 Proton effective conductivity (reconstructed ionomer)

As mentioned in section 3.5, a finer computational mesh is needed for computing the effective proton conductivity tensor (Table 1). This mesh is obtained by subdividing the FIB-SEM voxels by a factor of 6 in each direction. As indicated in Table 1, this leads to computational cubic finite volumes of side 0.833 nm for the proton transport computation. Due to the finer discretization, the grid block size $L_1$ is $250 \times 250 \times 250$ nm³, thus smaller than that for the oxygen transport computation. This corresponds to grid block of $300 \times 300 \times 300$ computational voxels. These blocks are partitioning a CCL sub-image of $1500 \times 1500 \times 1500$ nm³. The considered sub-image of $1500 \times 1500 \times 1500$ nm³ was selected based on the results presented in Fig.7 and Table 3 so as to be representative of the values reported in Table 3 as regards the oxygen effective diffusion tensor. Thus, the results presented in what follows for the scale $L_1$ are over 216 blocks of size $250 \times 250 \times 250$ nm³ partitioning the selected cubic image of size 1500 nm. The results are presented in dimensionless form in terms of the ratio $\frac{\epsilon_{I_R}}{\tau}$ where $\epsilon_{I_R}$ is the reconstructed ionomer volume fraction in the considered domain and $\tau$ is the tensor component tortuosity for easier comparisons with previous results from literature.
<figure><img src="image_9.png" /><figcaption>Fig.8 - Distribution of the dimensionless component of CCL proton effective conductivity tensor $\sigma_{\text{Peff}^{-I_R}}^*$ at scale $L_1$: a) diagonal components, b) diagonal components with off-diagonal components $\sigma_{\text{Peff}y^{-I_R}}^*$ and $\sigma_{\text{Peff}y^{-I_R}}^*$.</figcaption></figure>

Fig.8 presents the results obtained at scale $L_1$ ($L_1$ = 250 nm, i.e. cubic grid blocks of size 250 × 250 × 250 nm³). As can be seen, the variability of the tensor over the 216 blocks is significant, more than a factor of 10 between the less conductive and the most conductive blocks. The CCL is therefore significantly heterogeneous at this scale as regards the proton conductivity. Globally, the diagonal components of the proton conductivity tensor increase, as it could be expected, with the increasing reconstructed ionomer volume fraction ($\epsilon_{I_R}$). As for the oxygen transport, the through plane diagonal component ($\sigma_{\text{Peff}yy^{-I_R}}^*$) tends to lie in between the diagonal components in the in-plane directions with the same trend as regards these components, namely $\sigma_{\text{Peff}xx^{-I_R}}^* < \sigma_{\text{Peff}yy^{-I_R}}^* < \sigma_{\text{Peff}zz^{-I_R}}^*$. This will be more evident with the results at scale $L_2$ presented further below. A noticeable difference compared to the oxygen transport is the fact that some off-diagonal components, namely $\sigma_{\text{Peff}yz^{-I_R}}^*$ and $\sigma_{\text{Peff}zy^{-I_R}}^*$, are not negligible compared to the diagonal components. These off-diagonal components are compared with the diagonal components in Fig.8b. The other off-diagonal components are much smaller and therefore are not shown in Fig.8b. In other words, the CCL is anisotropic at this scale and the main axes of the proton conductivity tensor do not correspond to the cartesian coordinate systems associated with the FIB-SEM image.
<figure><img src="image_10.png" /><figcaption>a)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>b)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>c)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>d)</figcaption></figure>

**Fig.9 -** Distribution of the diagonal dimensionless components of CCL proton effective conductivity tensor $\sigma^*_{Peff-I_R}$ at scale $L_I$ (blocks of size 250 × 250 × 250 nm³): a) as a function of the coverage, b) as a function of the $I_R/C$ ratio, c) as a function of block porosity, d) illustration of the variability of the through-plane diagonal component $\frac{\sigma^*_{Peffyy-I_R}}{\sigma_P}$ over the 216 blocks.

The coverage, i.e. the fraction of the interfacial surface between the secondary pores and the solid phase corresponding to ionomer, varies from one block to the other as indicated in Fig.9a showing the variation of the diagonal components of the proton conductivity tensor as a function of coverage. The trend is, as it could be expected, the increase of the tensor diagonal components with an increasing coverage. The variability of the diagonal components is also illustrated in Fig.9 as a function of the $I_R/C$ ratio (Fig.9b) and the block porosity (Fig.9c). The $I_R/C$ ratio considered here correspond to the apparent I/C ratio considered in [26]. It is defined as follows. At the end of the ionomer reconstruction procedure, the segmented solid phase in FIB-SEM images is partitioned in reconstructed
ionomer voxels and complementary solid phase voxels. Thus,

$$V_{FIB-SEM-SSP} = V_{IR} + V_{CSP} \qquad (28)$$

where $V_{FIB-SEM-SSP}$ is the volume of the FIB-SEM segmented solid phase in the considered image, $V_{IR}$ is the volume of reconstructed ionomer in the image, i.e. the filament and coating ionomer, and $V_{CSP}$ is the volume of the FIB-SEM segmented solid phase not identified as filament or coating ionomer, i.e. the complementary solid phase. Then the apparent I/C ratio, i.e. $I_R/C$, is defined as $\frac{I_R}{C} = \frac{\rho_I V_{IR}}{\rho_c V_{CSP}}$, where $\rho_I$ is the ionomer density ($\rho_I = 2 \text{ g/cm}^3$) and $\rho_c$ is the carbon density ($\rho_c = 1.8 \text{ g/cm}^3$). This ratio is defined as apparent in [26] since, as explained previously, some ionomer, i.e. the inter-carbon grain ionomer, is located elsewhere in the segmented solid phase than in the coating or the filaments.

Fig.9c shows that the proton conductivity decreases with the increasing block porosity, consistent with Fig.8a, 9a and 9b since the ionomer volume fraction, the coverage and the $I_R/C$ ratio increases with a decreasing block porosity. The spatial heterogeneity of the through plane diagonal component ($\sigma_{Peffyy-I_R}^*$) is illustrated in Fig.8d. This figure suggests some layering (illustrated by the blue blocks on the right face), i.e. some possible spatial correlation between the blocks as regards the distribution of the proton conductivity. However, this aspect has not been investigated further and is left for a future work.

The results at the larger scale, i.e. $L_2$ (domain of $1500 \times 1500 \times 1500 \text{ nm}^3$ corresponding to the block assembly depicted in Fig. 9d), are summarized in Tables 4 and 5.

<table><caption>**Table 4** - Dimensionless components of CCL proton effective conductivity tensor $\sigma_{Peff-I_R}$, i.e. $(\frac{1}{\tau})$ at scale $L_2$. The through plane coefficient is $\sigma_{Peffyy-I_R-L_2}$ (for brevity the subscript $I_R$ is omitted in the table for the tensor components).</caption><thead><tr><th>$\frac{\sigma_{Peffxx-L_2}}{\sigma_p}$</th><th>$\frac{\sigma_{Peffxy-L_2}}{\sigma_p}$</th><th>$\frac{\sigma_{Peffxz-L_2}}{\sigma_p}$</th><th>$\frac{\sigma_{Peffyx-L_2}}{\sigma_p}$</th><th>$\frac{\sigma_{Peffyy-L_2}}{\sigma_p}$</th><th>$\frac{\sigma_{Peffyz-L_2}}{\sigma_p}$</th><th>$\frac{\sigma_{Peffzx-L_2}}{\sigma_p}$</th><th>$\frac{\sigma_{Peffzy-L_2}}{\sigma_p}$</th><th>$\frac{\sigma_{Peffzz-L_2}}{\sigma_p}$</th></tr></thead><tbody><tr><td>0.14</td><td>0.001</td><td>-0.004</td><td>-0.002</td><td>0.2</td><td>0.112</td><td>-0.004</td><td>0.116</td><td>0.23</td></tr></tbody></table>
<table><caption>**Table 5** - Dimensionless components of CCL proton effective conductivity tensor $\sigma_{\text{Peff}^{-I_R}}^*$, i.e. $(\frac{\epsilon_{I_R}}{\tau})$ at scale $L_2$. The through plane coefficient is $\sigma_{\text{Peffyy}^{-I_R-L_2}}^*$ (for brevity the subscript $I_R$ is omitted in the table for the tensor components)</caption><thead><tr><th>ε<sub>I_R</sub>σ<sub>Peffxx-L<sub>2</sub></sub>/σ<sub>p</sub></th><th>ε<sub>I_R</sub>σ<sub>Peffxy-L<sub>2</sub></sub>/σ<sub>p</sub></th><th>ε<sub>I_R</sub>σ<sub>Peffxz-L<sub>2</sub></sub>/σ<sub>p</sub></th><th>ε<sub>I_R</sub>σ<sub>Peffyx-L<sub>2</sub></sub>/σ<sub>p</sub></th><th>ε<sub>I_R</sub>σ<sub>Peffyy-L<sub>2</sub></sub>/σ<sub>p</sub></th><th>ε<sub>I_R</sub>σ<sub>Peffyz-L<sub>2</sub></sub>/σ<sub>p</sub></th><th>ε<sub>I_R</sub>σ<sub>Peffzx-L<sub>2</sub></sub>/σ<sub>p</sub></th><th>ε<sub>I_R</sub>σ<sub>Peffzy-L<sub>2</sub></sub>/σ<sub>p</sub></th><th>ε<sub>I_R</sub>σ<sub>Peffzz-L<sub>2</sub></sub>/σ<sub>p</sub></th></tr></thead><tbody><tr><td>0.017</td><td>0.00012</td><td>-0.0005</td><td>-0.00025</td><td>0.0248</td><td>0.0139</td><td>-0.0005</td><td>0.014</td><td>0.029</td></tr></tbody></table>

The results presented in Tables 4 and 5 confirm the trends discussed above for the scale $L_1$, i.e. $\sigma_{\text{Peffxx}-I_R-L_2} < \sigma_{\text{Peffyy}-I_R-L_2} < \sigma_{\text{Peffzz}-I_R-L_2}$, the off-diagonal components $\sigma_{\text{Peffzy}-I_R-L_2}$ and $\sigma_{\text{Peffyz}-I_R-L_2}$ are not negligible compared to the diagonal components contrary to the other off-diagonal components. From Table 4, it can be seen that the proton effective conductivity tensor is not symmetrical. However, the deviation to symmetry is weak.

Compared to the oxygen effective diffusion tensor (Tables 2 and 3), it can be seen that the off-diagonal components are thus less small compared to the diagonal components. In this respect, it can be noted that the volume fraction of the ionomer in the images (on the order of 0.1 on average) is significantly smaller than the secondary pore volume fraction (on the order of 0.5). The ionomer can therefore be expected to be more tortuous and less well connected than the secondary pores. This might explain the observed differences between the two tensors. Further investigations are however needed to clarify this point in more details.

A comparison with results from the literature [67-73] is proposed in Fig. 10 considering the results obtained at scale $L_1$, which is comparable to the size of the images generated in [3] for instance.

As can be seen, the experimental results from literature, which are the same as those considered in [3], spread over a quite significant range of values for a given ionomer volume fraction, except for the lower ionomer volume fractions for which the few available results are more grouped. Our results agree with the experimental measurement reported by Havranek et al. [69] and are close to the experimental data of Suzuki and al. [70] and Modestov et al. [71]. They are not too far from the experimental measurements of Du et al. [71]. However, the deviation is noticeable with the measurements of Yakolev et al. [68], Liu et al. [73] and Boyer et al. [67]. Several aspects are discussed in [3] to explain such a
variability between the results from different groups. These notably include the value of the ionomer conductivity ($\sigma_p$ in Eq.(9)) which can be different in thin ionomer films from the value for thicker films generally considered to obtain the factor ($\frac{\epsilon_I}{\tau}$) from the experimental measurement and the impact of liquid water. In this respect, we recall that the FIB-SEM image corresponds to a dry CCL.

<figure><img src="image_14.png" /><figcaption>Fig.10 - Comparison of the evolution of the dimensionless through plane effective proton conductivity obtained in the present work ($\frac{\sigma_{\text{Peff}}^{\text{I Ry-L1}}}{\sigma_p}$, blue dots) with other pore-scale simulations and experimental data from literature. See text for the meaning of correction 1 and correction 2.</figcaption></figure>

The possible impact of the coverage is also mentioned in [3]. However, in the case of our ionomer reconstruction procedure, the coverage is strongly correlated to $\epsilon_{IR}$ (as indicated by Figs.8a and 9a). Thus, a large variability of the coverage is not expected for a given $\epsilon_{IR}$. The likely impact of liquid water is also evoked in [9] where it is noted that the simulations on synthetic images tend to underestimate the proton conductivity compared to measurements. The results from simulations of He et al. [9] are reported in Fig.10. As noted in [9], they are
representative of the results obtained by simulations on synthetic images in other works [23, 62]. Interestingly, our results are in better agreement with the experimental data than those of [9], which underestimate the proton conductivity. This might be an indication than the coating ionomer is better represented in our approach. In this respect, it can be noted that in most simulations on synthetic images, only coating ionomer is considered. The ionomer possibly present in the solid phase away from the secondary pores / solid phase interface, i.e. the inter-carbon grain ionomer, is not considered. However, as mentioned before and detailed in [26], this fraction of ionomer can be quite significant. As reported in [26] for the coverage of 0.8 considered in the present work, 7% of the expected mass or volume of the ionomer is in the filaments, 50 % is in the ionomer coating and 41% is present as inter-carbon grain ionomer (for the whole image).

The sensitivity of the results to inter-carbon grain ionomer and liquid water (inside the primary pores) is explored in what follows.

### Sensitivity to proton transport in complementary solid phase.

As discussed previously, it is considered that a significant fraction of the ionomer is present in the solid phase complementary to the reconstructed ionomer (CSP in Fig. 2). This phase is considered to be a nanoporous medium made of carbon grains, platinum, ionomer and voids (primary pores possibly containing water). The proton transport governing equation in the CSP is similar to Eq. 4. It can be upscaled following the developments in Sec. 3 to obtain a macro-scale equation valid in the complementary solid phase, which is expressed as

$$ \nabla \cdot (\varepsilon_{ICSP} \sigma_{Peff-CSP} \nabla \varphi_I) = 0 \qquad (29) $$

where $\varepsilon_{ICSP}$ is the volume fraction of ionomer in the complementary solid phase and $\sigma_{Peff-CSP}$ is the proton effective conductivity in the complementary solid phase. The latter is expressed as $\sigma_{Peff-CSP} = \frac{\sigma_P}{\tau_{CSP}}$, where $\tau_{CSP}$ is the tortuosity in the complementary solid phase. For simplicity a homogeneous and isotropic conductivity tensor is assumed in the complementary solid phase, which is consistent with the somehow random distribution of the carbon grains. Eq.(29) is solved together with Eq.(9) in the reconstructed ionomer, namely

$$ \nabla \cdot (\sigma_P \nabla \varphi_I) = 0 \qquad (30) $$

Because of the continuity of fluxes and potential at the interface $I_R-CSP$, it is possible to solve a single problem over the computational domain illustrated in Fig.11, using the following equation

$$ \nabla \cdot (\Sigma_P \nabla \varphi_I) = 0 \qquad (31) $$
with $\Sigma_P = \sigma_P$ in the reconstructed ionomer voxel, $\Sigma_P = \varepsilon_{ICSP} \frac{\sigma_P}{\tau_{CSP}}$ in the complementary solid phase, and $\Sigma_P = 10^{-5}\sigma_P$ in the secondary pore space (penalization technique).

<figure><img src="image_15.png" /><figcaption>Fig.11 - Local proton conductivity distribution in the representative grid block ($\varepsilon = 0.52, \varepsilon_I = 0.126, \text{coverage} = 0.82$). Coating ionomer and filament in red. Secondary pore space in darker blue (within the framework of the used penalization technique a very low conductivity value is assigned to secondary pore voxels), complementary solid phase in lighter blue.</figcaption></figure>

The simulations are performed over a representative grid block selected among the 216 grid blocks considered to obtain the results plotted in Fig.10. The selected block is considered representative because its porosity and coverage are close to the values for the whole image. The volume fraction $\varepsilon_{ICSP}$ in Eq. (29) is determined using Eq. (A.21) (derived in the Appendix).

For the considered block, $\varepsilon_{ICSP} = 0.1$. The tortuosity $\tau_{CP}$ in the complementary solid phase is varied and the effective through plane proton conductivity is obtained as described in Section 3 with a ionomer total volume fraction $\varepsilon_I = \varepsilon_{IR} + \varepsilon_{ICG}$ where $\varepsilon_{IR}$ is the reconstructed ionomer volume fraction in the image as before and $\varepsilon_{ICG}$ is the volume fraction of inter-
carbon grain ionomer in the image ($\epsilon_{ICG} = \epsilon_{ICSP}(1 - \epsilon_{IR} - \epsilon) = 0.035$). The obtained results are shown in Fig. 12.

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

Fig.12 - Variation of dimensionless through-plane proton effective conductivity $\frac{\sigma_{effyy}^*}{\sigma_P}$ as a function of tortuosity $\tau_{CSP}$ in the complementary solid phase when proton transport in the ionomer in the complementary solid phase (inter-carbon grain ionomer) is considered (red square symbols) and when the proton transport is considered in both the inter-carbon grain ionomer and liquid water in the primary pores (blue triangle symbols). The black straight line labelled “initial value” corresponds to the case when only the reconstructed ionomer (filaments and secondary pore wall coating) is considered.

As can be seen from Fig.12 (square red symbols), the impact of the inter-carbon grain ionomer is noticeable for all the values of tortuosity tested. This tortuosity highly depends on the morphology of the ionomer network in the complementary solid phase. For a tortuosity comparable to the tortuosity of the ionomer coating (5 from Table 4) for instance, the increase in the effective conductivity shown in Fig.12 is by a factor of 1.75. To obtain the
results shown in Fig.12, we have implicitly assumed that the ionomer network in the complementary solid phase is percolating. However, it is argued in [74] that this network does not form a percolating network in the case of Nafion. On the other hand, it is considered in [44] that the primary pores are filled with water under the usual PEMFC operating conditions. This water can provide a pathway for the proton and the system formed by the ionomer and the water in the complementary solid phase can then percolate. Thus, let's assume that the primary pores are filled with water and contribute to the proton transport. To explore this case, we solve again Eq. (29) replacing $\varepsilon_{ICSP}$ in Eq. (27) by $\varepsilon_{IWC} = \varepsilon_{ICSP} + \varepsilon_{WC}$ where $\varepsilon_{WC}$ is the volume fraction of primary pores in the CSP. Varying again the tortuosity leads to the results shown with blue symbols in Fig.12. As it could be expected, the consideration of the transport in the water leads to a significant further increase in the proton effective conductivity. It can be noted that we have assumed for these simulations that the proton conductivity in the water was the same as in the ionomer, which perhaps leads to overestimate the impact of the water transport.

In order to further illustrate the impact of the proton transport in the complementary solid phase, the tortuosity is estimated using the classical Bruggeman relationship [63], namely $\tau_{CP} = \varepsilon_{IWC}^{-0.5}$, with $\varepsilon_{IWC} = \varepsilon_{ICSP}$ when only the inter-carbon grain ionomer is considered and $\varepsilon_{IWC} = \varepsilon_{ICSP} + \varepsilon_{WC}$ when water is considered with the inter-carbon grain ionomer. For the chosen representative grid block, this leads to the results reported in Fig.10.

In this figure, “correction 1” refers to the case where only the inter-carbon grain ionomer is considered whereas “correction 2” refers to the case where both the inter-carbon grain ionomer and water are considered. As can be seen, the consideration of the inter-carbon grain ionomer alone does not dramatically change the comparison with the experimental results. This is different when transport in water in the primary pores is considered with a significant shift toward a larger proton conductivity in line with the consideration in the literature [3, 9] on the likely impact of liquid water on the proton transport in the CCL. As also mentioned in [9], proton transport could also occur along the carbon or ionomer surfaces and thin liquid water films in the interfacial region with the secondary pores in addition to the transport in the primary pores [75, 76]. Although the results shown in Fig.10 are encouraging, they also illustrate that predicting the proton conductivity from simulations on a CCL microstructure and/or modelling is still a challenge, as also pointed out in [3, 9].

## 4.3 Impact of oxygen transport in ionomer

From the results obtained for the proton conductivity, i.e. the tortuosity of the reconstructed
ionomer (Table 4), one can compare the oxygen transport in the gas phase and the ionomer so as confirm that the oxygen transport in the ionomer is negligible compared to the oxygen transport by diffusion in the gas phase. The oxygen transport in the ionomer is also modelled using a diffusion equation [13], namely,

$$ \nabla \cdot (\mathbf{D}_{IO_{2eff}}^* \cdot \nabla C_{IO_2}) = 0 \qquad (32) $$

where $\mathbf{D}_{IO_{2eff}}^*$ is the oxygen effective diffusion tensor in the ionomer and $C_{IO_2}$ is the oxygen concentration in the ionomer. Boundary conditions apply at the interface V-I$_R$ (the inter-carbon grain ionomer is not considered in this sub-section), namely Eq. (2), with an equation corresponding to the continuity of fluxes. Considering only the through-plane transport and assuming that the oxygen transport in the gas phase occurs in parallel with the oxygen transport in the reconstructed ionomer leads to define the global oxygen effective diffusion coefficient as

$$ \frac{D_{G-O_2 effyy}^*}{D_{O_2-bulk}} = \frac{D_{O_2 effyy}^*}{D_{O_2-bulk}} + RTH_e \frac{D_{O_2 effyy}^*}{D_{O_2-bulk}} = \frac{\epsilon}{\tau_{O_2yy}} + RTH_e \frac{D_{IO_2-bulk}}{D_{O_2-bulk}} \frac{\epsilon_{I_R}}{\tau_{Pyy}} \qquad (33) $$

Eq.(33) can be also expressed as

$$ \frac{D_{G-O_2 effyy}^*}{D_{O_2 effyy}^*} = 1 + RTH_e \frac{D_{IO_2-bulk}}{D_{O_2-bulk}} \frac{\tau_{O_2yy}}{\tau_{Pyy}} \frac{\epsilon_{I_R}}{\epsilon} \qquad (34) $$

For representative values of the literature ($H_e = 5 \times 10^{-5}$ mol Pa$^{-1}$ m$^3$ [76]), $\frac{D_{IO_2-bulk}}{D_{O_2-bulk}} \approx 10^{-5}$ [13, 77]) and the data on the ratios $\frac{\epsilon_{I_R}}{\tau_{Pyy}}$ and $\frac{\epsilon}{\tau_{O_2yy}}$ indicated in Tables 3 and 5, it is found that the second term in the r.h.s. of Eq.(34), i.e. $RTH_e \frac{D_{IO_2-bulk}}{D_{O_2-bulk}} \frac{\tau_{O_2yy}}{\tau_{Pyy}} \frac{\epsilon_{I_R}}{\epsilon}$, is on the order of $10^{-7}$.

This confirms that the through plane transport of the oxygen in the ionomer is negligible compared to the through plane oxygen transport in the pore space.

It can be noted here that the above is based on the use of Henry's law to relate the oxygen concentration in the ionomer and in the gas. This is only a useful approximation, sufficient in order to compare the through plane oxygen transport in the ionomer and the gas phase. More refined descriptions have been studied in the literature to relate both concentrations, e.g. [78,79].

## 5. Conclusion.

Characterizing the effective transport properties of the various porous layers of Polymer Electrolyte Membrane Fuel Cells (PEMFC) is important as inputs for the models aiming at optimizing the PEMFC performances and analysing aging issues. This is particularly
important as regards the cathode catalyst layer (CCL), since the CCL is the locus of the oxygen reduction reaction. In this work, the effective oxygen diffusion in the secondary pores of a CCL and the effective proton conductivity were studied from direct numerical simulations on three-dimensional images of the CCL microstructure. The latter were obtained by combining FIB-SEM imaging and numerical methods for identifying the ionomer coating and filaments in the images. Contrary to most of previous works, which only focus on the through plane oxygen effective diffusion and effective proton conductivity components, the full tensors were characterized in this work. The results show that the obtained tensors are not isotropic. It is also found that some off-diagonal components are not negligible in the case of the proton effective conductivity tensor, thus indicating that the principal axis are not colinear to the image reference frame.

Regarding the oxygen effective diffusion through-plane component, our results are comparable to those obtained in previous works on synthetic images. They confirm that Knudsen diffusion has a significant impact and must be considered together with Fick diffusion.

The results show that the effective proton conductivity is heterogeneous at the scale of the computational domains typically considered in many previous works due to the variability in the ionomer distribution at this scale. Contrary to previous works based on numerical simulations on synthetic images, which tend to underestimate the proton effective conductivity through-plane component, our results are globally in better agreement with experimental data. This might be due to the better representativity of the CCL microstructure obtained via the FIB-SEM imaging compared to synthetic images. This is also seen as an indicator of the relevance of our ionomer reconstruction procedure. Nevertheless, the spread of the available experimental data is broad, which makes difficult to reach a firm conclusion regarding our method of reconstruction of the ionomer coating in the 3D images. This ionomer reconstruction method leads to distinguish the ionomer filaments, the ionomer coating and the inter-carbon grain ionomer, which corresponds to ionomer present in the complementary solid phase (CSP), i.e. the solid phase not identified in the reconstruction procedure as ionomer coating or ionomer filaments. The ionomer in the CSP should correspond to the inter-agglomerate ionomer in [3] and/or the free ionomer network discussed in [72]. The contribution of this fraction of ionomer was evaluated considering the proton transport in the CSP in addition to the transport in the ionomer coating and the filaments. The simulations indicate that this transport has an impact but not sufficiently for changing significantly the comparison with the experimental data. The impact of liquid water
in the primary pores on the proton transport in conjunction with the transport in the ionomer in the CSP was also explored. In agreement with discussions in previous works [3, 9], the simulations indicate a significant impact, that could potentially reconcile the simulations and the experimental data. In spite of the encouraging results obtained via the combination of FIB-SEM imaging and the ionomer reconstruction procedure presented in [26], all this shows that further work is still needed to better characterize the morphology of the ionomer in the CCL and the proton transport mechanisms.

**Appendix:** Inter-carbon grain ionomer volume fraction in complementary segmented solid phase.

This inter-carbon grain ionomer volume fraction is denoted by $\varepsilon_{ICSP}$. The complementary segmented solid phase (subscript CSP) is the segmented solid phase remaining after identification of coating ionomer and filament ionomer.

Let $V$ be the volume of the considered domain (volume of a block of size $250 \times 250 \times 250$ nm³). Then

$$V = V_{\text{void}} + V_{\text{Pt}} + V_{\text{C}} + V_{\text{I}} \quad (A-1)$$

Where $V_{\text{void}}$ is the volume of pores, $V_{\text{Pt}}$ is the volume of platinum, $V_{\text{C}}$ is the volume of carbon and $V_{\text{I}}$ is the volume of ionomer.

We also have

$$V_{\text{void}} = V_{\text{void-macro}} + V_{\text{void-micro}} \quad (A-2)$$

$V_{\text{void-macro}}$ is the volume of macropores (the one identified in the FIB-SEM image, also referred to as the secondary pores) and $V_{\text{void-micro}}$ is the volume of micro-meso pores (the pores in the Pt-carbon particles, also referred to as the primary pores).

$$V_{\text{I}} = V_{\text{IR}} + V_{\text{ICG}} \quad (A-3)$$

where $V_{\text{IR}}$ is the volume of reconstructed ionomer (the identified coating and filament ionomer) and $V_{\text{ICG}}$ is the volume of inter-carbon grain ionomer.

We have also the following relationships,

$$V_{\text{void-macro}} = \varepsilon V \quad (A-4)$$

Where $\varepsilon$ is the porosity identified from the FIB-SEM image ($\varepsilon \approx 0.55$ for the whole image).

$$V_{\text{SSP}} = V_{\text{C}} + V_{\text{void-micro}} + V_{\text{Pt}} + V_{\text{IR}} + V_{\text{ICG}} \quad (A-5)$$

Where $V_{\text{SSP}}$ is the volume of the segmented solid phase in the FIB-SEM image. For the considered block of size $250 \times 250 \times 250$ nm³ and porosity 0.55, $V_{\text{SSP}}$ and $V_{\text{IR}}$ are known.
We also known the expected I/C ratio (0.8), denoted by

$$R_{I/C} = \frac{m_I}{m_C} = \frac{\rho_I V_I}{\rho_C V_C} \qquad (A-6)$$

where $\rho_I$ is the ionomer density (2 g/cm³) and $\rho_C$ is the carbon density (1.8 g/cm³) and the platinum carbon ratio (0.462),

$$R_{Pt/C} = \frac{m_{Pt}}{m_{Pt}+m_{C}} = \frac{\rho_{Pt}V_{Pt}}{\rho_{Pt}V_{Pt}+\rho_{C}V_{C}} \qquad (A-7)$$

Where $\rho_{Pt}$ is the platinum density (21.45 g/cm³).

From Eqs.(A-6) and (A-7), the volume $V_I$ and $V_{Pt}$ can thus be expressed as a function of the carbon volume,

$$V_I = \frac{\rho_C}{\rho_I} R_{I/C} V_C \qquad (A-8)$$

$$V_{Pt} = \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} V_C \qquad (A-9)$$

By definition,

$$\varepsilon_{ICSP} = \frac{V_{IICG}}{V_C + V_{void-micro} + V_{IICG} + V_{Pt}} \qquad (A-10)$$

Which can be expressed as

$$\varepsilon_{ICSP} = \frac{V_I - V_{I_R}}{V_C + V_{void-micro} + V_I - V_{I_R} + V_{Pt}} \qquad (A-11)$$

Combining Eq.(A-11) with Eqs.(A-8) and (A-9) we obtain,

$$\varepsilon_{ICSP} = \frac{\frac{\rho_C}{\rho_I} R_{I/C} V_C - V_{I_R}}{V_C + V_{void-micro} + \frac{\rho_C R_{I/C}}{\rho_I} V_C - V_{I_R} + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} V_C} \qquad (A-12)$$

Which can be expressed as,

$$\varepsilon_{ICSP} = \frac{\frac{\rho_C}{\rho_I} R_{I/C} V_C - V_{I_R}}{V_C \left( 1 + \frac{\rho_C R_{I/C}}{\rho_I} + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} \right) + V_{void-micro} - V_{I_R}} \qquad (A-13)$$

The next step is to express $V_{void-micro}$ as a function of $V_C$

Let $\varepsilon_{micro}$ denotes the porosity of the carbon particles ($\varepsilon_{micro} \approx 0.25$). It is known that a fraction of the platinum particles is in the carbon particles micro-porosity. Here we assume that the platinum particles are taken into account in the known micro-porosity of the carbon particles (in other words, the microporosity would be larger than considered below if the platinum particles inside the carbon particles were removed). Then, we have

$$V_C + V_{Pt} = (1 - \varepsilon_{micro}) V_{Pt-carbon-particles} \qquad (A-14)$$

$$V_{void-micro} = \varepsilon_{micro} V_{Pt-carbon-particles} \qquad (A-15)$$

Where $V_{Pt-carbon-particles}$ is the volume of the Pt-carbon-particles. Combining Eqs. (A-14) and (A-15) yields,
$$V_{\text{void-micro}} = \frac{\varepsilon_{\text{micro}}}{(1-\varepsilon_{\text{micro}})} (V_C + V_{Pt}) \qquad (\text{A-16})$$

Which is combined with Eq.(A-9) to obtain

$$V_{\text{void-micro}} = \frac{\varepsilon_{\text{micro}}}{(1-\varepsilon_{\text{micro}})} \left( 1 + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} \right) V_C \qquad (\text{A-17})$$

Substituting this expression in Eq.(A-13) yields,

$$\varepsilon_{ICSP} = \frac{\frac{\rho_C}{\rho_I} R_{I/C} V_C - V_{I_R}}{V_C \left( 1 + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} + \frac{\varepsilon_{micro}}{(1-\varepsilon_{micro})} \left( 1 + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} \right) \right) - V_{I_R}} \qquad (\text{A-18})$$

The sole unknown in Eq.(A-18) is $V_C$. To obtain $V_C$ we start from Eq.(A-5) into which we substitute Eqs. (A-8), (A-9) and (A-17), ($V_{SSP} = V_C + V_{\text{void-micro}} + V_{Pt} + V_{I_R} + V_{I_{ICG}}$)

$$V_{SSP} = V_C \left[ 1 + \frac{\varepsilon_{micro}}{(1-\varepsilon_{micro})} \left( 1 + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} \right) + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} + \frac{\rho_C}{\rho_I} R_{I/C} \right] \qquad (\text{A-19})$$

Which leads to express $V_C$ as

$$V_C = \frac{V_{SSP}}{\left[ 1 + \frac{\varepsilon_{micro}}{(1-\varepsilon_{micro})} \left( 1 + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} \right) + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} + \frac{\rho_C}{\rho_I} R_{I/C} \right]} \qquad (\text{A-20})$$

Combining Eq.(A-18) and (A-20) leads finally to

$$\varepsilon_{ICSP} = \frac{\frac{\rho_C}{\rho_I} R_{I/C} V_C - V_{I_R}}{V_{SSP} - V_{I_R}} = \frac{\left( \frac{\rho_C}{\rho_I} R_{I/C} \right) \left( \frac{V_{SSP}}{1 + \frac{\varepsilon_{micro}}{(1-\varepsilon_{micro})} \left( 1 + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} \right) + \frac{\rho_C R_{Pt/C}}{\rho_{Pt}(1-R_{Pt/C})} + \frac{\rho_C}{\rho_I} R_{I/C} \right) - V_{I_R}}{V_{SSP} - V_{I_R}} \qquad (\text{A-21})$$

### CRediT authorship contribution statement

**Mohamed Ahmed-Maloum:** Investigation, Methodology, Software, Writing – Original draft. **Joel Pauchet:** Writing - Review & Editing. **Michel Quintard:** Conceptualization, Methodology, Supervision, Writing - Review & Editing. **Marc Prat:** Conceptualization, Methodology, Supervision, Writing - Review & Editing.

### 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:** This project has received funding from the Fuel Cells and Hydrogen 2 Joint Undertaking (now Clean Hydrogen Partnership) under Grant Agreement No 875025. This Joint Undertaking receives support from the European Union's Horizon 2020 Research
References

1. Karan K. PEFC catalyst layer: Recent advances in materials, microstructural characterization, and modeling. Current Opinion in Electrochemistry 2017; 1: 27-35.

2. Wang YD, Meyer Q, Tang K, McClure JE, White RT, Kelly ST, Crawford MM, Iacoviello F, Brett DJL, Shearing PR, Mostaghimi P, Zhao C, Armstrong RT. Large-scale physically accurate modelling of real proton exchange membrane fuel cell with deep learning. Nature Communications 2023; 14:745.

3. Olbrich W, Kadyk T, Sauter U, Eikerling M, Gostick J. Structure and conductivity of ionomer in PEM fuel cell catalyst layers: a model-based analysis. Scientific Reports 2023; 13:14127.

4. Djilali N. Computational modelling of polymer electrolyte membrane (PEM) fuel cells: Challenges and opportunities. Energy 2007; 32 (4): 269-280.

5. Wu HW. A review of recent development: Transport and performance modeling of PEM fuel cells. Applied energy 2016;165 (1): 81-106.

6. Tzelepis S, Kavadias KA, Marnellos GE, Xydis G. A review study on proton exchange membrane fuel cell electrochemical performance focusing on anode and cathode catalyst layer modelling at macroscopic level. Renewable and Sustainable Energy Reviews 2021; 151: 111543.

7. Chen Q, Niu Z, Li H, Jiao K, Wang Y. Recent progress of gas diffusion layer in proton exchange membrane fuel cell: Two-phase flow and material properties. International Journal of Hydrogen Energy 2021 46(12): 8640-8671.

8. Kosakian A, Wei F, Zhou J, Jung S, Sharman J, Secanell M. Impact of liquid-water accumulation and drainage cycles on fuel-cell performance and stability. Electrochemica Acta 2023; 469: 143221.

9. He Y, Bai M, Hao L. Pore-Scale Simulation of Effective Transport Coefficients in the Catalyst Layer of Proton Exchange Membrane Fuel Cells. Journal of The Electrochemical Society 2023; 170: 044501.

10. Chen L, Zhang R, He P, Kang Q, He YL, Tao WQ. Nanoscale simulation of local gas transport in catalyst layers of proton exchange membrane fuel cells. Journal of Power Sources 2018; 400: 114-125.

11. Fathi H, Raoof A, Mansouri SH, van Genuchten MT. Effects of Porosity and Water Saturation on the Effective Diffusivity of a Cathode Catalyst Layer. J. Electrochem. Soc. 2017; 164: F298.

12. Lange KJ, Sui PC, Djilali N. Pore Scale Simulation of Transport and Electrochemical Reactions in Reconstructed PEMFC Catalyst Layers. J. Electrochem. Soc. 2010; 157: B1434-B1442.

13. Barreiros Salvado M, Schott P, Guétaz L, Gerard M, David T, Bultel Y. Towards the understanding of transport limitations in a proton-exchange membrane fuel cell catalyst layer: Performing agglomerate scale direct numerical simulations on electron-microscopy-based geometries. Journal of Power Sources 2021; 482: 228893.

14. Jung S, Sabharwal M, Jarauta A, Wei F, Gostick J, Gingras M, Secanell M. Estimation of relative transport properties in porous transport layers using pore-scale and pore-network simulations. Journal of The Electrochemical Society 2021; 168: 064501.

15. Sabharwal M, Gostick JT, Secanell M. Virtual liquid water intrusion in fuel cell gas diffusion media. Journal of the Electrochemical Society 2018; 165 (7): F553.
16. Agaesse T, Lamibrac A, Buechi F, Pauchet J, Prat M. Validation of pore network simulations of ex-situ water distributions in a gas diffusion layer of proton exchange membrane fuel cells with X-ray tomographic images. Journal of Power Sources 2016; 331: 462-474.

17. Tranter TG, Gostick JT, Burns AD, Gale WF. Pore Network Modeling of Compressed Fuel Cell Components with OpenPNM. Fuel Cells 2016; 16: 504.

18. El Hannach M, Pauchet J, Prat M. Pore network modeling: Application to multiphase transport inside the cathode catalyst layer of proton exchange membrane fuel cell, Electrochimica Acta 2011; 56: 10796-10808.

19. Sadeghi MA, Khan ZA, Agnaou M, Hu L, Litster S, Kongkanand A, Padgett E, Muller DA, Friscic T, Gostick J. Predicting PEMFC performance from a volumetric image of catalyst layer structure using pore network modeling, Applied Energy 2024; 353:122004.

20. Wang CY. Fundamental Models for Fuel Cell Engineering, Chemical Reviews 2004; 104 (10): 4727-4766.

21. Ahmed-Maloum M, David T, Guetaz L, Duru P, Pauchet J, Quintard M, Prat M. Computation of oxygen diffusion properties of the gas diffusion medium -microporous layer assembly from the combination of X-ray microtomography and focused ion beam three dimensional digital images. Journal of Power Sources 2023; 561: 232735.

22. Sabharwal M, Pant LM, Patel N, Secanell M. Computational Analysis of Gas Transport in Fuel Cell Catalyst Layer under Dry and Partially Saturated Conditions. Journal of the Electrochemical Society 2019; 166: F3065-F3080.

23. Lange KJ, Sui PC, Djilali N. Pore Scale Simulation of Transport and Electrochemical Reactions in Reconstructed PEMFC Catalyst Layers. J. Electrochem. Soc. 2010; 157: B1434-B1442.

24. Kim SH, Pitsch H. Reconstruction and Effective Transport Properties of the Catalyst Layer in PEM Fuel Cells, J. Electrochem. Soc. 2009; 156: B673-B681.

25. Grunewald JB, Mistry AN, Verma A, Goswami N, Mukherjee PP, Fuller TF. Perspective—mesoscale physics in the catalyst layer of proton exchange membrane fuel cells. J. Electrochem. Soc. 2019; 166: F3089-F3092.

26. Ahmed-Maloum M, David T, Guetaz L, Morin A, Pauchet J, Quintard M, Prat M. Characterizing PEM fuel cell catalyst layer properties from high resolution three-dimensional digital images, Part I: A numerical procedure for the ionomer distribution reconstruction.

27. Yeong CLY, Torquato S. Reconstructing random media. Physical review E 1998; 57 (1): 495-506.

28. Inoue G, Park K, So M, Kimura N, Tsuge Y. Microscale simulations of reaction and mass transport in cathode catalyst layer of polymer electrolyte fuel cell. International J. of Hydrogen Energy 2022; 47: 12665-12683.

29. Cetinbas FC, Ahluwalia RK, Kariuki NN, Myers DJ. Agglomerates in Polymer Electrolyte Fuel Cell Electrodes: Part I. Structural Characterization. J. Electrochem. Soc. 2018; 165: F1051.

30. Zhang X, Ostadi H, Jiang K, Chen R. Reliability of the spherical agglomerate models for catalyst layer in polymer electrolyte membrane fuel cells. Electrochimica Acta 2014; 133: 475-483.

31. Zhang X, Gao Y, Ostadi H, Jiang K, Chen R. Method to improve catalyst layer model for modelling proton exchange membrane fuel cell. Journal of Power Sources 2015; 289: 114-128.

32. Sabharwal M, Pant LM, Putz A, Susac D, Jankovic J, Secanell M. Analysis of Catalyst Layer Microstructures: From Imaging to Performance. Fuel Cells 2016; 16 (6): 734-753.
33. Kaneko T, Yoshimoto Y, Hori T, Takagi S, Ooyama J, Terao T, Kinefuchi I. Relation between oxygen gas diffusivity and porous characteristics under capillary condensation of water in cathode catalyst layers of polymer electrolyte membrane fuel cells. International J. of Heat and Mass Transfer 2020; 150: 119277.

34. Nan N, Wang J. FIB-SEM three-dimensional tomography for characterization of carbon-based materials. Advances in Materials Science and Engineering 2019; 8680715.

35. Cetinbas FC, Ahluwalia RK, Kariuki N, De Andrade V, Fongalland D, Smith L, Sharman J, Ferreira P, Rasoulie S, Myers DJ. Hybrid approach combining multiple characterization techniques and simulations for microstructural analysis of proton exchange membrane fuel cell electrodes. J. Power Sources 2017; 344: 62-73.

36. OpenFOAM, http://www.openfoam.com/

37. Gostick J, Khan ZA, Tranter TG, Kok MDR, Agnaou M, Sadeghi MA, Jervis R. PoreSpy: A python toolkit for quantitative analysis of porous media images. Journal of Open Source Software 2019; 4(37): 1296.

38. Girod R, Lazaridis T, Gasteiger HA, Tileli V. Three-dimensional nanoimaging of fuel cell catalyst layers. Nat. Catal. 2023; 6: 383-391.

39. Lopez-Haro M, Guétaz L, Printemps T. et al. Three-dimensional analysis of Nafion layers in fuel cell electrodes. Nat. Comm. 2014; 5: 5229.

40. Soleymani AP, Reid M, Jankovic J. An Epoxy-Free Sample Preparation Approach to Enable Imaging of Ionomer and Carbon in Polymer Electrolyte Membrane Fuel Cells. Adv. Funct. Mater. 2023; 33: 2209733.

41. Chabot F, Lee J, Vandenberghe F, Guétaz L, Gebel G, Lyonnard S, Porcar L, Rosini S, Morin A. Detailed Catalyst Layer Structure of Proton Exchange Membrane Fuel Cells from Contrast Variation Small-Angle Neutron Scattering. ACS Appl. Energy Mater. 2023; 6(3): 1185-1196.

42. Liu H, L. Ney L, Zamel N, Li X. Effect of catalyst ink and formation process on the multiscale structure of catalyst layers in PEM Fuel Cells. Appl. Sci. 2022; 12(8): 3776.

43. Soboleva T, Zhao X, Malek K, Xie Z, Navessin T, Holdcroft S. On the Micro-, Meso-, and Macroporous Structures of Polymer Electrolyte Membrane Fuel Cell Catalyst Layers, ACS Appl. Mater. Interfaces 2010; 2(2): 375-384.

44. Olbrich W, Kadyk T, Sauter U, Eikerling M. Modeling of wetting phenomena in cathode catalyst layers for PEM fuel cells. Electrochimica Acta 2022; 431: 140850.

45. Whitaker S. The Method of Volume Averaging. Springer Science & Business Media; 2013.

46. Chen L, Zhang R, Kang Q, Tao WQ. Pore-scale study of pore-ionomer interfacial reactive transport processes in proton exchange membrane fuel cell catalyst layer. Chemical Engineering Journal 2020, 391: 123590.

47. Soulaine C, Debenest G, Quintard M. Upscaling multi-component two-phase flow in porous media with partitioning coefficient. Chem Eng Sci. 2011; 66: 6180-6192.

48. Knudsen M. Die Gesetze der Molekularströmung und der inneren Reibungsströmung der Gase durch Röhren. Ann. Phys. (Leipzig) 1909; 28: 75-130.

49. Becker J, Wieser C, Fell S, Steiner K. A multi-scale approach to material modeling of fuel cell diffusion media. International Journal of Heat and Mass Transfer 2011; 54 (7): 1360-1368.

50. Pollard WG, Present RD. On gaseous self-diffusion in long capillary tubes. Phys. Rev. 1948; 73: 762-774.
51. Young JB, Todd B. Modelling of multi-component gas flows in capillaries and porous solids. Int. J. Heat Mass Transfer 2005, 48:5338-53

52. Krishna R, Wesselingh JA. The Maxwell-Stefan Approach to Mass Transfer. Chemical Engineering Science 1997; 199252: 861-911.

53. Quintard M, Bletzacker L, Chenu D, Whitaker S. Nonlinear, multicomponent, mass transport in porous media. Chemical Engineering Science 2006; 61: 2643-2669.

54. Guo J, Quintard M, Laouafa F. Dispersion in Porous Media with Heterogeneous Nonlinear Reactions. Transport in Porous Media 2015; 109: 541-570.

55. Jinnouchi R, Kudo K, Kodama K, Kitano N, Suzuki T, Minami S, Shinozaki K, Hasegawa N, Shinohara A. The role of oxygen-permeable ionomer for polymer electrolyte fuel cells, Nat. Com. 2021; 12, 4956.

56. Paul DK, McCreery R, Karan K. Proton transport property in supported Nafion nanothin films by electrochemical impedance spectroscopy. J Electrochem Soc. 2014; 161 (14): F1395-F1402.

57. Ono Y, Nagao Y. Interfacial structure and proton conductivity of Nafion at the Pt-deposited surface. Langmuir 2016; 32:352-358.

58. Zenyuk IV, Litster S. Modeling ion conduction and electrochemical reactions in water films on thin-film metal electrodes with application to low temperature fuel cells. Electrochim Acta 2014; 146:194-206.

59. Muzaffar T, Kaydyk T, Eikerling M. Physical modeling of the proton density in nanopores of PEM fuel cells catalyst layers. Electrochim Acta 2017; 245 (10): 1048-1058.

60. Guibert R, Horgue P, Debenest G, Quintard M. A Comparison of Various Methods for the Numerical Evaluation of Porous Media Permeability Tensors from Pore-Scale Geometry. Mathematical Geosciences 2016; 48: 329-347.

61. Bamberger A. Approximation des coefficients d'opérateurs elliptiques stables pour la G-convergence. Technical report 1977.

62. Chen L, Wu G, Holby EF, Zelenay P, Tao WQ, Kang Q. Lattice Boltzmann pore-scale investigation of coupled physical and electrochemical processes in C/Pt and non-precious metal cathode catalyst layers in proton exchange membrane fuel cells. Electrochim. Acta 2015; 158: 175-186.

63. Tjaden B, Cooper SJ, Brett DJL, Kramer D, Shearing R. On the origin and application of the Bruggeman correlation for analysing transport phenomena in electrochemical systems. Current Opinion in Chemical Engineering 2016: 12, 44-51.

64. Yu Z, Carter RN. Measurement of effective oxygen diffusivity in electrodes for proton exchange membrane fuel cells. J. Power Sources 2010: 195: 1079.

65. Yu Z, Carter RN, J. Zhang. Measurements of pore size distribution, porosity, effective oxygen diffusivity, and tortuosity of PEM fuel cell electrodes. Fuel cells 2012; 12: 557.

66. Chen L, Wu G, Holby EF, Zelenay P, Tao WQ, Kang Q. Lattice Boltzmann pore-scale investigation of coupled physicalectrochemical processes in C/Pt and non-precious metal cathode catalyst layers in proton exchange membrane fuel cells. Electrochim. Acta 2015; 158: 175.

67. Boyer C, Gamburzev S, Velev O, Srinivasan S, Appleby A. Measurements of proton conductivity in the active layer of PEM fuel cell gas diffusion electrodes. Electrochim. Acta 1998; 43: 3703-3709.

68. Yakovlev YV, Lobko YV, Vorokhta M, Nováková J, Mazur M, Matolínová I, Matolín V. Ionomer content effect on charge and gas transport in the cathode catalyst layer of proton-exchange membrane fuel cells. J. Power Sources 2021; 490: 229531.

69. Havránek, A, Wippermann K. Determination of proton conductivity in anode catalyst layers of the direct methanol fuel cell (DMFC). J. Electroanal. Chem. 2004; 567: 305-315.
70. Suzuki T, Murata H, Hatanaka T, Morimoto Y. Analysis of the catalyst layer of polymer electrolyte fuel cells. R & D Rev. Toyota CRDL 2003; 39(3): 33-38.

71. Du C, Shi P, Cheng X, Yin G. Effective protonic and electronic conductivity of the catalyst layers in proton exchange membrane fuel cells. Electrochem. Commun. 2004; 6: 435-440.

72. Modestov AD, Kapustin AV, Avakov VB, Landgraf IK, Tarasevich MR. Cathode catalyst layers with ionomer to carbon mass ratios in the range 0-2 studied by electrochemical impedance spectroscopy, cyclic voltammetry, and performance measurements. J. Power Sources 2014; 272: 735-742.

73. Liu Y, Jorne J, Gu W. Thin shell model for proton conduction in PEM fuel cell cathodes. J. Electrochem. Soc. 2010; 157(7): B1068.

74. Ren H, Teng Y, Meng X, Fang D, Huang H, Geng J, Shao Z. Ionomer network of catalyst layers for proton exchange membrane fuel cell. J. of Power Sources 2021; 506: 230186.

75. Choi P, Jalani NH, Datta R. Thermodynamics and proton transport in nafion: II. Proton diffusion mechanisms and conductivity. J. Electrochem. Soc. 2005; 152: E123.

76. Artemov VG, Uykur E, Kapralov PO, Kiselev A, Stevenson K, Ouerdane H, Dressel M. Anomalously high proton conduction of interfacial water. J. Phys. Chem. Lett. 2020; 11: 3623.

77. Randrianarizafy B, Schott P, Chandesris M, Gerard M, Bultel Y. Design optimization of rib/channel patterns in a pemfc through performance heterogeneities modelling. Int. J. Hydrogen Energy 2018; 43(18):8907-8926.

78. Shen S, Cheng X, Wang C, Yan X, Ke C, Yin J, Zhang J. Exploration of significant influences of the conditions on the local O2 transport in proton exchange membrane fuel cells (PEMFCs). Phys. Chem. Chem. Phys. 2017; 19: 26221.

79. Xu H, Zheng Z, Hu J, Li H, Cheng X, Guo Y, Luo L, Yin J, Yan X, Wu R, Shen S., Zhang J. An improved agglomerate model with adsorption isotherm for oxygen transport in CCLs with low Pt loading. Int. J. of Heat and Mass Transfer 2024; 222: 125173.