Joel Pauchet

Arno Morin

Thomas David

Pascal Schott

Keywords: PEMFC, MPL, GDL, GDM, X-ray tomography, FIB-SEM, diffusion tensor, numerical simulations

To my parents, to my brothers and sisters, to my cousins and friends » Je tiens tout d'abord à exprimer ma sincère gratitude envers mes directeurs de thèse, Marc Prat et Michel Quintard, pour leur précieuse contribution humaine et scientifique tout au long de ces 3 années. Leur vaste expérience en direction de thèses et en recherche a été d'une valeur inestimable pour l'achèvement de ce projet. Je tiens à remercier Marc pour sa disponibilité et son aide cruciale tout au long de cette période, et particulièrement dans les phases finales de ma thèse, ainsi que Michel pour son soutien continu, ses conseils avisés et ses retours enrichissants. Leur encadrement a grandement contribué à mon développement en tant que chercheur. Je souhaite également remercier chaleureusement les collaborateurs du CEA, notamment Laure Guetaz,

« The scientific discovery lies in the creation of a fortunate and fruitful hypothesis; it is given by the genius of the scientist who created it. »

Claude Bernard (1813 -1878)

Résumé

Les performances des piles à combustible à membrane d'échange de protons (PEMFC) sont étroitement liées à la microstructure de la couche de diffusion de gaz (GDL) comme à celle de la couche catalytique à la cathode (CCL) et à leurs propriétés de transport. Dans ce contexte, l'objectif principal de cette thèse est de caractériser la microstructure et les propriétés des

GDLs et CCLs grâce à l'utilisation combinée de techniques d'imagerie 3D et de simulations numériques.

Dans un premier temps, nous étudions la couche de diffusion (GDL). Celle-ci est une structure composite composée d'un milieu de diffusion fibreux (GDM) couplé à une couche microporeuse (MPL). Nous évaluons les propriétés effectives de diffusion de cet assemblage composite à l'aide de simulations numériques réalisées sur des images tridimensionnelles des microstructures et de l'assemblage GDM/MPL obtenues en combinant tomographie par rayons X (X-Ray) et microscopie électronique à balayage couplée à une technique d'ablation par faisceau d'ions focalisés (FIB-SEM). En utilisant une méthodologie multi-échelle, nous caractérisons le tenseur de diffusion effective de la MPL. Les mécanismes de diffusion de Knudsen et de Fick sont tous deux pris en compte. Pour le calcul du coefficient de diffusion à travers le plan de l'assemblage GDM-MPL une approche hybride est utilisée. Elle combine une description continue pour la MPL avec une prise en compte explicite des fissures de la MPL et des fibres du GDM. Nous examinons également l'impact des fissures de la MPL sur les propriétés de diffusion de l'assemblage à l'aide d'un modèle de résistance à la diffusion considérant l'interpénétration de la MPL dans le GDM, ainsi que l'influence de la compressibilité du GDM.

Dans un second temps, nous nous concentrons sur la couche catalytique cathodique (CCL) en nous appuyant sur des images de la microstructure obtenues par FIB-SEM et par microscopie électronique en transmission. L'imagerie FIB-SEM permet de différencier les pores de la phase solide dans la limite de la résolution considérée (5 nm) mais ne permet pas de différencier le ionomère, le platine et le carbone au sein de la phase solide segmentée.

L'étude commence donc par une étape de reconstruction dans laquelle la distribution du ionomère est reconstruire dans l'image FIB-SEM segmentée en distinguant le ionomère présent sous forme de filaments et le ionomère formant un revêtement sur les particules de carbone. L'identification du ionomer de revêtement repose sur l'idée que les régions concaves sont des sites préférentiels d'adhésion du ionomère sur les particules de carbone lors de l'étape de séchage de la fabrication de la CCL. La procédure développée permet de reconstruire le revêtement du carbone par le ionomère avec des épaisseurs de revêtement variables correspondant aux données expérimentales disponibles. De plus, les particules de platine, qui ne sont identifiables dans les images FIB-SEM, sont réparties à l'aide d'une procédure numérique sur la surface de la phase carbone.

Dans une dernière étape, la microstructure complète de la CCL reconstruite est utilisée pour caractériser par simulations numériques les propriétés effectives associées à cette couche, telles que le tenseur de diffusion de l'oxygène ou de conductivité protonique. Enfin, la reconstruction du ionomère dans l'image effectuée permet de déterminer la distribution spatiale de l'angle de contact à différentes échelles, de l'échelle du grain de platine à l'échelle des pores en mettant en oeuvre une approche liant l'angle de contact à l'épaisseur locale du revêtement en ionomère. La connaissance pour chaque pore de sa taille et de son angle de contact effectif permet aussi l'évaluation de la courbe de pression capillaire de la CCL.

Abstract

The performance of proton exchange membrane fuel cells (PEMFC) is closely linked to the microstructure of the gas diffusion layer (GDL) as well as that of the cathode catalytic layer (CCL) and to their transport properties. In this context, the main objective of this thesis is to characterize the microstructure and properties of GDLs and CCLs through the combined use of 3D imaging techniques and numerical simulations.

In a first step, the diffusion layer (GDL) is studied. This is a composite medium composed of a fibrous diffusion medium (GDM) coupled to a microporous layer (MPL). We evaluate the effective diffusion properties of this composite assembly using numerical simulations carried out on three-dimensional images of the microstructures and the GDM/MPL assembly obtained by combining X-ray tomography (X-Ray) and scanning electron microscopy coupled with a technique of focused ion beam ablation (FIB-SEM). Using a multi-scale methodology, we characterize the effective diffusion tensor of the MPL. Both Knudsen and Fick diffusion mechanisms are considered. A hybrid approach, combining a continuous description for the MPL with an explicit consideration of the cracks of the MPL and the fibers of the GDM, is used for the calculation of the diffusion coefficient of the GDM-MPL assembly in the through plane direction. We also examine the impact of MPL cracks on the diffusion properties of the assembly using a diffusion resistance model considering the interpenetration of MPL in the GDM, as well as the influence of the compressibility of the GDM.

In a second step, the cathode catalyst layer (CCL) studied from images of the microstructure obtained by FIB-SEM and transmission electron microscopy (TEM). FIB-SEM imaging makes it possible to differentiate the pores from the solid phase within the limit of the resolution considered (5 nm) but does not allow to differentiate the ionomer, platinum and carbon within the segmented solid phase. The study therefore begins with a reconstruction step in which the distribution of the ionomer is reconstructed in the segmented FIB-SEM image by distinguishing the ionomer present in the form of filaments and the ionomer coating the carbon particles. The ionomer coating is reconstructed considering that the concave regions are preferential sites of adhesion of the ionomer on the carbon particles during the drying step of the fabrication of the CCL. The developed procedure makes it possible to reconstruct the ionomer coating with variable ionomer coating thicknesses corresponding to the available experimental data. In addition, the platinum particles, which are not identifiable in the FIB-SEM images, are distributed on the surface of the carbon phase using a numerical procedure.

In a last step, the full reconstructed microstructure of the CCL is used to compute effective properties of this layer, such as the oxygen diffusion tensor or the proton conductivity tensor using the same numerical simulation tools developed for the GDL. The ionomer and platinum distribution are also used to determine the local contact angle in the images at various scales, from the platinum particle scale to the pore scale with the help of a theoretical approach relating the contact angle, the local ionomer coating thickness and the interaction between the ionomer and the platinum or carbon substrates. Knowing the size and effective contact angle of each pore also enables us to evaluate the capillary pressure curve of the CCL.

General introduction

In the urgent context of climate change, the imperative of decarbonizing various industries looms large as a crucial step towards mitigating global warming. Recent developments have thrust hydrogen into the spotlight as a promising agent for decarbonization, particularly when sourced from clean energy production methods. This versatile hydrogen resource can be employed in multiple ways, from combustion in thermal engines or furnaces to a more environmentally sustainable approach such as the chemical conversion of hydrogen into electricity through the utilization of fuel cell technology.

Among the diverse range of fuel cell technologies, Proton Exchange Membrane (PEM) fuel cells have emerged as a specific and influential category. These fuel cells harness the principle of proton transport across a specialized membrane that separates two electrodes, generating a difference in electric potential and, consequently, an electric current. Then, a stack of these cells is employed to amplify the electrical voltage output.

However, despite its immense potential, commercial deployment of PEM fuel cell technology encounters significant challenges. Foremost among them is the barrier of cost, a hurdle that has thus far impeded widespread adoption. A significant portion of this cost can be attributed to the catalyst materials used, most notably platinum. Thus, optimizing the loading of catalyst materials becomes imperative, alongside the need to comprehensively understand and address other performance constraints. These constraints typically manifest in the form of limitations related to gas and proton transport to these catalyst surfaces, as well as chemical kinetics.

These transport limitations can be linked to various layers, such as gas diffusion layers (GDLs) and catalyst layers (CLs). Therefore, it is essential to comprehend the microstructure of these distinct layers to optimize their development.

The primary objective of this research endeavor is to quantify these various transport limitations through numerical simulations, with a particular emphasis on processes involving oxygen diffusion, proton transport, and the management of water.

Within the scope of this thesis, our inquiry commences with the characterization of the structural and transport properties of the different layers. Firstly, we investigate the fibrous gas diffusion mediums (GDMs) used in the electrodes to facilitate the ingress of gases such as oxygen and hydrogen while concurrently facilitating the removal of generated water. To enhance the efficacy of these layers, they are often coupled with microporous layers (MPLs) characterized by a finer pore size distribution. Our investigation is specifically directed towards understanding the impact of MPLs on the overall gas diffusion resistance within the fuel cell system.

Subsequently, our research turns its focus towards the active layer, commonly referred to as the catalyst layer (CL), where a complex amalgamation of materials, including carbon, platinum catalyst, and ionomer, plays a pivotal role. Characterizing the structure and properties of this layer becomes indispensable in comprehending the performance limitations associated with proton transport, water management, and gas diffusion. Employing advanced 3D imaging techniques, we scrutinize the microstructure of the catalyst layer. However, owing to inherent limitations in resolution and contrast of 3D images, we have developed a novel numerical approach to identify and delineate the various phases that may elude detection through conventional imaging methods. Then, the newly reconstructed structure serves as a foundation for computing and comprehending the diverse transport properties governing fuel cell operation. The oxygen and proton transport in the CL is evaluated numerically. In addition, the local wetting property is investigated by computing contact angle local distribution.

Before presenting the developed works, the following chapter (chapter 1) is dedicated to providing more general information on the context and PEMFCs. As for the current introduction, the English language in Chapter 1 has been corrected and improved with the assistance of ChatGPT. This assistance has not been employed for the other chapters. It's worth noting that Chapters 2 to 5 are presented in a scientific publication format.

Chapter 1

General context

Since the advent of the industrial revolution, there has been a significant surge in global production, resulting in a dramatic exponential rise in carbon dioxide (CO2) emissions. For a long time, the world had not fully realized the far-reaching consequences of these emissions.

It was only in the 1970s that the first analysis of the global warming emerged, thanks to the pioneering work of the scientists John Tyndall and Svante Arrhenius. Despite these early warnings, the profound impact of global warming largely went unrecognized.

In 1988, a pivotal moment in international efforts to address climate change occurred with the establishment of the Intergovernmental Panel on Climate Change (IPCC). This marked a significant milestone in global cooperation and climate research. In 1992, the United Nations Framework Convention on Climate Change (UNFCCC) was adopted, laying the groundwork for international collaboration in mitigating and adapting to climate change. Subsequently, the Kyoto Protocol in 1997 and the Paris Agreement in 2015 set forth ambitious goals for limiting global temperature increases, with the latter targeting a maximum of 1.5 degrees Celsius.

Achieving these ambitious climate targets necessitates substantial transformations across various industries and lifestyles. Among the most carbon-intensive sectors is transportation, where a comprehensive shift away from hydrocarbon-based combustion engines is imperative. Recent years have witnessed the emergence of diverse solutions to address this challenge. Electrification through batteries has made significant strides in the automotive sector, though it grapples with issues such as charging infrastructure and limited autonomy, particularly for heavy-duty vehicles contributing significantly to transport emissions.

Another innovative solution has arisen, combining the environmental benefits of batteries and electrification with the extended autonomy characteristic of thermal engines. This solution, central to the subject of this thesis, is the fuel cell. Fuel cells are electrochemical systems, akin to batteries, but they differ in being open systems that draw energy from an external storage source, enhancing their autonomy and lifespan.

The roots of fuel cell technology trace back to the 19th century, when Sir William Grove pioneered the concept of water electrolysis [1][2]. In 1839, he developed a functional fuel cell utilizing zinc and platinum electrodes submerged in an electrolyte solution, generating electricity through the reaction of hydrogen and oxygen. This foundational work laid the groundwork for the proton-exchange membrane (PEM) fuel cell, which stands as one of the most extensively studied and utilized variants of fuel cell technology.

Figure 1 illustrates a comparison between various battery technologies and fuel cells. It is evident that fuel cells offer higher energy density but lower output power density compared to batteries. This underscores the advantage of fuel cells in applications such as aviation or aerospace, where energy density is a critical factor. Additionally, as previously mentioned, factors like autonomy and durability are important considerations.

Looking ahead, it's clear that fuel cells and batteries will likely complement each other in hybrid systems. Batteries can serve for energy storage and boosting power density, while fuel cells can handle power generation.

Applications:

PEM fuel cells are widely recognized as the most efficient among various types of fuel cells.

Despite their long history, there has been a remarkable surge in the adoption of PEM fuel cell technology in recent years, particularly in the realm of transportation. This surge is prominently observed in applications such as automotive and aviation.

In the automotive sector, PEM fuel cells play a pivotal role in hydrogen fuel cell vehicles (FCVs) [4], offering a clean and highly efficient means of propulsion. FCVs excel in terms of extended driving ranges and significantly shorter refueling times when compared to electric vehicles relying on battery technology. As depicted in Figure 1, several recent car and bus models incorporate this technology, with companies like Toyota and Hyundai at the forefront of advancements in this market.

Fig. 2. Recent automotive project using PEM fuel cell propulsion [5] Furthermore, fuel cells are currently under consideration as a means to decarbonize aviation transport. Two possibilities are being explored: one involves integrating fuel cells into a hybrid system that combines a hydrogen combustion engine and a hydrogen fuel cell, primarily targeted for long-distance aircraft. Airbus is actively pursuing this approach within its ZEROe program. For shorter-distance aircraft, electric engines powered by hydrogen fuel cells present a promising opportunity. ZeroAvia has already developed an initial model, with successful test flights conducted.

In addition to aviation, PEM fuel cells have recently garnered attention as a power source for boats and maritime cargo vessels. The Energy Observer serves as a test boat that integrates an autonomous energy system relying on renewable energy sources in conjunction with electrolyzers and fuel cells. PEM fuel cells find utility beyond mobile applications; notably, they are valuable in stationary power applications where they can effectively address the intermittency challenge associated with renewable energies. Excess energy production can be harnessed by utilizing water electrolyzers to generate hydrogen, which can then be stored. During periods of lower energy production, PEM fuel cells can seamlessly take the reins, converting this stored hydrogen into electricity. This approach contributes to a more stable and reliable energy supply. Furthermore, fuel cell systems have the potential to supply both electricity and heat for residential use, exhibiting remarkable efficiency and resulting in reduced energy expenses.

Additionally, fuel cells have diverse applications, particularly in the aerospace and defense sectors, where their extended autonomy proves invaluable for critical systems. They serve as dependable backup power sources for essential infrastructures such as data centers and telecommunications equipment, ensuring uninterrupted operation even during power disruptions.

PEM fuel cell structure:

The proton exchange membrane (PEM) fuel cell operates as a system comprising a stack of individual cells. Each cell is constructed with two electrodes, as illustrated in Figure 1, and is sandwiched between two layers known as the Gas Diffusion Layer (GDL). The primary function of the GDL is to facilitate the flow of reactant gases to the Membrane Electrode Assembly (MEA). The MEA consists of two critical catalyst layers: the cathode catalyst layer (CCL) and the anode catalyst layer (ACL).

Fig. 4. Schematic of a PEM fuel cell stack [7].

In the fuel cell's operation, hydrogen gas undergoes oxidation within the anode catalyst layer (oxidation reaction), and then protons are transported across the membrane to react with oxygen within the cathode catalyst layer (reduction reaction), as visually depicted in the figure. The reduction electro-chemical reaction produces water vapor as a byproduct, and it is essential that this condensed water is efficiently removed by the GDL to prevent flooding and maintain optimal performance. The circulation of electrons from the anode to the cathode generates the electric current.

Fig. 5. Scheme of conventional cell [8] Gas Diffusion Layer (GDL):

In certain convention as in Fig. 5, the gas diffusion layer (GDL) is attributed to the fibrous medium. Other conventions use this term for the assembly of fibrous medium with the microporous layer (MPL). Two avoid the confusion, the following terminology is adopted: the fibrous medium is called the gas diffusion medium (GDM) and the GDL is the MPL-GDM assembly.

GDM: This fibrous medium consists of a stack of fibers, each with a diameter typically in the range of 10 µm. The GDM generally exhibits high porosity and anisotropy, meaning that it allows gases to flow through it more readily in certain directions than in others. The fibers are often aligned primarily in the in-plane direction, which is crucial for optimizing gas diffusion and water management within the cell.

Manufacturers such as SGL and Freudenberg are prominent in producing Gas Diffusion Medium (GDM) tailored to diverse fuel cell applications. These materials undergo fabrication

processes that include roll-to-roll techniques and other methods to control properties like porosity, fiber alignment, and structural integrity, ensuring they meet the stringent requirements of fuel cell designs. This layer has a typical thickness on order of hundreds of micrometers.

MPL: The microporous layer (MPL), which is often associated with GDMs in some configurations, is a thin layer applied onto the fibrous medium. The MPL serves to improve catalyst utilization, enhance water management, and further optimize reactant gas distribution within the fuel cell. Hence, a GDL with an MPL is effectively constituted by a combination of the Gas Diffusion Medium (GDM) and the microporous layer (MPL). The overlapping between the GDM and the MPL plays a key role in the diffusion process. In general an important fraction of MPL thickness overlap with the fibrous medium, as visible in Fig. 6. The typical thickness of the MPL is in order of tens micrometers.

Catalyst Layer (CL):

The catalyst layer is the most complex and expensive layer in the cell. Therefore optimizing his performance is crucial. As represented in Fig. 5 and Fig. 6, it is composed by carbon particle agglomerates as support, platinum particle (Pt) as catalyst, and an ionomer as a binder. Protons are conducted by the ionomer to the catalyst surface where the reaction takes place.

Cell performance limitations

The performance of the PEM fuel cell depends on the cell local performance. In general this is measured via the polarization curve as shown in Fig. 7, which represents the cell voltage as a function of current density. Decreasing of voltage is observed as the current intensity increases. This is due to different losses occurring during the cell operation. First losses are due to the activation energy, then conduction losses are observed, and at last diffusion losses occur at higher current density. It can be observed that diffusion losses are the most critical at high current density. Fig. 7. PEM fuel cell polarization curve [10] Considering this decomposition, we can identify the multiphysics processes taking place within the cell. When coupling this with the cell decomposition into different layers (GDM, MPL, CCL), three main processes can be identified: ORR kinetics in the CCL, proton conduction in the CCL, and oxygen diffusion in different layers.

Recent studies in the literature aim to understand the limitations of these different processes and quantify the contribution of each layer. In [12], simulations were conducted to assess the contributions with varying Pt loading. It was found that more than 50% of the performance is attributed to the CCL (Fig. 8a). The distribution of this percentage depends on the Pt loading.

It was observed that ORR kinetic losses increase with higher Pt loading, indicating that higher number of platinum particles become inactive at high Pt loading. Proton transport resistance fraction decreases in this case to 4%. However, oxygen transport limitation fraction remains relatively constant, at about 28%, distributed as 8% due to the CCL and 20% due to the GDL.

Furthermore, it was noted that approximately 23% of the limitation remains unattributed. This a) data from [12]; b) data from [13] If we consider only the transport limitation processes, [13] has simulated their contributions.

Figure 8.b shows that a significant fraction is attributed to oxygen diffusion in the GDL and oxygen transport to the catalyst surface (near Pt). The objective of this work is to characterize by numerical simulations the corresponding effective transport properties using 3D imaging of the different layer.

PEM fuel cell modeling:

The cost and challenges associated with experimental setup and measurements have led many researchers to turn to modeling tools. These numerical models can also serve as valuable complements to experiments for understanding and quantifying the most important parameters in fuel cells.

Starting from 0D and 1D models [14][15] to 2D and 3D computational models [16][17], the progression in the power of numerical simulations has been impressive. In 3D simulations, a significant challenge is modeling the three-dimensional structure of the cell. Some effective 3D models consider equivalent continua for the different layers in the cell. However, in this case, effective properties are required as inputs. To obtain these properties, empirical analytical formulas are sometimes used, or they are derived from experiments in other cases.

However in a multiscale approach, these effective parameters can also be computed directly from simulations, as developed in this work.

Pore-scale models are more accurate than continuum models but are also more computationally costly. Typically, pore-scale models are used to inform the macroscopic continuum models. Some intermediate models, such as pore network models [18][19][20], simplify the complex structures and shapes into a discrete pore network composed of a distribution of pore domains and links between them.

The structure used in pore-scale modeling can be reconstructed numerically [21][22][23].

However, for more accuracy, the real structure distribution can be obtained by 3D imaging technology. This is essential for comprehending the intricate microstructure of PEM fuel cells. Various techniques with varying capacities and resolutions are available, as illustrated in Figure 2. The choice of technique depends on the specific scale of investigation, with each one tailored to its intended purpose.

X-ray computed tomography is generally used to scan the gas diffusion medium, as its resolution matches the scale of this layer. However, for imaging smaller details such as the micro-porosity of the MPL or the catalyst layer, this technology is less suitable. Nano X-ray tomography can provide more details but remains limited. Technologies such as Focused Ion Beam Scanning Electron Microscopy (FIB-SEM) can be highly useful for scanning the details of microstructures, especially as regards the microporous layer and the catalyst layer. One drawback of this method is the destruction of the sample during the image acquisition process.

Additionally, for example it doesn't always allow for the separation of ionomer coating from the carbon phase in the catalyst layer, which is a challenge that will be investigated in this thesis. Other highly resolved methods, such as transmission electron microscopy, can be used to separate close-contrast materials, but only in rather small samples.

Further-FC project

The present work has been developed within the framework of Further-FC (Further Understanding Related to Transport limitations at High current density towards future Electrodes for Fuel Cells). Further FC is a European funded project including several collaborators in Europe and beyond as listed in Fig. 10. The overarching goal of this project is to advance the production and adoption of fuel cell vehicles. This is achieved by enhancing the specifications and design of the Membrane Electrode Assembly (MEA). The aim is to develop MEAs that not only perform at a higher level but also exhibit greater durability, thus contributing to the overall improvement of PEMFC technology for fuel cell vehicles. The project focuses on the catalyst layer but other layers have been also investigated as in the present work.

Within this project, our work has developed in special collaboration with CEA-Liten. CEA-Liten notably obtained the FIB-SEM 3D imaging of the CCL and MPL. Additionally, the other partners were involved in discussions during various project meetings Fig. 10. Further-FC project organization [11].

Outline of PhD thesis manuscript

This PhD thesis manuscript is organized as follows:

This present chapter provides an introduction to the global context, challenges, and objectives of the research. The investigation and results are then presented in chapters 2 through 5.

In Chapter 2, we present our investigation on the gas diffusion layer. The microstructure of the gas diffusion layer is scanned using X-ray tomography, while the microstructure of the MPL is analyzed using FIB-SEM. By coupling these two imaging technologies, we develop a numerical multiscale approach to compute the oxygen diffusion coefficient at different scales and evaluate the various diffusion resistances at the macroscale associated with the MPL and the diffusion medium (GDM).

In Chapter 3, we shift our focus to the catalyst layer, more precisely the cathode catalyst layer (CCL). Due to the pore submicronic size, FIB-SEM imaging is used to obtain 3D digital images of the CCL microstructure, which consists of carbon, platinum, and ionomer.

However, FIB-SEM does not allow for the separation of these different phases. Therefore, in this chapter, we develop a reconstruction model to complement FIB-SEM with missing details such as ionomer and platinum spatial distributions.

In Chapter 4, the 3D images with the reconstructed ionomer distribution in the CCL microstructure from Chapter 3 are used to study the effective proton conductivity tensor and effective oxygen diffusion tensor in the CCL from direct numerical simulations on the reconstructed structure. The results are compared with experimental data from the literature.

Chapter 5 is devoted to a study of the wetting property of the cathode catalyst layer using again the reconstructed images from Chapter 3. The local contact angle spatial distribution in the images at a scale consistent with the size of the platinum particles is obtained. The wettability is also characterized at the pore scale via the concept of effective contact angle.

The obtained results are used to evaluate the capillary pressure curve.

Chapter 6 is a global conclusion and recapitulation of results and discussions.

Introduction

Polymer Electrolyte Membrane Fuel Cells (PEMFCs) are one of the most promising technologies as power source for use in a wide range of applications in the context of the energy transition and the deployment of clean energy systems. Nevertheless, improvements in terms of performance, durability and cost effectiveness must be achieved. This can be expected from improved design of an elementary cell and its components. The properties of the various components are therefore of paramount importance in this context. In relation with the water management, a key issue in the operation of PEMFC [1], it has been shown that inserting a microporous layer (MPL) between the gas diffusion medium (GDM) and the catalyst layer (CL) improves the performance [2][3][4][5][6][7][8][9]. As pointed out in [10], the improvement brought by the MPL is, however, not fully understood. A possible positive impact of the MPL is to prevent the liquid water forming by condensation in the GDM [11,12] to reach the catalyst layer and to block the reaction sites. It has also been argued that the MPL might reduce the saturation in the GDM by reducing the number of liquid pathways through the GDM [13,14,15,16]. The MPL might also act as a capillary barrier for the liquid water forming in the catalyst layer, which would be favorable for the membrane hydration, e.g. [17]. Inserting a MPL might also be beneficial in reducing the thermal and electrical contact resistances between the GDM and the catalyst layer [18]. Also, cracks are generally present in the MPL. However, the role of the cracks in the performance improvement due to the MPL is still unclear. Cracks might provide preferential liquid water pathways limiting the number of liquid water entry points in the GDM, which could be beneficial to the oxygen access toward the cathode catalyst layer. On the other hand, it has been reported that cracks might degrade the performance [19,20]. In this context, characterizing the MPL, the GDM and the GDM-MPL assembly in terms of structure and properties is a natural step in order to better understand the impact of the MPL and the GDM-MPL assembly. As reviewed for instance in [21], numerous studies have been performed to characterize the properties of various GDLs with or without MPL. However, as pointed out in [22], the studies allowing to differentiate the impact of the GDM and the MPL are much scarce. A reason is that the mean pore size in the MPL is in the submicronic range (on the order of a hundred of nanometers)

whereas the pore size in the GDM are much bigger (on the order of a few tens of micronmeters). As exemplified in [22], analyzing the GDM-MPL assembly then imposes to combine various techniques when the characterization is sought from high-resolution 3-D imaging techniques, a powerful and increasingly used approach to characterize porous materials [23], [24]. In this work, the GDM-MPL assembly is characterized by combining focused ion beam -scanning electron microscopy (FIB-SEM) [24] and X-ray computed tomography (XCT) [23]. The MPL is characterized via FIB-SEM thanks to the high resolution achievable with this technique with voxel size of 5 nm whereas the GDM fibrous material and the cracks in the MPL are characterized via XCT with 1 μm voxel size. Using FIB-SEM combined or not with XCT to characterize the MPL or the GDM-MPL assembly is not a novelty. In [25], MPL digital images were obtained via FIB-SEM and MPL permeability and tortuosity were determined from Lattice Boltzmann Model (LBM) simulations. However, the numerical simulations were restricted to the MPL. Furthermore, the cracks were not characterized, nor considered in the simulations. A more comprehensive characterization study of the GDM-MPL assembly combining FIB-SEM and XCT was developed by Wargo et al. [22]. However, Knudsen diffusion was not considered and only the through-plane direction was considered in the characterization of the diffusive transport.

FIB-SEM imaging was also used in [26] to characterize the MPL. Only the MPL was characterized. The MPL effective diffusion coefficient was determined but only in the through-plane direction and the impact of the cracks was not considered. FIB-SEM and XCT were also combined in [27] to characterize the GDM and the MPL microstructures. The tortuosity was characterized using a random walk method. However, the cracks in the MPL were not considered and no direct numerical simulations considering the Knudsen diffusion were performed. In [28], MPL digital 3-D images were obtained using a dedicated nano-scale X-ray computed tomography technique. Only a crack-free MPL subdomain was considered.

By contrast to these studies, the MPL, GDM and the GDM-MPL assembly with full consideration of the MPL cracks are considered in the present work. The effective oxygen diffusion tensor of the MPL is computed considering Knudsen and Fick diffusion. The through-plane diffusive property of the GDM-MPL assembly is characterized with a special attention devoted to the impact of MPL cracks. All the computations are performed considering dry media, i.e. in the absence of liquid water in the pores.

The paper is organized as follows. The imaging of the MPL and the GDM-MPL assembly via FIB-SEM and XCT is presented in Section 2. The numerical methodology used to compute the effective diffusion properties is presented in Section 3. Results of simulations are presented in Section 4. A discussion is proposed in Section 5. Section 6 consists of the main conclusions of the study.

MPL FIB-SEM and GDM-MPL XCT imaging

A Sigracet ® SGL 22BB sample was selected to perform the study and analyze the MPL, the GDM and the GDM-MPL assembly.

MPL microstructure via FIB-SEM

3D FIB-SEM data was acquired on a Zeiss Crossbeam 550. For the MPL, a volume of 10 µm by 10 µm by 7.5 µm was imaged with pixel size of 5 nm and inter-slices distance of 5 nm too, resulting in cubic voxel of 5x5x5 nm 3 . The FIB milling was done at 30kV with a current of 700 pA, and electron column was operated at 1.2 kV and 1 nA. Images from InLens and SE2 detectors were acquired simultaneously. Image processing before segmentation included slice-to-slice alignment, homogenization of contrast/brightness in all three directions, correction of curtaining effect and resampling in z direction to correct for non-homogeneous inter-slice distances during milling and imaging. These steps were done thanks to homemade python scripts based on usual image processing python libraries (e.g. numpy, scikit-image…). Segmentation was done thanks to machine learning based on random forest classifier with ilastik software [29] using the 'pixel classification' workflow. The annotation and training phase was done on a sub-volume of 192 x 196 x 400 pixels.

The 3D segmented image of the MPL microstructure obtained via FIB-SEM is shown in Figure 1a. In order to characterize the heterogeneity of the sample, the image is partitioned in adjacent cubes of size L. The MPL porosity is characterized by computing the void volume fraction in each cube. This operation is performed varying the cube side length between 250 nm and the image full size (7.5 µm x 10 µm x 10 µm). The obtained results are shown in Fig. 1b and illustrate the heterogeneous nature of the MPL. The porosity varies from one cube to the other for all cubic block sizes considered suggesting that a representative averaging volume (REV) can be hardly identified in the considered range of sizes. However, the variance of the porosity results for each cube size decreases monotonously with increasing cube size. This is in favour of the emergence of some REV, which should not be far from the full sample size, given the observed variance decreasing rate. The MPL porosity for the entire image is 0.71, which is significantly larger than the MPL porosity of 0.42 reported in [22] for the SGL 10BC and the value of 0.5 reported in [28] for the SGL 24BC. The value of 0.71 is in very good agreement with the value of 0.69±0.03 reported in [10] for the microporosity of the same GDL, i.e. a SGL 22 BB. In accordance with [10], it is concluded here that the MPL microporosity considered in some previous studies [30][31][32][33] for Sigracet ® materials was too low. In this respect, it is interesting to notice that the characterization technique used in [10] is X-ray tomographic microscopy, and not FIB-SEM. The fact that both techniques leads to comparable results is interpreted as a good signal of result consistency.

The MPL pore size distribution is characterized with the image analysis toolbox Porespy [34] using the Maximum Sphere Inscription (MSI) method. The mean pore diameter is about 150 nm, which is about half the mean pore size reported in [10] for the SGL 10BC from the chord length distribution and about twice the average pore size (90 nn) reported in [28] for the SGL 24BC but similar to the mean pore size reported in [26] for a Johnson Matthey MPL [25].

GDM-MPL assembly via XCT

The body of the MPL matrix being characterized through the FIB-SEM techniques, µm-size structures (cracks, GDM) have been characterized using XCT. To this end, samples were 3Dscanned using a Phoenix Nanotom tomograph, at high spatial resolution, with an isotropic voxel size of 1 µm. Two triangular pieces of the GDM-MPL assembly were cut with a chisel and inserted in a glass capillary tube with a 2 mm inner diameter. Two pieces of GDM-MPL assembly were used in order to have a high-enough contrast during the alignment procedure of the capillary tube, source and flat panel detector, which relies on live visualizations of Xray radiographs. Only the top part of the samples, which stood slightly above the capillary tube top end, was exposed to the X-ray beam and thus 3D scanned. Image acquisition was performed at a 50 kV voltage and 275 µA amperage. 1440 radiographies were taken with an exposure time of 1000 ms and 4 frames averaging for a total scan time of around 2 hours.

Correction of ring artefacts and reconstruction was performed using Phoenix Datos|X software and tomographic slices were exported for further image processing. Avizo software was then use to denoise the tomographic slides (using a 3D non local mean filtering algorithm) and to segment them, i.e. to separate the different constitutive materials.

The resolution of the XCT image (~1 μm) does not allow to capture the sub-micronic pores in the MPL but allows to capture the GDM-MPL microstructure with an overall size adequate to characterize the GDM fibrous material, the "diffuse" MPL-GDM interface, and the crack network in the MPL. The 3D image of the GDM-MPL assembly obtained via XCT is illustrated in Fig. 2. The MPL is shown in red and the GDM in green in Fig. 2.

Cracks in the MPL are clearly visible in Fig. 2. The crack pattern is similar to those reported in previous works, e.g. [10,[34][35][36][37]. In agreement with previous works, i.e. [35], the majority of these cracks are traversing the MPL thickness as illustrated in Figure 2e. The mean crack width in the image shown in Fig. 2f is 14 µm whereas the mean in-plane distance between cracks is 196 µm. These lengths, illustrated in Fig. 2f, correspond to maximum inscribed circles in the cracks or in the regions between cracks determined using Porespy [34]. The corresponding distributions are shown in Fig. 2f.

The porosity of the fibrous material, i.e. the GDM, is estimated as 0.88, which is consistent with previous works, e.g. [22], for uncompressed GDL whereas the porosity of the GDM-MPL assembly is 0.81. This value is obtained from the XCT image assigning the value 0.71 (obtained from the FIB-SEM image) to the MPL voxels in the XCT image. As shown in Fig. 2f, the pore sizes in the GDM determined using Porespy [34] are of a few tens of micrometers, which is also consistent with previous works, e.g. [22], and about 2-3 orders of magnitude larger than in the MPL (Fig. 1c). Fig. 2f also shows the variation of the macro-porosity in the through-plane direction. The macro-porosity is the porosity determined in successive slices of thickness 1 µm (the XCT image voxel size) from the XCT images. The MPL matrix is considered as a solid phase in these images. Hence, the porosity of the MPL matrix is not included in the macro-porosity.

The latter corresponds to the cracks in the MPL and the void space between the fibers in the GDM. Note that the XCT image includes a zone free of MPL (zone corresponding to z approximately smaller than 30 µm in Fig. 2f.). The plateau for 50𝜇𝑚 ≤ 𝑧 ≤ 90𝜇𝑚 in Fig. 2f corresponds to the cracks in the MPL. This plateau is consistent with the fact that the cracks are weakly tortuous crossing the MPL. The zone further on the right in Fig. 2f corresponds to the GDM.

Oxygen diffusion numerical computation

The objective is to characterize the diffusion property of the GDM-MPL assembly from numerical simulations on the 3D-images. Due to the quite large difference in the pore sizes between the GDM and the MPL a fully direct approach over the computational domain illustrated in Fig. 2c is not doable (using the FIB-SEM voxel (~of size 5 nm) as elementary computational control volume, the direct computation over the domain shown in Fig. 2c would correspond to a much too large computational domain with more than 10 15 grid points!). For this reason, the solution is sought via a multiscale approach. The effective diffusion tensor is first computed for the MPL matrix, i.e. the MPL crack free region, which corresponds to the 3D image obtained via FIB-SEM (Fig. 1a). The obtained data are then used to perform the computation over the GDM-MPL assembly image obtained via XCT (Fig. 2c), thus considering the MPL on the GDM surface and the fraction of the MPL that has penetrated between the GDM fibers. At the scale of the GDM-MPL assembly, the MPL is therefore considered as a continuum material containing a discrete network of cracks. Since the computations are based on the images and the segmented images are obtained via a series of operations (raw data acquisition, image processing, segmentation), the question arises as to whether the computational results are sensitive to these operations. However, owing to the computational cost of the simulations, a sensitivity study has not been performed. The final images obtained at the end of the segmentation step were judged as sufficiently representative of the real structures. The fact that the obtained porosity for instance is consistent with recent work in the literature, e.g. [10] as regards the MPL, is also an indication of the representativeness of the images. Furthermore, it can be noted that the transport phenomenon considered in the present study, i.e. diffusion, is generally considered as the less sensitive to geometrical details (roughness, corners, etc) among the transport processes of interest in a PEMFC.

Pore scale diffusion transport governing equation

In the absence of chemical reaction and other mechanisms that could enhance the multicomponent character of the diffusion process, it is acceptable to describe the oxygen transport by diffusion at the pore scale through the Fickian diffusion equation

𝜕𝑐 𝜕𝑡 = ∇. (𝐷(𝑥)∇𝑐) ( 1
)
where c is the oxygen concentration and D is the local diffusion coefficient. Given the small pore sizes in the MPL, Knudsen diffusion [38] must be considered in conjunction with Fick's diffusion. As in several previous works, e.g. [39] and references therein, our approach is based on the Bosanquet relation [40] combining the molecule-wall and molecule-molecule interactions according to the formula,

𝐷(𝑥) = 1 𝐷 𝑏𝑢𝑙𝑘 -1 +𝐷 𝐾𝑛 -1 (2)
where 𝐷 𝑏𝑢𝑙𝑘 is the bulk diffusion coefficient and 𝐷 𝐾𝑛 is the Knudsen diffusion coefficient.

The latter is expressed as

𝐷 𝐾𝑛 = 𝑑 𝑝 3 √ 8𝑅𝑇 𝜋𝑀 (3)
where R, M, T, d p are the ideal gas constant, molar mass, temperature and pore diameter, respectively. A first hint on the significance of Knudsen diffusion can be gained from the computation of the Knudsen number,

𝐾𝑛 = ℓ 𝑑 𝑝 (4)
where

ℓ = 𝑅𝑇 √2𝜋𝑑 𝑚𝑜𝑙 2 𝑁 𝐴 𝑃
is the mean free path (P is the pressure, 𝑑 𝑚𝑜𝑙 is the molecular diameter and 𝑁 𝐴 is the Avogadro number). As illustrated in Figure 3a, the Maximum Sphere Inscription (MSI) method is used to determine the local pore diameter 𝑑 𝑝 using Porespy [34].

The method, e.g. [41], consists in determining for each pore voxel the maximum sphere radius that can include this voxel. This leads to segment the pore space into a set of several spheres of different characteristic length (diameter). Fig. 3b shows the variation of the Knudsen number so obtained for standard PEMFC operating conditions (P ~1.5 bar, T ~ 80°C). As can be seen, the Knudsen number is lower than 1 except in a few pores in the MPL, where Kn can be on the order of 1.2. According to [42], the Bosanquet formula Eq. ( 2) can be used with some confidence in this range of Knudsen number, based on Molecular Dynamics simulations for one dimensional channels.

Nevertheless, doubts subsist about the accuracy of the approach for complex pore structures.

Addressing this problem would require heavy MD simulations, which are beyond the scope of the present work. This must be kept in mind when discussing the results. Only significant differences can be considered as non-questionable. It is important to notice that the nonnegligible values of the Knudsen number in the MPL make the coefficient 𝐷(𝑥) in Eq.( 1)

dependent on the considered location in the pore space.

Upscaled diffusion transport governing equation in MPL matrix (outside the cracks)

According to the volume averaging method [43], the upscaled diffusion transport governing equation reads,

𝜀 𝜕𝐶 𝜕𝑡 = ∇. (𝜀𝐃 𝑒𝑓𝑓 . ∇𝐶) ( 5
)
Where 𝜀 is the MPL local porosity and 𝐃 𝑒𝑓𝑓 is the effective diffusion tensor; 𝐶 can be interpreted as the intrinsic average concentration over the averaging volume, see [43] for more details. Within the framework of the volume averaging method, the derivation of Eq.( 5)

is associated with the use of a small parameter, which is the ratio of the averaging volume size to the size of the porous medium domain (the MPL here) or equivalently with the existence of a representative elementary volume (REV) whose size is small compared to the porous medium size. In this respect, based on the results on the porosity shown in Fig. 1b, the validity of Eq.( 5) can be questioned since no clear REV can be inferred from these results. This is a common situation encountered when dealing with heterogeneous materials [44,45].

Following these approaches, resolution of diffusion problems at the MPL scale is carried out by a grid-block technique, which can be described as follows: Eq. ( 5) is resolved over a coarse numerical grid with effective diffusion tensors calculated for every grid-block based on a finer description of the structure. Because of the lack of true REV size, one can play with the coarse grid-block size to control the accuracy of the simulation. There are many advantages to this method since one does not try to homogenize the entire MPL domain.

Sub-grid-block size heterogeneity effects are incorporated in every gridblock effective tensor, while coarse-scale heterogeneity effects are kept, thus, for instance, allowing to reproduce some abnormal (non-Fickian) diffusion effects induced by these heterogeneities.

Incidentally, but this is very important on a practical basis, fine calculations are only carried out on the subgrid-blocks, thus limiting the need for huge computational resources. This gridblock approach can be reproduced sequentially, if needed because of computational requirements. In our case, as for the porosity characterization (Section 2), the grid-block step consists in partitioning the MPL matrix digital image in adjacent cubic blocks of equal size L. In the next section, we explain how the effective diffusion tensor is computed for a given cubic block.

Effective diffusion tensor computation over a block.

A somewhat commonly used approach to determine an effective property in one direction (often the through-plane direction in PEMFC component characterization) is to impose Dirichlet conditions, i.e.uniform concentrations over the inlet and outlet computational domain limiting surfaces perpendicular to the considered direction and spatially periodic boundary conditions in the lateral directions, e.g. [22,26].

However, the MPL 3D image (Fig. 1a) or a sub-block extracted from the image are not spatially periodic and do not necessarily correspond to a spherical effective diffusion tensor.

The image (or the sub-block) can be rendered spatially periodic [46] but this significantly increases the size of the image and also alter the anisotropy properties. For this reason, following [46], spatially periodic boundary conditions are not used in what follows.

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 will be estimated in order 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 (as sketched in Fig. 4 with c = c 1 at the outlet and c = c 2 at the inlet) whereas two types of lateral conditions are considered: i) zero-flux conditions, ii) Bamberger conditions [47]. As explained in [46], the Bamberger condition consists in imposing a constant concentration The steady-state version of Eq.( 1) together with the boundary conditions discussed previously is solved using the LaplacianFoam solver of the open source CFD software OPENFOAM. Note that the fact that oxygen diffusion does not occur in the solid phase is taken into account via a penalization technique consisting in imposing a very low value of the diffusion coefficient in the solid phase (~10 -9 m 2 /s, which is 10 4 smaller that the diffusion coefficient in the pores).

Once the pore scale concentration field c is obtained, the average diffusive flux over the computational domain is computed as

〈𝒋〉 = 1 𝑉 ∫ (-𝐷(𝒙)∇𝑐) 𝑑𝑉 𝑉 𝑓 (6)
where 𝑉 is the volume of the computational domain and 𝑉 𝑓 is the volume occupied by the pore space in V.

Then, the effective diffusion tensor must be determined according to Step 2 described above.

This step requires that the effective diffusion tensor is estimated so as to reproduce the micro-scale average fluxes. In principle, one should solve Eq. ( 5) and use some optimization method to estimate 𝐃 𝑒𝑓𝑓 . However, approximate methods are available to avoid such an optimization method. They are based on approximate solutions of diffusion over the block in the case of a homogenoeus effective diffusion tensor. Let us start with the Bamberger conditions. There is a simple analytical solution that corresponds to homogeneous diffusion in a slab of infinite lateral extent. The macro-scale concentration gradient, ∇C, is uniform, as well as the macro-scale flux, J. The components of the effective diffusion tensor are then determined by expressing that the average flux given by Eq. ( 6) must be equal to the macroscopic diffusion flux, namely 𝑱 = -𝜀𝐃 𝑒𝑓𝑓 . ∇𝐶 = 〈𝒋〉 (7) For the configuration illustrated in Fig. 4, the macroscopic gradient ∇𝐶 corresponds to

∇𝐶 = 𝑐 1 -𝑐 2 𝐻 𝑥
𝒆 𝑥 , where 𝒆 𝑥 is the unit vector in the direction of x and H x is the dimension of the computational domain in the x direction. From Eq. (7), this leads to the following relationships

𝐷 𝑒𝑓𝑓𝑥𝑥 = - 𝐻 𝑥 𝜀(𝑐 1 -𝑐 2 )
〈𝑗 𝑥 〉 (8)

𝐷 𝑒𝑓𝑓𝑦𝑥 = - 𝐻 𝑥 𝜀(𝑐 1 -𝑐 2 )
〈𝑗 𝑦 〉 (9)

𝐷 𝑒𝑓𝑓𝑧𝑥 = - 𝐻 𝑥 𝜀(𝑐 1 -𝑐 2 ) 〈𝑗 𝑧 〉 (10)
The same procedure is applied in the y and z directions, which allows computing all the components of tensor 𝐃 𝑒𝑓𝑓 ,

( 𝐷 𝑒𝑓𝑓𝑥𝑥 𝐷 𝑒𝑓𝑓𝑥𝑦 𝐷 𝑒𝑓𝑓𝑥𝑧 𝐷 𝑒𝑓𝑓𝑦𝑥 𝐷 𝑒𝑓𝑓𝑦𝑦 𝐷 𝑒𝑓𝑓𝑥𝑧 𝐷 𝑒𝑓𝑓𝑧𝑥 𝐷 𝑒𝑓𝑓𝑧𝑦 𝐷 𝑒𝑓𝑓𝑧𝑧 ) 𝐿 1 = - 1 𝜀 𝐿 1 (𝑐 1 -𝑐 2 ) ( 𝐻 𝑥 〈𝑗 𝑥 〉 1 𝐻 𝑦 〈𝑗 𝑥 〉 2 𝐻 𝑧 〈𝑗 𝑥 〉 3 𝐻 𝑥 〈𝑗 𝑦 〉 1 𝐻 𝑦 〈𝑗 𝑦 〉 2 𝐻 𝑧 〈𝑗 𝑦 〉 3 𝐻 𝑥 〈𝑗 𝑧 〉 1 𝐻 𝑦 〈𝑗 𝑧 〉 2 𝐻 𝑧 〈𝑗 𝑧 〉 3 ) (11)
where indices 1, 2 and 3 correspond to the application of the inlet-outlet concentration difference (𝑐 1 -𝑐 2 ) in the x, y and z directions, respectively. It can be noticed that here 𝐻 𝑥

=𝐻 𝑦 = 𝐻 𝑧 = 𝐿 1 .
There is no such an analytical solution for the computations with the zero flux lateral conditions. However, the concentration and flux fields far from the lateral boundaries are close to the Bamberger solution. Therefore, one will use the same estimation procedure as described by Eqs. (11). It must be noticed that these estimations do not impose that the effective tensor is symmetric. This is not an inconvenience, as long as the effective tensor is definite positive (see discussion in [46]).

Sequential methodology

Given the MPL image overall size (5.6 billion voxels), it is not numerically efficient to apply the procedure presented in Section 3.3 on the whole image owing to the long computational times and memory size requirements. As mentioned before, a sequential approach is adopted.

This avoids very long computation times and, also, allows analyzing the distribution of effective coefficients at intermediate Darcy-scales. This is useful to discuss the existence of a REV (Representative Elementary Volume), analyze spatial heterogeneities, etc. In this work a two scale sequential procedure is developed. As displayed in Fig. 4e, the whole image is partitioned into blocks of a chosen size L 1 . The effective tensor 𝐃 𝑒𝑓𝑓-𝐿 1 is computed for each block of size L 1 according to the procedure described in Section 3.3 (Eqs. (6)(7)(8)(9)(10)(11)). The diffusion problem at the scale L 2 of the whole image is then formulated as

∇. (𝜀 𝐿 1 (𝒙)𝐃 𝑒𝑓𝑓-𝐿 1 (𝒙). ∇𝐶) = 0. ( 12
)
with appropriate boundary conditions. In this approach, the effective diffusion properties determined for one block and is 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. ( 12) is solved numerically using the same boundary conditions as presented in section 3.3 considering successively the directions x, y and z to apply the inlet-outlet concentration difference (𝑐 1 -𝑐 2 ) and the effective tensor components for the whole image (scale L 2 ) are determined using relationships similar to Eq.( 11), namely

( 𝐷 𝑒𝑓𝑓𝑥𝑥 𝐷 𝑒𝑓𝑓𝑥𝑦 𝐷 𝑒𝑓𝑓𝑥𝑧 𝐷 𝑒𝑓𝑓𝑦𝑥 𝐷 𝑒𝑓𝑓𝑦𝑦 𝐷 𝑒𝑓𝑓𝑥𝑧 𝐷 𝑒𝑓𝑓𝑧𝑥 𝐷 𝑒𝑓𝑓𝑧𝑦 𝐷 𝑒𝑓𝑓𝑧𝑧 ) 𝐿 2 = - 1 〈𝜀〉 𝑳 𝟐 (𝑐 1 -𝑐 2 ) ( 𝐻 𝑥 〈𝑗 𝑥 〉 1 𝐻 𝑦 〈𝑗 𝑥 〉 2 𝐻 𝑧 〈𝑗 𝑥 〉 3 𝐻 𝑥 〈𝑗 𝑦 〉 1 𝐻 𝑦 〈𝑗 𝑦 〉 2 𝐻 𝑧 〈𝑗 𝑦 〉 3 𝐻 𝑥 〈𝑗 𝑧 〉 1 𝐻 𝑦 〈𝑗 𝑧 〉 2 𝐻 𝑧 〈𝑗 𝑧 〉 3 ) (13)
where 𝐻 𝑥 = 10 µm, 𝐻 𝑦 =10 µm, 𝐻 𝑧 =7.5 µm (Fig. 1a) and

〈𝒋〉 = 1 𝑉 ∫ (-𝜀 𝐿 1 (𝒙)𝐃 𝑒𝑓𝑓-𝐿 1 (𝒙). ∇𝐶 )𝑑𝑉 𝑉 (14)
As for the first step (Section 3.3), Eq. ( 12) is solved using the LaplacienFoam solver in OpenFOAM. A structured hexahedral mesh is used proportional to the voxel matrix. A subdivision factor allows refining or coarsening the mesh, depending of the size of the blocks. Computations are performed using Olympe super-calculator from the regional computational center "Calmip". The typical CPU time is about 400h to obtain the final effective diffusion tensor at scale L 2 .

Results

The multiscale modeling approach imposes to calculate first the MPL properties (outside the cracks) in order to use them in the computation of the GDM-MPL assembly. Therefore, the results are presented in two parts. The MPL diffusion properties are presented first and the impact of anisotropy and Knudsen diffusion are evaluated. Then, the diffusion properties of the GDL with and without MPL are presented. Finally, the results for the GDL-MPL assembly are presented and the impact of MPL penetration into the GDL and MPL cracks are investigated.

MPL diffusion properties

In this section, the methodology described in Sec. 3 is used.

MPL effective diffusion tensor at block scale L1

As illustrated in Fig. 4e, the MPL effective diffusion tensor is calculated for two block sizes, namely 500 nm and 1500 nm.

The distribution of the effective diffusivity components for both block sizes is shown in Fig. 5 with or without consideration of Knudsen diffusion. It can be first noticed that the porosity varies significantly from one block to the other. The diagonal components globally increase with the porosity but the spreading of the results for a given intermediate porosity is significant when the Knudsen effect is not considered. This is another illustration of the lack of REV at this scale. The effective diffusivity tensor non-diagonal components are quite small compared to the diagonal components (less than 10% for the larger block size). The values of the tensor diagonal components are slightly different for the considered block sizes, indicating an anisotropy effect. This anisotropy effect might be due in part to the shadowing effect in the FIB-SEM technique [48]. Since this anisotropy effect is small, this point has not been investigated further. Nevertheless, it can be noted that the anisotropy of the microstructure could be also characterized directly from the images by computing the chord length distributions, e.g. [49]. Interestingly, the relatively small value of the off-diagonal components suggests that the principal axes are aligned with the coordinate axes. The spreading of the results is a bit less for the largest block size but remains noticeable.

As can be seen from Fig. 5, the impact of Knudsen diffusion is significant. One notice that the spreading of the results for a given porosity is much less. Also, the off-diagonal components are quite negligible compared to the diagonal components. Fig. 5e gives a more detailed view of the first component distribution of the effective diffusivity tensor, i.e. 𝐷 𝑒𝑓𝑓𝑥𝑥-𝐿 1 , for both block sizes. This plot better illustrates the impact of Knudsen diffusion with a mean value with Knudsen diffusion that is about half the value obtained when only Fick diffusion is considered. Also, it can be seen that there is almost no overlap between both distributions for the largest block size considered. This figure also illustrates the impact of the block size on the spreading of the distribution. The larger the block size, the smaller is the spreading. In this paragraph, the computed relative diffusion coefficient at this scale, i.e.

𝐷 𝑒𝑓𝑓𝑥𝑥-𝐿 1 𝐷 𝑏𝑢𝑙𝑘 is compared with the experimental results reported in [27] considering the same pressure (P = 201 kPa) and temperature (T = 22°C) as in [27]. As explained in [27], the experimental results were also used to extract the relative diffusion coefficient when the Knudsen diffusion is not considered. As can be seen from Fig. 6, our simulation results are in a quite good agreement with the experimental results of Inoue and al. [27] for both cases, with and without consideration of the Knudsen diffusion. It can be noticed that the MPLs considered in [27] are from Gore TM and not from SGL. Nevertheless, the results are quite comparable in Fig. 6. In the absence of Knudsen diffusion, the relative diffusion coefficient is classically expressed in terms of tortuosity, i.e. . In the fuel cells literature, the tortuosity 𝜏 is often estimated using Bruggeman's relationship, 𝜏 = 𝜀 -0.5 [50], noting that the exponent value 0.5 applies for a medium made of spherical particles [50]. In line with [50], it can be seen that Bruggeman relationship with the exponent -0.5 overestimates the relative effective diffusion coefficient.

Nevertheless, it can be seen from Fig. 6 that the exponent can be adjusted so as to represent reasonably well the data. When the Knudsen diffusion is considered, the exponent actually depends on the considered pressure, temperature and gas.

It can be also noticed that our simulation results with Knudsen diffusion are in quite good agreement with those of Guo et al. [51] obtained from a lattice-Boltzmann model on a numerically reconstructed MPL microstructure. The following relationship is used in [51] to relate the through plane effective diffusion coefficients with or without consideration of Knudsen diffusion,

𝐷 𝑒𝑓𝑓𝑥𝑥-𝐾𝑛 𝐷 𝑒𝑓𝑓𝑥𝑥-𝐹𝑖𝑐𝑘 = 1 1+𝜒〈𝐾𝑛〉 ( 15
)
where 𝜒 is a numerical factor and 〈𝐾𝑛〉 is the average Knudsen number defined as the ratio of the gas mean free path to the average pore diameter in the considered block. Adjusting parameter 𝜒 (𝜒 = 2.8) leads to the data shown with green stars in Fig. 6, in quite good agreement with the pore scale computation data obtained from the direct simulations with consideration of the Knudsen diffusion. In agreement with [51], Eq.( 15) provides an effective analytical model to estimate simply the impact of Knudsen diffusion on the gas effective diffusion coefficient. It can be noted however that the adjusted value of the numerical factor 𝜒 is greater than the value reported in [51] (𝜒 = 0.9). This is indication that the value of this coefficient depends on the considered nanoporous materials since the MPL microstructure studied in [51] is not obtained from FIB-SEM but via a numerical generation procedure. 1 summarizes the result obtained at the scale of the whole MPL FIB-SEM image (L 2 scale) for both block sizes. Some slight differences can be observed between the results for both block sizes but the results are quite close, indicating a quite weak dependence of the results with the block size at this scale.

MPL effective diffusion tensor at scale L2 Table

Consistently with the results at scale L 1 , the off-diagonal terms can be considered as negligible. A slight anisotropy is still present. However, based on the results shown in Table 1, the hypothesis of isotropy commonly considered in the literature is quite reasonable. 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.6 compared to the computation not considering the Knudsen diffusion.

MPL-GDM assembly

A rendering of the GDM-MPL assembly 3D image obtained via XCT is shown in Fig. 7. In order to determine the porosities and thicknesses of the assembly components, various limiting surfaces were defined as sketched in Fig. 7.

In Fig. 7, ℓ 1 is the thickness of the GDM region free of MPL: 3.

As indicated in Table 2, the MPL reduces the porosity of the assembly by 8% to 10% depending on the presence of cracks or not compared to the GDM alone. The MPL thickness is about 76 μm in average, which agrees with the apparent MPL thickness of SGL 22BB characterized in [10]. The thickness of the GDM alone is about 220 μm, but there is an important overlap between the MPL and the GDM.

GDM through plane effective diffusion coefficient

A digital 3D-image of the GDM alone is obtained from the image of the GDM-MPL assembly by changing the MPL voxels into gas voxels. Then, the through plane effective diffusion coefficient of the GDM alone is computed numerically solving the diffusion equation in the pore space using the same method as for the MPL effective diffusion tensor computation. Only the through plane component of the tensor is considered in what follows.

This yields a relative diffusion coefficient

𝐷 𝑒𝑓𝑓𝑥𝑥-𝐺𝐷𝑀 𝐷 𝑏𝑢𝑙𝑘 = 1 𝜏 ≈ 0.
8. This value is in good agreement with the literature for the porosity (𝜀 = 0.88) of the considered GDM. For example, a semi-analytical model based on the fiber organization was proposed in [50]. This model reads,

𝐷 𝑒𝑓𝑓 𝐷 𝑏𝑢𝑙𝑘 = 1 𝜏 = 1 6 ( 4.7-1.7𝜀 2.7-1.7𝜀 + 9 6.4-3.4𝜀 ) (16)
Applying Eq. ( 16) yields 1 𝜏 = 0.88, which is, however, a bit larger than the value obtained from our simulations. As discussed in [52], another analytical model commonly used for fibrous medium reads, (17) where 𝜁 and 𝜀 𝑝 are two parameters fitted from experiments. According to [52], 𝜀 𝑝 = 0.11, 𝜁 = 0.875 for application to GDM. This leads to 1 𝜏 = 0.89 that is again a bit greater than the value obtained from the simulations. This difference with the analytical model prediction is discussed in [27], where it is shown that the analytical or semi-analytical models overestimate the relative effective diffusion coefficient compared to the experiments. This is explained by the effect of the binder in the fibrous medium. The results presented in [52] can be described using a Bruggeman like relationship with an exponent -3, that is 𝜏 = 𝜀 -3 , which for 𝜀 = 0.88 gives 1 𝜏 = 0.68, lower than our value of 0.8. In this regard, it should be noticed that the GDMs considered in [52] are not from SIGRACET. They might contain more binder. Also, the GDL is compressed in [52], which is not the case of our GDL. This important point is further discussed in Section 5. Finally, as indicated in Table 2, the value of 0.8 from our simulations is in good agreement with the LBM simulations on reconstructed GDMs presented in [53], where a Bruggeman like relationship with an exponent 1.61 was found to represent the data obtained from the numerical simulations, i.e. 𝜏 = 𝜀 -1.61 . For our data, a best fit leads to a slightly greater exponent, 𝜏 = 𝜀 -1.74 .

𝐷 𝑒𝑓𝑓 𝐷 𝑏𝑢𝑙𝑘 = 1 𝜏 = ( 𝜀-𝜀 𝑝 1-𝜀 𝑝 ) 𝜁

GDM-MPL assembly through plane effective diffusion coefficient

In order to evaluate the impact of cracks in the MPL, two cases are considered: with and without cracks. To obtain a GDM-MPL assembly without cracks, the crack voxels are converted into MPL voxel according to the following procedure. Two planes are determined by computing the average position of the more advanced (least advanced respectively) MPL voxels in the through plane direction in the image. These planes are illustrated in Fig. 7b and correspond to the horizontal planes separated by thickness ℓ 2 + ℓ 3 . As illustrated in Fig. 7b, all crack voxels located between these two planes are converted into MPL voxel to obtain the GDM-MPL assembly without cracks. Results of simulations are reported in Table 2. The calculation of the through plane effective diffusion coefficient of the GDL with MPL in Table 2 first shows the impact of the MPL. The effective diffusion coefficient is reduced by a factor 0.82 with the MPL with cracks and 0.75 with the crack free MPL compared to the value for the GDM alone. A somewhat similar significant reduction in the effective diffusion coefficient was reported in [38] from simulations over reconstructed GDM and MPL with a decrease of 40 % in the effective diffusion coefficient. As discussed in [39], the decrease in the effective diffusion coefficient compared to the GDM alone depends on the ratio between the MPL effective diffusion coefficient and that for the GDM as well as the size of the GDM-MPL overlapping region. Since the thickness of the GDM-MPL assembly is different from the thickness of the GDM alone, it is more insightful to discuss the results in term of overall diffusion resistance [27]. The diffusion flux through a layer of thickness ℓ is expressed as,

𝐽 = ℰ𝐷 𝑒𝑓𝑓 ∆𝑐 ℓ = ℰ𝐷 𝑏𝑢𝑙𝑘 1 𝜏 ∆𝑐 ℓ = ∆𝑐 𝑅 ( 18
)
Where ∆𝑐 is the concentration difference across the layer. Eq. ( 18) defines the layer diffusion resistance as 𝑅 = ℓ𝜏 ℰ𝐷 𝑏𝑢𝑙𝑘 . Using the data in Table 2, the GDM-MPL assembly diffusion resistance is 60 % larger in the presence of cracks and 80 % in the absence of cracks compared to the GDM alone. The crack free GDM-MPL diffusion resistance is therefore about 10% larger than for the GDM-MPL assembly with cracks. This is a clear indication that the cracks facilitate the oxygen access to the catalyst layer (when they are free of liquid water). Table 3. Expression of diffusion resistances used in the resistance models sketched in Fig. 7.

Cracked MPL Uncracked MPL R 1 ℓ 1 𝜏 𝐺𝐷𝑀 ℰ 𝐺𝐷𝑀 𝐷 𝑏𝑢𝑙𝑘 ℓ 1 𝜏 𝐺𝐷𝑀 ℰ 𝐺𝐷𝑀 𝐷 𝑏𝑢𝑙𝑘 R 2 ℓ 2 𝜏 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 (1 -ℰ 𝑐 )ℰ 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 𝐷 𝑏𝑢𝑙𝑘 ℓ 2 𝜏 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 ℰ 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 𝐷 𝑏𝑢𝑙𝑘 R 3 ℓ 3 𝜏 𝑀𝑃𝐿 (1 -ℰ 𝑐 )ℰ 𝑀𝑃𝐿 𝐷 𝑏𝑢𝑙𝑘 ℓ 3 𝜏 𝑀𝑃𝐿 ℰ 𝑀𝑃𝐿 𝐷 𝑏𝑢𝑙𝑘 R 4 (ℓ 2 + ℓ 3 ) ℰ 𝑐 𝐷 𝑏𝑢𝑙𝑘 -
An interesting question is the impact of the overlapping region between the GDM and the MPL [39]. Let-us assume no overlap. Then from the definition of the diffusion resistance, the diffusion resistance of the assembly is simply the sum of the MPL and GDM resistances,

(ℓ 𝐺𝐷𝑀 +ℓ 𝑀𝑃𝐿 )𝜏 𝐺𝐷𝑀-𝑀𝑃𝐿 ℰ 𝐺𝐷𝑀-𝑀𝑃𝐿 = ℓ 𝐺𝐷𝑀 𝜏 𝐺𝐷𝑀 ℰ 𝐺𝐷𝑀 + ℓ 𝑀𝑃𝐿 𝜏 𝑀𝑃𝐿 ℰ 𝑀𝑃𝐿 (19)
where ℓ 𝐺𝐷𝑀 = ℓ 1 + ℓ 2 and ℓ 𝑀𝑃𝐿 = ℓ 2 + ℓ 3 . Using the data from Tables 1 and2 and the thickness values indicated in Fig. 7, this gives 1/𝜏 𝐺𝐷𝑀-𝑀𝑃𝐿 = 0.64 in the case of an uncracked MPL, which is greater than the computed value of 0.6 (Table 2, uncracked MPL). The total diffusive resistance of the assembly in the absence of overlap is greater than the resistance with overlap ( 𝑅 𝐺𝐷𝑀-𝑀𝑃𝐿-𝑛𝑜-𝑜𝑣𝑒𝑟𝑙𝑎𝑝 𝑅 𝐺𝐷𝑀 = 1.76 against 1.6 with the overlap without cracks in the MPL, Table 2), mainly due to the greater thickness of the assembly in the absence of overlap.

Referring to Fig. 7b, the consideration of the overlap zone in the resistance model leads to the following relationship,

(ℓ 1 +ℓ 2 +ℓ 3 )𝜏 𝐺𝐷𝑀-𝑀𝑃𝐿 ℰ 𝐺𝐷𝑀-𝑀𝑃𝐿 = ℓ 1 𝜏 𝐺𝐷𝑀 ℰ 𝐺𝐷𝑀 + ℓ 2 𝜏 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 ℰ 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 + ℓ 3 𝜏 𝑀𝑃𝐿 ℰ 𝑀𝑃𝐿 (20)
where 𝜏 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 is the tortuosity of the overlap zone. The latter can be estimated from the relationship,

𝜏 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 = 𝐷 𝑏𝑢𝑙𝑘 𝐷 𝑒𝑓𝑓-𝑜𝑣𝑒𝑟𝑙𝑎𝑝 = 𝐷 𝑏𝑢𝑙𝑘 𝐷 𝑒𝑓𝑓-𝑀𝑃𝐿 𝐷 𝑒𝑓𝑓-𝑀𝑃𝐿 𝐷 𝑒𝑓𝑓-𝑜𝑣𝑒𝑟𝑙𝑎𝑝 = 𝜏 𝑀𝑃𝐿 𝜏 𝐺𝐷𝑀 (21)
whereas the porosity in the overlap region is given by ℰ 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 = ℰ 𝐺𝐷𝑀 ℰ 𝑀𝑃𝐿 . Using the latter relationship and Eq.( 21) and applying Eq.( 20) leads to 1/𝜏 𝐺𝐷𝑀-𝑀𝑃𝐿 = 0.53. This is smaller than the computed value (0.6 in Table 2 for the case of the crack free MPL) but clearly illustrates the impact of the overlap on the diffusion process [54] from the comparison with the value 0.64 in the absence of overlap.

In the case of the MPL with cracks the relative through plane effective coefficient is larger as shown in Table 2. Considering the cracks in the resistance model is not difficult as sketched in Fig. 7b. This leads to the relationship,

(ℓ 1 +ℓ 2 +ℓ 3 )𝜏 𝐺𝐷𝑀-𝑀𝑃𝐿 ℰ 𝐺𝐷𝑀-𝑀𝑃𝐿 = ℓ 1 𝜏 𝐺𝐷𝑀 ℰ 𝐺𝐷𝑀 + 1 𝜀 𝑐 (ℓ 2 +ℓ 3 ) +(1-𝜀 𝑐 )( ℓ 2 𝜏 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 ℰ 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 + ℓ 3 𝜏 𝑀𝑃𝐿 ℰ 𝑀𝑃𝐿 ) -1 (22)
where 𝜀 𝑐 represents the crack volume fraction in the MPL. From image analysis, 𝜀 𝑐 ≈ 16%.

Applying Eq.( 22) yields 𝜏 𝐺𝐷𝑀-𝑀𝑃𝐿 = 0.67, which is in excellent agreement with the numerical simulation computed value 0.66 reported in Table 2.

Impact of GDM compression

The GDM-MPL assembly is compressed in a PEMFC whereas the 3D digital image of the GDM-MPL assembly obtained via XCT illustrated in Figs. 2 and 7 corresponds to an uncompressed assembly. The compression of the MPL can be expected to be negligible compared to the GDM. It has been shown that the compression has an impact on the GDL transport properties, e.g. [55], [56] and references therein. Thickness compression rates in the range [15% -47%] are for example considered in [54]. In this regard, it would be certainly interesting to obtain images of the compressed GDM-MPL assembly via XCT, as for instance in [27] or [56], and repeat the numerical simulations. However, since we do not have such images, insights on the impact of the compression can be gained by adapting Eq. ( 22) noticing that the property of interest here, i.e. the through plane diffusion coefficient, only depends on two parameters, the porosity and the tortuosity. The impact of compression on both parameters can indeed be relatively easily evaluated from relationships available in the literature. Considering that the MPL is not affected by the compression, the compression impacts the GDM thickness, porosity and tortuosity outside the overlap region where the GDM properties are assumed not affected by the compression. Following [57], the porosity of the GDM can be estimated from the relationship

𝜀 𝐺𝐷𝐿 = 𝜀 0 -𝛼 1-𝛼
, where 𝜀 0 is the porosity of the uncompressed GDM and 𝛼 is the GDM thickness compression rate defined as 𝛼 = ℓ 0 -ℓ ℓ 0

, where ℓ 0 is the thickness of the uncompressed medium and ℓ is that of the compressed medium. As discussed in Section 4.2.1, the GDM tortuosity can be described by a power law, namely 𝜏 𝐺𝐷𝐿 = 𝜀 𝐺𝐷𝐿 -1.74 . Accordingly, Eq. ( 22) can be expressed as

(ℓ 1 (1-𝛼)+ℓ 2 +ℓ 3 )𝜏 𝐺𝐷𝑀-𝑀𝑃𝐿 ℰ 𝐺𝐷𝑀-𝑀𝑃𝐿 = ℓ 1 (1 -𝛼) ( 𝜀 0 -𝛼 1-𝛼 ) -2.74 + 1 𝜀 𝑐 (ℓ 2 +ℓ 3 ) +(1-𝜀 𝑐 )( ℓ 2 𝜏 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 ℰ 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 + ℓ 3 𝜏 𝑀𝑃𝐿 ℰ 𝑀𝑃𝐿 ) -1 (23)
where

ℰ 𝐺𝐷𝑀-𝑀𝑃𝐿 = ℓ 1 (1-𝛼)( 𝜀 0 -𝛼 1-𝛼 )+ℓ 2 ℰ 𝑜𝑣𝑒𝑟𝑙𝑎𝑝 +ℓ 3 ℰ 𝑀𝑃𝐿 (ℓ 1 (1-𝛼)+ℓ 2 +ℓ 3 )
.

As illustrated in Fig. 8, the equivalent resistance model with compression shows that the impact of cracks on the assembly through plane diffusion coefficient compared to the case without cracks is similar at different compression ratio, a little bit more important at high compression. For a typical compression rate of 20%, the effective diffusion coefficient decreases by about 6%. Based on the results shown in Fig. 8, it can be concluded that the impact of the compression is comparable in magnitude to the impact of the cracks Fig. 8. GDM-MPL assembly through plane effective diffusion coefficient for different compression rates from equivalent resistance model (Eq.( 23)).

Conclusion

A multiscale approach was developed to study the gas diffusion properties of the GDM-MPL assembly from numerical simulations on 3D -digital images obtained via FIB-SEM (MPL) and X-ray computed tomography (GDM, GDM-MPL assembly) considering both Knudsen diffusion and the MPL cracks. The MPL matrix effective diffusion tensor was characterized and shown to be slightly anisotropic. The commonly used assumption of isotropy is therefore quite reasonable for the MPL sample size considered in the present study.

Knudsen diffusion is important in the MPL matrix. It was shown that the effective diffusion coefficients (diagonal terms of the tensor) decrease by 40% compared to the computation considering the Fick diffusion only.

The simulation procedure was compared against experimental data confirming, as suggested in [50], that Bruggeman model with the exponent 0.5 overestimates the effective diffusion coefficient.

Computation of the through plane effective diffusion coefficient of the GDM-MPL assembly

shows that the MPL has a quite significant impact on the GDM-MPL diffusion property compared to the GDM alone. Consistently with previous works, the MPL-GDM overlap region has a quite noticeable impact on the diffusion resistance of the assembly. Somewhat surprisingly, MPL cracks have a relatively weak impact on the assembly diffusive resistance, with a reduction by about 10%.

An equivalent resistance model was used to characterize the GDM-MPL assembly through plane diffusion coefficient. This model was used to explore the impact of the GDM compression on the diffusion property of the GDM-MPL assembly and it was found that the impact of the compression is on the same order of magnitude as the impact of the cracks in the MPL.

Nevertheless, since the focus was on the diffusion properties only, further studies are obviously necessary to assess from modelling the impact of MPL cracks on PEMFC performances. Other transport phenomena at play in the PEMFC must be considered as well as the coupling between the various layers in the PEMFC. More generally, it should be clear that modelling a full MEA implies to consider many other phenomena, such as gas diffusion in CL, its compression, membrane deformation, water management and so on. In this respect, the work illustrates the potential of multiscale imaging of porous materials combined to numerical simulations for characterizing the various porous components and ultimately developing comprehensive simulations in order to both better understand the various phenomena occurring in a PEMFC and offer insights for improving the performance and durability of PEMFCs.

Introduction

Polymer Electrolyte Membrane Fuel Cells (PEMFC) have emerged as an alternative power source to combustion engines and batteries. However, improvements in terms of durability, performance and cost effectiveness are still necessary. A PEM fuel cell is an assembly of various layers among which the cathode catalyst layer (CCL) is of central importance since this is the place of the oxygen reduction reaction (ORR), which is the most limiting reaction.

A CCL is a complex layer in which carbon black, Pt particles and ionomer are combined to form a nano-porous material in which electrons, protons, water and oxygen are transported.

To optimize the use of the costly platinum in the CCL and more generally the CCL performances, it is crucial to better understand and reduce as much as possible the different resistances that limit the efficiency of the various transport processes occurring in the CCL.

In particular, this concerns the ionomer which is a pathway for the ionic and gaseous reactants toward the platinum catalyst surfaces and acts also as a binder for the platinum loaded carbon particles. As a result of the CCL fabrication process from a dispersion of the above ingredients in a solvent [1], the carbon particles form agglomerates in the CCL [2][3].

A low carbon agglomerate coverage by the ionomer or a poor connectivity of the ionomer in the CCL is detrimental to the performance [4]. The ionic transport is affected, as well as the wettability properties of the CCL and thus the impact of liquid water in the CCL [5].

Moreover, the local accumulation of ionomer introduces an additional gas diffusion resistance [1]. Therefore, it is essential to understand the ionomer distribution in the CCL, the latter being considered as the most structurally complex component of the CCL [6,7].

The knowledge on the ionomer distribution in the CCL has significantly evolved over the years. In some earlier modelling studies, the pore space in the CCL is assumed completely filled with the ionomer [8,9]. However, experimental studies [10,11] at about the same time indicated that this assumption was questionable. Since these pioneering studies, several experimental works have been developed to better understand the structure of the CCL and the distribution of the ionomer in the CCL. It is now well established that the ionomer does not fully fill the pores in the CCL. Investigations based on cryogenic scanning electron microscopy (cryo-SEM) shows that the ionomer structure in the catalytic layer depends on the water/alcohol ratio (W/A) in the solvent [12,13], i.e. on the composition of the ink. The ink refers to the colloidal dispersion used to fabricate the CCL [1]. For a high W/A ratio, more ionomer is free in the solvent, while for a low W/A ratio, the free ionomer fraction is lesser. Free ionomer refers to ionomer not in the coating of the Pt/C particles. Thus, as discussed in [12,14,15], not all the ionomer present in the ink ends as ionomer coating the Pt/C phase. Information about free ionomer is also reported in [16,17,18], where it is also concluded that only a fraction of the ionomer initially present in the ink is coated on the Pt/C support. The remaining of the ionomer is deposited on the surfaces onto which the CCL is in contact during the fabrication process [1] or remains in the catalyst layer. The values of the corresponding fractions are not well known and depend on the carbon type and ionomer carbon weight ratio (denoted by I/C). In [16], it is reported based on electron tomography analysis that between 60 to 76% of the ionomer is in the coating, depending on the ionomercarbon ratio. In [19], where the investigation is based on nano-scale resolution X-ray computed tomography, it is considered that only 50% of the ionomer is in the coating of the carbon support. Both studies concern the case of Vulcan XC72 carbon black. The study in [20] for porous Ketjenblack using Cryo-electron tomography analysis shows that the segmented ionomer quantity depends on the ink preparation process. Also, the studies in [12] or [21] emphasize the importance of ionomer coverage on the electrocatalyst as a key parameter governing the properties of the catalyst layer, where the coverage is the surface fraction of the Pt/C phase surface coated with ionomer. This is in contrast with several numerical studies where, as mentioned later in the article, a full coverage is considered. The fact that the coverage is generally significantly less than 100% has been confirmed by other studies. The value of the coverage depends on the catalyst support [22]. For instance, a coverage of 75% on Ketjenblack is reported in [19] for an I/C ratio of 0.6 while values of 50 % and 80 % on Vulcan black are reported in [16] for an I/C ratio of 0.2 and 0.5 respectively.

Even for higher I/C ratios the coverage is significantly less than 100 % (0.76 for the I/C ratio of 1.2 in [20] for example based on cryo-electron tomography). The carbon particle agglomeration process also depends on the W/A ratio. Larger agglomerates are formed at high (W/A) ratio whereas the carbon agglomerates are smaller in the opposite case [12].

Contrary to an assumption often made in numerical models (some of them are briefly discussed below), the experimental studies have also shown that the thickness of the ionomer coating on the carbon surface is not uniform. In [16], the ionomer film thickness is reported to vary between 5 nm and 15 nm with an average thickness of about 7 nm. Similar values are reported in [23] from AFM measurements. In [19], the average film thickness is found to vary with the I/C ratio and a wider distribution with thicknesses as high as 30 nm or lower than 5 nm is reported. In [20], the film thickness on the Pt/C phase surface varies between about 2 nm and 20 nm with an asymmetric distribution skewed to the left, i.e. with a majority of thickness data in the range [2-6 nm]. Values as low as 2 to 3 nm are also reported in [24] from an analysis of small-angle neutron scattering profiles.

In addition to experimental characterizations, another type of approach consists in studying CCL properties or the transport phenomena occurring in a CCL by reconstructing the CCL microstructure via computational reconstruction methods [25]. In [26], the CCL is reconstructed by placing randomly carbon spheres of a given radius in the computational domain. An ionomer film coats the spheres. The ionomer film thickness is uniform and the coverage is 100%. Thus, there is no fraction of carbon particle surface in direct contact with the gas phase in the pore space. A similar assumption regarding the ionomer distribution is made in [27] where a reconstruction method based on a different technique for distributing the spherical particles in the computational domain is presented. Here again the ionomer is present as a film of uniform thickness covering all the carbon particle external surface. As mentioned above, considering an ionomer film uniform thickness and a full coverage is not consistent with the experimental observations. In more recent studies, i.e. [28,29,30], a procedure for generating a coating of variable thickness and partial coverage is implemented.

However, information on the obtained film thickness distribution and coverage in the reconstructed CCL images is not provided.

The method presented in [27] is used in a more recent study [31] where all the ionomer corresponding to a given I/C ratio is distributed in voxels adjacent to the solid phase (the spherical particles). The voxels adjacent to more solid voxels are filled with ionomer preferentially. This can be regarded as a preferential filling of the concave zones at the Pt/C phase surface as discussed later in the article. As a result, the ionomer film on the carbon particles is not of uniform thickness. However, no information on the resulting film thickness distribution is provided. Also, based on the results presented later in the present paper, it can be considered that the volume of ionomer in the macropores, i.e. the pores in between the carbon particles, is overestimated in this approach.

A different reconstruction method is used in [32] and also in [33]. The interesting idea is to perform the reconstruction by mimicking the CCL experimental fabrication procedure [1].

However, it can be noted that the size of 200 nm x 100 nm x 100 nm of the reconstructed sample in [32] is quite small compared to the CCL thickness (6 µm in the case of the CCL considered in the present study). The fabrication process seems also very roughly mimicked.

For instance, the drying step, which is a key step in the CCL fabrication process [1], is not at all considered. Also, the image of the reconstructed structure seems to indicate a full coverage and a too large presence of ionomer in the pores. However, it can be noted that the ionomer is not distributed as a thin film uniformly covering the surface of the Pt/C structure, which is more realistic that the uniform film thickness assumption considered in [26,27] for instance. In [33], simulations are performed over a significantly larger and thus presumably more representative reconstructed sample of size 4000 nm x 500 nm x 500 nm with a voxel size of 5 nm. However, as discussed later in the present paper, such a resolution is too low compared to the expected ionomer film thickness, which can be as low as about 2 nm [24].

The idea of mimicking the fabrication process to reconstruct the CCL is also present in a series of works where various assumptions regarding the ionomer distribution are considered [4,[34][35][36][37]. In [4], the CCL is reconstructed by first generating aggregates of spherical carbon particles. Then these aggregates are positioned in the computational domain to form the CCL. Regarding the ionomer, two cases are distinguished: the case of a uniform coating of the Pt/carbon phase by the ionomer and the case of a non-uniform coating. The latter is related to the drying step in the CCL fabrication process [1]. It is argued that the ionomer adheres to the Pt/carbon phase with a certain curvature resulting from the impact of capillary effects during the drying step. It is thus considered that the ionomer coating structure can be reproduced by a curvature model. However, no information on the corresponding curvaturebased ionomer distribution algorithm is given. Nevertheless, the study shows that the CCL effective transport properties are sensitive to the assumptions made to reconstruct the ionomer in the CCL. These two options (coating with uniform thickness or coating with uniform curvature) are also considered in [34,35] where the impact of the two options on the gas diffusion performance [34] and the ionomer conductivity [35] is further investigated. An attempt for modelling the drying step is proposed in [35] using the same concept as in [4,34] of agglomerate packing made of particle aggregates. However, the impact of drying on the ionomer distribution is not considered and an ionomer film with a constant thickness is assumed. The curvature-based ionomer distribution model is used again in [36] in a simulation study of Ketjenblack aggregates. Simulations at the scale of the CCL with consideration of the electro-chemical reaction are presented in [37] with again two options for the coating ionomer, uniform thickness or uniform curvature. The simulations show that the results depend on the ionomer adhesion model.

Although simulations on numerically reconstructed microstructures can be insightful, study of a porous material combining imagery techniques and pore scale numerical simulations has become a widely used approach, referred to as Digital Rock Physics in reservoir engineering [38]. This type of approach is expected to allow consideration of more representative microstructures of the considered materials and is also used to study the PEMFC porous components, e.g. [39] and references therein. The approach is often based on the use of highresolution X-ray Computed Tomography (XRCT). However, for nanoporous material such as a CCL, focused ion beam-scanning electron microscopy (FIB-SEM) [40] is an attractive option since the voxel resolution can be of only a few nanometers (5 nm in the present study). Both techniques, i.e. XRCT and FIB-SEM, are used to study PEMFC porous components and can actually be combined as in the case of the GDL as the difference in pore size between the microporous layer and the fibrous diffusion medium is significant [39]. In the present study, FIB-SEM is used to obtain a digital image of the considered CCL. FIB-SEM has been used in previous studies on CCL. In [41], it is reported that the ionomer cannot be distinguished from the carbon particle agglomerates. To introduce the ionomer in the digital image, a thin ionomer film of uniform thickness is coated numerically. However, it is not perfectly clear whether the ionomer is coated on the segmented solid phase, thus in the pores, or within the segmented solid phase, i.e. along its surface but inside. It is likely that the second option is chosen. In this approach the ionomer film thickness is based on the ratio between the average agglomerate size and the ionomer film thickness reported in the literature. It can be noted that the resolution of 10 nm in the FIB-SEM milling is large compared to the expected minimum value of ionomer thin film (about 2 nm according to [24]). A higher resolution of 5 nm is used in [42]. The ionomer film thickness is uniform and four CCL structures are obtained and compared by considering four different values of the film thickness, namely zero, 5nm, 10 nm or 15 nm. In [43], a slice thickness of 20 nm is considered, which is again large compared to the expected ionomer film thickness. However, no attempt is made to reconstruct the ionomer distribution in the FIB-SEM image. This is also the case in [44] where Cryo FIB-SEM with a slice resolution of 6 nm is used to study the oxygen diffusivity from numerical simulations on the FIB-SEM image. In [45], where additional references on FIB-SEM in relation with catalyst layers can be found, the FIB-SEM image is obtained with a resolution of 2 nm. The segmented solid phase in the FIB-SEM image is assimilated to the Pt/C phase. An ionomer film of uniform thickness (~8 nm) is then distributed along the surface of the segmented solid phase. The ionomer voxels are therefore selected in the segmented pore space and not in the segmented solid phase. Also, it can be noted that the simulations presented in [45] are performed over small sub-volumes (100 x100 x 100 nm 3 ) of the imaged CCL.

From this review of the literature, three main conclusions emerged. Firstly, the ionomer morphological distribution, i.e. its spatial distribution within the CCL microstructure, is a key factor for understanding the transport phenomena in the CCL and to optimize its performance from numerical simulations at pore scale. Secondly, the ionomer is often distributed in the reconstructed CCLs, either from imaging techniques or numerical reconstruction techniques, under the form of an ionomer thin film of uniform thickness. This type of assumption is questionable and a more realistic ionomer distribution is highly desirable. Thirdly, the ionomer cannot be distinguished in the FIB-SEM images because of the lack of contrast between the ionomer and the carbon support since these have a similar mass density.

As a result, the objective of the present work is to propose a numerical approach to reconstruct the ionomer distribution in a CCL microstructure three-dimensional image obtained by FIB-SEM. The assumption of uniform coating ionomer film thickness is not considered and a coating of variable thickness is obtained. As in some of the abovementioned works, namely [4,[34][35][36][37], the ionomer distribution is considered to depend on curvature. However, contrary to [4,[34][35][36][37] where the curvature of the ionomer is assumed to be constant, the curvature considered in the present work is the local curvature of the solid surface which is not constant but varies from one place to the other over the solid phase surface. As in [4,[34][35][36][37], the motivation for considering the local curvature as a parameter controlling in part the coating ionomer distribution is related to the drying step in the CCL fabrication process [1]. Drying is controlled in part by capillary effects [46] leading to the formation of curved liquid -gas interfaces. The resulting invasion process by the gas phase is hierarchical with preferential invasion of the bigger pores. Then it is considered that the suspension during this invasion process tends to be trapped by capillarity in the concave regions at the surface of the Pt/C phase. Ionomer would deposit in these concave regions as the result of evaporation of the dispersion trapped in these cavities. The final result would be that the more concave the surface cavity, the more ionomer is deposited. These qualitative considerations thus introduce a correlation between the local curvature of the solid phase surface and the ionomer presence and the coating ionomer local thickness, which thus varies depending on the local curvature. Also, contrary to several previous works on CCL microstructure reconstruction, the ionomer is not considered to be present only as coating ionomer. In addition to coating ionomer, a fraction of the ionomer is considered to be present as filaments partially occupying the pores. This type of ionomer structures can be observed in TEM image as reported in the present article and also in [20] or [47]. Also, as evidenced by the experimental studies and contrary to the assumption considered in some numerically reconstructed CCL, the coverage of the carbon support by the ionomer is not full and can vary depending on the considered I/C ratio. The coverage is therefore an important parameter in the presented reconstruction method.

Finally, as regards the methodological aspects, the presented CCL reconstruction method can be regarded as hybrid [48] since it combines an imagery technique (FIB-SEM) with a numerical reconstruction approach to reconstruct in the image an element, the ionomer, which cannot be directly seen in the FIB-SEM segmented image.

The paper is organized as follows. The imaging of the CCL via FIB-SEM and TEM is presented in Section 2. Section 3 provides information of the pore size distribution and size distribution in the solid phase. In Section 4, arguments leading to the conclusion that the ionomer is not deteriorated by the FIB but should be considered as part of the segmented solid phase are presented. Section 5 is the main part of the article with a presentation of the ionomer reconstruction algorithm and illustrations of the obtained ionomer distribution. A method for distributing the platinum particle in the reconstructed CCL is presented in Section 6. A discussion is proposed in Section 7. Section 8 consists of the main conclusions of the study.

Characterization of the CCL by FIB-SEM, Transmission Electron Microscopy (TEM)

The CCL considered in the present study is a commercially available catalyst layer (TEC10E50E; Tanaka Kikinzoku Kogyo K.K., Japan) with a Pt content of 46.2 wt%, I/C ratio of 0.8 and Ketjenblack high surface area (HSA) carbon.

FIB-SEM

3D FIB-SEM data was acquired on a Zeiss Crossbeam 550. A volume of 10.95 µm by 13.25 µm by 5.24 µm was imaged with pixel size of 5 nm and inter-slices distance of 5 nm too, resulting in cubic voxels of 5 nm 3 . The FIB milling was done at 30kV with a current of 700pA, and electron column was operated at 1.2 kV and 1 nA. Images from backscattered and secondary detectors were acquired simultaneously.

Image processing before segmentation included slice-to-slice alignment, homogenization of contrast/brightness in all three directions, correction of curtaining effect and resampling in z direction to correct for non-homogeneous inter-slice distances during milling and imaging.

These steps were done thanks to homemade python scripts based on usual image processing python libraries (e.g. numpy, scikit-image…).

As sketched in Fig. 1a, segmentation was done thanks to machine learning based on random forest classifier with ilastik software [49] using the 'pixel classification' workflow. The annotation and training phase was done on a sub-volume of 200x200x200 pixels. Fig. 1b illustrates the obtained segmented 3D image. The volume faction of the void phase in the image is 0.55. This volume fraction, denoted by 𝜀, is referred to as the CCL macro-porosity.

This macro-porosity does not include the nanopores in the carbon particles [50] which cannot be seen in the image for the considered resolution. In what follows, the porosity of the carbon particles is referred to as the micro-porosity. Also, due to the lack of contrast in SEM images between the ionomer and the other components, the ionomer distribution remains unknown.

Reconstructing the ionomer in the image and especially the ionomer coating is the major objective of the present article. The value of 0.55 for the CCL macro-porosity is comparable to the measurements reported in [51] but greater than the value of 0.4 reported in [43], both for a Vulcan black CCL. It is slightly lower than the values of 0.57-0.59 reported in [52] for graphitized Ketjenblack (I/C = 1,) and comparable to the values reported in [53] for the same carbon support as in our CCL (~ 0.55 for a I/C ratio of 1). It is shown in [53] that the final porosity depends on the conditions in the drying step of the fabrication process. Thus, some variability in the final porosity can be expected from one CCL to the other even for the same composition of the ink.

Transmission Electron Microscopy (TEM)

In order to get insights into the distribution of the ionomer in the CCL, high resolution images of small regions in the CCL were obtained by high resolution transmission electron microscopy (HRTEM). As discussed and shown for instance in [16,20], TEM allows to visualize the ionomer and the Pt loaded carbon particles. As illustrated in Fig. 2, inspection of the TEM images suggests to distinguish two different structures of ionomer in the images. The first one is referred to as ionomer filaments. The ionomer filaments partially occupy the void space and form protrusions on the surface of the carbon phase. This type of ionomer structure is illustrated in Fig. 2a. The second one corresponds to ionomer coating the carbon phase. This type of ionomer structure is illustrated in Fig. 2b and is referred to as coating ionomer film. As also illustrated in Fig. 2b, the coating ionomer tends to accumulate in the concave regions (referred to as cavities in Fig. 2) of the carbon phase surface. This point will be considered further in the modelling section.

Size distributions in void and solid phases

The Maximum Sphere Inscription (MSI) algorithm available in the image analysis toolbox This algorithm gives for each voxel, the diameter of the maximum sphere inscribed either in the solid phase or the void phase to which the considered voxel belongs. Therefore, application of the algorithm to each voxel in the solid phase gives a distribution of inscribed sphere diameters, which is referred to as the size distribution in the solid phase. Similarly, application of the algorithm to each voxel in the void phase gives a distribution of inscribed sphere diameters, which is referred to as the pore size distribution.

Figure 3b shows the size distribution in the solid phase, with a dominant value of 45 nm, which is larger than the expected average size of a carbon particle known to be around 30 nm [1,20]. The obtained larger value is regarded as a first indication that ionomer is present in the segmented solid. However, this larger size could be due in part to the arrangement of carbon particles in aggregates. It is reported in previous studies [35,36,55] that the aggregates are typically formed by only two or three particles on average, which is consistent with the dominant value of 45 nm observed here.

The obtained pore size distribution is shown in Figure 3a. As can be seen, the pore diameters vary between 15 and 200 nm with an average value around 60 nm. This result agrees with previous studies [4,56]. The local porosity in the segmented image is computed over cubic averaging volumes paving the image. The obtained distribution of local porosity is shown in Fig. 3c as a function of the cubic averaging volume size. For the cubic averaging volume of 1500 nm, the local porosity varies only by ±7% from one cubic volume to the other. This indicates that a representative elementary volume (R.E.V.) can be identified as regards the porosity. This is in contrast with the case of the micro-porous layer (MPL) reported in [39],

for which a R.E.V. could not be identified.

Impact of focused ion beam milling on ionomer

A first question is to identify whether the ionomer is still present in the solid phase segmented from the FIB-SEM image acquisition procedure since the question arises as to whether the FIB milling deteriorates the ionomer. To address this question, the thickness of the CCL is expressed from the CCL composition and its porosity as described below. The thickness is first expressed as:

𝑙 = 𝑚 𝑡 𝜌 𝑡 (1)
where 𝑚 𝑡 is the CCL mass per unit surface area and 𝜌 𝑡 is the CCL apparent density.

The mass 𝑚 𝑡 can be expressed as

𝑚 𝑡 = 𝑚 𝑃𝑡 + (1-𝑃𝑡/𝐶) 𝑃𝑡/𝐶 𝑚 𝑃𝑡 + 𝐼/𝐶 (1-𝑃𝑡/𝐶) 𝑃𝑡/𝐶 𝑚 𝑃𝑡 (2)
where 𝑚 𝑃𝑡 (mg/cm 2 ) is the platinum loading (0.2 mg/cm 2 for the considered CCL), Pt/C is the platinum / carbon weight ratio (defined as the mass of platinum over the sum of the mass of carbon and platinum), I/C is the ionomer -carbon ratio (defined as the mass of ionomer per unit mass of carbon). The CCL apparent density can be expressed as:

𝜌 𝑡 = 𝜀 𝜌 𝑎𝑖𝑟 + (1 -𝜀)(𝜙 𝐼 𝜌 𝐼 + 𝜙 𝐶 𝜌 𝐶 + 𝜙 𝑃𝑡 𝜌 𝑃𝑡 ) (3)
where 𝜙 𝐼 , 𝜙 𝐶 and 𝜙 𝑃𝑡 are the volume fractions of ionomer, carbon and platinum relative to the volume of the solid phase respectively, 𝜌 𝐼 is the ionomer density (2 g/ cm 3 ), and 𝜌 𝑃𝑡 is the platinum density (21.45 g/cm 3 ). In Eq.( 3), 𝜌 𝐶 is the carbon apparent density. The latter is less than the carbon density since the carbon particles are porous and can be expressed as

𝜌 𝐶 = (1 -𝜀 𝑚 ) 𝜌 𝐶 𝑚
, where 𝜀 𝑚 is the micro-porosity, i.e. the porosity of the carbon particles, and 𝜌 𝐶 𝑚 is the density of non-porous carbon (1.8 g/cm 3 ).

The volume fractions can be expressed as a function of the apparent carbon density 𝜌 𝐶 , the ionomer carbon ratio 𝐼/𝐶 and the platinum carbon ratio 𝑃𝑡/𝐶 as: (6) For the values of the I/C ratio (0.8) and Pt/C ratio (0.462) of the considered CCL, using the above equations leads to the results shown in Fig. 4b, where the CCL thickness computed from Eq.( 1) is plotted as a function of the CCL macro-porosity and carbon particle porosity. The latter cannot be determined from the FIB-SEM image and is considered here as a parameter. For the measured CCL thickness of 6 µm (Fig. 4a), the results shown in Fig. 4b indicate that the CCL macro-porosity should be between 50% and 60%, depending on the value of the carbon micro-porosity. This is in good agreement with the value of the macroporosity of 0.55 determined from the segmented image (Fig. 1b). Furthermore, it can be seen that the macro-porosity of 0.55 is obtained for a micro-porosity 𝜀 𝑚 of about 0.25 in Fig. 4b. This value is in quite good agreement with the value reported in the literature for the HSA carbon particles. This value is obtained for Ketjenblack carbon from N 2 adsorption measurement in [56]. The same value is considered in [57]. In [24], a liquid water fraction of 23% is quantified inside the carbon particles. This fraction can be assimilated to the microporosity. The conclusion is therefore that the FIB milling does not deteriorate the ionomer. This conclusion is also supported by the following consideration. The macro-porosity without ionomer, denoted by 𝜀 ′ , can be expressed from the macro-porosity with ionomer as 𝜀 ′ = 𝜀 + 𝜙 𝐼 (1 -𝜀). From the latter equation, the CCL porosity without ionomer should be between 65% and 76%, which is a high value compared to the value of 0.55 obtained from the segmented image. In conclusion, it is considered in what follows that the segmented solid phase in the image includes the platinum loaded carbon particles and the ionomer. Also, it can be seen from Fig. 4c that the platinum volume fraction in the solid phase is negligible compared to the carbon and ionomer volume fractions. The next step is to identify the ionomer and thus to differentiate the ionomer from the carbon phase in the solid phase segmented from the FIB-SEM images.

𝜙 𝐼 = 𝐼/𝐶

Ionomer filaments and coating ionomer film

A step-by-step algorithm is developed to discriminate in the FIB-SEM image the ionomer from the carbon phase in the segmented solid phase. The algorithm relies on the observations made on the TEM images presented in section 2.2 (Fig. 2). As indicated in section 2.2, the TEM images (Fig. 2) suggest to distinguish the ionomer filaments (Fig. 2a) from the coating ionomer films (Fig. 2b). The former, i.e. the filaments, correspond to pieces of ionomer in the void space forming protrusions from the carbon phase surface. Also, as exemplified in Fig. 2b, the TEM images show that the coating ionomer tends to accumulate in the concave regions at the surface of the carbon phase.

Ionomer filaments

The algorithm proceeds by first identifying the ionomer filaments. To this end, the results obtained from the Maximum Sphere Inscription (MSI) algorithm applied to the segmented solid phase are used. As discussed in Section 3 and illustrated in Figs. 3b and5a, application of the MSI algorithm to the segmented solid phase gives an inscribed sphere diameter distribution in the range [10 nm -90 nm] with an average size of 45 nm. The latter is greater than the expected mean size of carbon particles (30 nm) [19], which, as mentioned in Section 3, is regarded as an additional confirmation that the ionomer is not deteriorated by the FIB.

Ionomer filaments are detected from the MSI distribution considering that the ionomer filaments are smaller than carbon particles. To this end, an upper thresholding size is introduced, which is assumed to correspond to the filament maximum size. A filament minimum size is also introduced. This minimum size is first limited by the resolution of 5 nm. In addition, to avoid considering isolated voxels, the ionomer filament minimum size is slightly increased by considering that ionomer filament voxels must belong to voxel groups of at least three solid connected voxels (two voxels are considered as connected if they share a common face). This minimum size is consistent with TEM images (Fig. 2a is an example) showing some filaments with an apparent diameter on the order of 10 nm. Figure 5 illustrates the ionomer filament voxels (in green) identified for various upper thresholding sizes in the MSI distribution, i.e. the solid voxels associated with inscribed spheres in the solid phase smaller than the considered size limit. The ionomer filament distributions obtained for the upper thresholding size of 18 nm and 25 nm are shown over a smaller sample in Fig. 5b and5c. Ionomer shapes similar to those observed by TEM (Fig. 2a) can be observed in both cases. Table 1 gives the volume fractions corresponding to various values of the upper thresholding size. As shown in Fig. 3b and indicated in Table 1, the majority of the solid phase voxels is in inscribed spheres larger than 25nm. This observation is consistent with the average carbon particle size of 30 nm considered in many previous studies [1,20,48,56,58]. It is reported in [16] that only 60 to 76% of the ionomer is under the form of coating ionomer, depending on the ionomer-tocarbon ratio. Therefore, the filament ionomer cannot exceed 25 to 40% of the ionomer volume fraction. For Ketjenblack, it is assumed that not all the ionomer is in the coating or the filaments. This complementary fraction is discussed further in Sections 5.3 and 7.

Therefore, the ionomer fraction in the filaments should be even less than the 25-40% fraction mentioned above. In addition, the presence of some smaller carbon particles (20-25nm) is not impossible [48]. For all these reasons, an upper thresholding size of 18 nm is considered as the most appropriate in what follows. Although slightly lower, this size is also consistent with the ionomer thickness upper size of (20-24 nm) reported in [20]. As indicated in Table 1, this leads to an ionomer filament volume fraction of 7%.

Coating ionomer film. Computation of the solid phase surface local curvature

In this section, the approach for reconstructing the coating ionomer film in the image is presented. As discussed for instance in [1,4,53], the coating of the carbon particles by the ionomer results from a rather complex drying process starting from a liquid dispersion containing the ionomer, Pt/C particles and a multi-component solvent. As discussed in the Introduction, it is considered that this drying process leads to the preferential accumulation of the ionomer in the cavities at the surface of the solid phase as it is also observed in the TEM images (Fig. 2b), see also [20,47]. In other words, it is assumed that the concave zones at the solid surface are preferential sites for the ionomer deposit. As a result, the coating model presented in this work is based on the computation of the local curvature along the surface of the segmented solid phase. In this approach all the voxels corresponding to ionomer filaments identified in the previous step are excluded from the segmented solid phase.

The computation of the solid surface local curvature begins with the definition of a phase indicator 𝜙 taking the value 1 in the solid phase and 0 in the void phase. In other words, after excluding the filament ionomer voxels, the solid voxels are labelled by 1, and the pore voxels The curvature is computed from the unit normal vector to the interface as,

𝐶 = ∇. 𝒏 (9)
where 𝒏 is the interface unit normal pointing out from the solid phase into the void phase.

The latter is computed from the interface indicator 𝛼 as,

𝒏 =

∇𝛼 ‖∇𝛼‖+∆ (10) where ∆ is a constant to avoid the division by 0 [60]. Note that a curvature is computed for each voxel with the procedure. However the zone of interest where the curvature varies is obviously the interfacial region between the solid and the void. The procedure is implemented in OpenFOAM [60]. Due to the size of the image, a sequential approach is used, similar to that presented in our previous work on the computation of the effective diffusion tensor of a microporous layer [39]. For the present study, a cubic block with a side length of 1500 nm is considered as a representative elementary volume (Fig. 3c). A smaller cubic sample with a side length of 250 nm is extracted from this block for visualizations and illustrations. Figure 6b shows the computed local curvature after application of a smoothing procedure used to limit resolution effects. This smoothing procedure is illustrated in Fig. 6a.

As sketched in Fig. 6a, the smoothed curvature is the curvature averaged over a sphere with a radius of 3 voxels. As can be seen from Fig. 6b, the curvature distribution is heterogeneous and cavities with high concave curvature are present. This is in contrast with the model presented in [36], which assumes an iso-curvature after ionomer coating. However, it can be observed in Fig. 6b that the curvature can be considered as relatively constant locally, i.e. over small regions of the surface. This suggests that the iso-curvature assumption can be relevant over sufficiently small regions.

Coating ionomer film. Coverage modeling

The next step is to identify interface voxels as ionomer interface voxels. If one considers that the ionomer deposit occurs preferentially in the zones of higher concave curvature, one can then identify solid phase voxels at the solid phase surface corresponding to coating ionomer by setting a threshold value in the curvature distribution. It can be noted that these considerations imply that the coating by the ionomer does not smooth out the cavities at the surface of the Pt/C phase. Thus, our approach relies on the assumption that the ionomer coated carbon phase surface has about the same shape as the underlying carbon phase surface. This is consistent with the expected small thickness of the ionomer coating (on the order of a few nm, as mentioned in the introduction). Note that the curvature is negative in the concave zones due to the chosen normal orientation (and thus the more concave the surface, the more negative is the local curvature). Thus, all the interface voxels with a lower curvature than the threshold are identified as belonging to the ionomer coating whereas all the interface voxels with a greater curvature are identified as carbon phase voxels. The fraction of the interface voxels identified as coating ionomer voxels defines the coverage. Note that this step does not require to introduce explicitly an ionomer film thickness. This step is only based on the consideration of the interface voxels.

The coating ionomer distribution at the surface of the segmented solid phase obtained with this approach varying the curvature threshold value is illustrated in Fig. 7. According to this approach, convex zones are coated by the ionomer only for the higher values of the curvature threshold, which correspond to high values of the coverage. This agrees with the ionomer distribution shown in [20], where convex zones are free of ionomer as well as with the reconstructed distributions shown in [3] or [35].

Quantitative results are shown in Fig. 8. The obtained value of the coverage over a selected CCL sub-image of 1500 x 1500 x 1500 nm 3 is plotted as a function of curvature threshold in Fig. 8a. Two cases are distinguished in Fig. 8a: the case where the ionomer filaments are considered to determine the coverage and the case where the ionomer filaments are not considered. When the ionomer filaments are considered only the interface voxels not in contact with filament voxels are considered in the identification of coating ionomer from the local curvature. It is recalled that the considered volume fraction of the filament ionomer relative to the total volume of ionomer is 7% based on the discussion presented in section 5.1. Therefore, the ionomer coverage corresponding to ionomer filaments is not surprisingly quite low in Fig. 8a. It approximately represents 5% of the solid surface voxels. As a result, the impact of the filaments on the coverage is only visible at low levels of coverage. For higher levels of coverage, the zones covered by the ionomer filaments are superimposed with those with coating ionomer and both cases become indistinguishable in Fig. 8a. The question then arises as to what is the coverage in the considered CCL. Our approach establishes a link between the local curvature of the segmented solid phase and the coverage but does not tell us what the local curvature threshold and thus the coverage should be. A somewhat naïve view is to consider that all the ionomer initially present in the ink should be in the filaments and the coating. Then consider the results plotted in Fig. 8b showing the variation of the volume fraction of the ionomer (considering or not the filaments) identified via our approach as a function of the coverage, i.e. the local curvature threshold, for two values of the ionomer film thickness, namely 5 nm and 7.5 nm, consistently with the values considered in previous studies [15,19,23,41], considering here that the ionomer film thickness is uniform. According to Fig. 8b, the coverage would be about 100% for the film thickness of 7.5 µm and greater than 100% ! for the film thickness of 5 µm if one assumes that all the ionomer is in the filaments and the coating (data for the ionomer fraction of 100% in Fig. 8b). This is not consistent with the experimental observations indicating that the coverage is only partial [16,20,22,61]. Another way of illustrating this point is to compute the apparent I/C ratio. The latter is defined as the mass of ionomer identified with our approach, i.e. the ionomer in the coating and the filaments, divided by the mass of carbon in the considered domain. The evolution of the apparent I/C ratio as a function of the coverage is depicted in Fig. 8b (right vertical axis). This ratio is less than the expected I/C ratio of 0.8 for the considered CCL (except for the particular case of the average film thickness of 7.5 µm and the coverage of 100%). The conclusion is therefore that the ionomer can be divided in three fractions: the filaments, the coating and the not detected or missing ionomer. Table 2. gives the distribution of ionomer between the three forms discussed in this paper, namely the filament ionomer, the coating ionomer and the not detected ionomer for the value of the I/C ratio of 0.8, an average ionomer coating film thickness of 5 nm and the filament filtering size of 18 nm for various values of the coverage.

Coverage

As can be seen from Table 2, the missing ionomer represents a quite significant fraction of the ionomer for the range of coverage ([0.5 -0.9]) expected from the literature, i.e. [16], [19, [20]. The missing ionomer and the fact that the coating ionomer represents only a fraction of the ionomer in the CCL are further discussed in Section 7 in relation with the literature and the results shown in Fig. 8.

Coating ionomer film. Ionomer film thickness distribution

An ionomer average thickness of about 5 nm for Ketjenblack support is reported in the literature [36]. However, the TEM image in Fig. 2b, see also [20], indicates that the coating ionomer film thickness is not uniform with more coating ionomer in zones of greater concave curvature at the carbon surface. The next step is therefore to consider the variability of the coating ionomer film thickness. This is important not only in the prospect of computing the proton conductivity but also the CCL wettability properties [5] for instance. The local thickness is correlated to the local curvature, i.e. the lower the curvature, the greater is the thickness (it is recalled that the curvature is negative in the concave region for the chosen unit normal vector orientation at the surface, thus the most negative values of the curvature correspond to the most concave regions). Due to resolution constraints, a discrete distribution of the coating ionomer thickness is considered. The thickness can take the following values: 2.5 nm, 5 nm, 7.5 nm and 10 nm with an average value of 5 nm. This range is consistent with the ionomer film thicknesses reported in [20] from cryo-electrontomography analysis. The consideration of the ultra-thin thickness of 2.5 nm is made possible by using a numerical artificial resolution consisting in subdividing the segmented FIB-SEM image in voxels of 2.5 x 2.5 x 2.5 nm 3 . It can be noted that this thickness is consistent with results reported in [24], where an ionomer film thickness as small as approximately 2.7 nm is identified from small-angle neutron scattering (SANS) experiments. The ionomer film thickness histogram in Fig. 9a shows the considered correlation between thickness and curvature whereas Fig. 9b shows 3D renderings of the obtained ionomer film thickness distribution in a sub-region of the CCL of 250 x 250 x 250 nm 3 . A different colour is assigned to each discrete thickness in Fig. 9b. The various 3D images in Fig. 9b correspond to different side views of the same sample. As can be seen, the thickest ionomer film (10 nm in thickness) is present in the more concave zones while the thinnest ionomer film (2.5 nm in thickness) is present only on convex surfaces. The obtained distribution is qualitatively in agreement with the segmented distribution obtained from cryo-electron tomography analysis in [20]. The results shown in Fig. 9b are for a coverage of 82%. The reconstruction algorithm can be of course used for other values of the coverage.

Fig. 10 shows the obtained ionomer network with (right column) and without (left column)

the Pt/C phase. It can be seen in Fig. 10 that the coating ionomer distribution is not homogeneous over the CCL but varies from one region to the other. Nevertheless, as illustrated in Fig. 10, the ionomer is well connected. This feature is further investigated by identifying the coating ionomer percolating cluster using the cluster labeling algorithm available in Porespy [54], where the percolating cluster is the cluster connecting two opposite sides of the considered domain. The coating ionomer percolating cluster is illustrated in Fig. 11 for four values of the coverage. The ratio of the percolating coating ionomer volume to the total volume of the coating ionomer, referred to as the percolating volume fraction, is also indicated in Fig. 11 for the four coverage values. As can be seen from Fig. 11, the percolating coating ionomer forms a poorly connected structure only for a quite low coverage (0.23 in the example shown in Fig. 11). For such a low value of the coverage, the coating ionomer is mostly distributed in the non-percolating clusters (the percolating volume fraction is only of 9%). However, for the expected values of the coverage in a CCL, i.e. greater than at least 0.6, it can be seen from Fig. 11 that the situation is completely the opposite. Almost all the coating ionomer is in the percolating cluster. This is in line with the experimental observations reported in [20]. In other words, the coating film ionomer forms a very wellconnected percolating pathway for the H+ transport within the CCL for the usual values of the coverage expected in a CCL.

Platinum particles

The volume of platinum in the CCL is quite small compared to the volume of ionomer or carbon as illustrated in Fig. 4c. However, platinum particles play a key role in the oxygen reduction reaction in the CCL. Also, the platinum particles have an impact on the local contact angle [5], and thus on the liquid water distribution in the CCL. Therefore, it is interesting to insert platinum particles in the reconstructed CCL images, notably in the prospect of performing numerical simulations of transport phenomena coupled to electrochemistry in the CCL [45].

To this end, platinum particles are added in the segmented image as a final step of the CCL microstructure reconstruction, i.e. after segmentation of the solid phase in filament ionomer, coating ionomer film and carbon phase. It is known that the major fraction of platinum particles is inside the carbon particle porosity, about 75% according to [3]. Since the reconstruction of the porosity inside the carbon particles is beyond the scope of the presented reconstruction method, only the 25 % remaining platinum particles present on the carbon phase external surface are considered. From TEM image analysis, it is known that the platinum particle size distribution can be approximated by a Gaussian law, with an average diameter around 4 nm and a variance of 3nm [45]. This leads to consider a truncated Gaussian distribution of Pt spherical particle sizes in the range [2 nm -10 nm]. From this information, the platinum particles are distributed randomly on the external surface of the carbon phase assuming that particle centres are level with the carbon phase surface. This means that about 50% of the volume of each platinum particle is embedded in the carbon phase at the carbon phase surface. Illustration of the obtained distribution on the carbon phase is shown in Fig. 12a whereas Fig. 12b shows the distribution together with the reconstructed filament and coating ionomer.

Note that due to the size of the Pt particles compared to the FIB-SEM image voxel size of 5 nm, the considered image sample is sub-divided in much smaller voxels of size 0.5 nm to distribute the Pt particles over the carbon phase surface.

As can be seen from the comparison between Fig12a and Fig. 12b, some platinum particles can be entirely coated by ionomer (in concave regions of the carbon phase surface) whereas other particles are not in contact at all with ionomer (convex regions). Intermediate cases, i.e. platinum particles in partial contact with ionomer are also visible in Fig. 12b. Thus, the thickness of the ionomer fully or partially coating the platinum particles varies, which is consistent with literature, e.g. [5].

Discussion

According to our approach, all the ionomer is not in the coating or the filaments (the filaments representing a small fraction compared to the coating according to Table 2). This is consistent with results reported in [16] or [20] showing that the coating ionomer cannot be expected to represent the full amount of ionomer, in the CCL. As reported in [16,20], the ionomer coverage depends on the ink composition. Different models of coverage, depending on the I/C ratio, are presented in the literature, e.g. [35]. For an I/C ratio of 0.8, numerical simulations of ionomer coating on reconstructed carbon particle aggregates leads to a coverage of approximately 0.82 for Ketjenblack [36] and 0.95 for Vulcan XC72 [35].

However, these values overestimate the coverage determined in experimental studies [16,20,22,62]. From [62], the fraction of the coating ionomer is 50 %. Then, if only 50% of the ionomer initially present in the ink is assumed to form the coating, the simulation results presented in [36] agrees reasonably well with the experimental data. As can be seen from Table 2, this value of 50% is consistent with our results. As reported in [36], the simulation results (with the 50% correction) can be expressed according to the logarithmic fitted law 𝜃 = 0.18 ln(𝐼/𝐶) + 0.74 relating the coverage, denoted by 𝜃 and the I/C ratio. A coverage of 0.78 for an I/C ratio of 0.7 is reported in [19], which agrees with the logarithmic fitted law from [36]. As illustrated in Fig. 8b, this result agrees well with our data (consistently with the average film thickness of 5 nm reported in [36] for Ketjenblack). Also, it can be seen that the agreement is less good with the simulations reported in [35] for the case of Vulcan XC72. This is attributed to the approach used in [35], which assumes a constant coating ionomer fraction of 50% of the total ionomer, which is not accurate since the coating ionomer volume fraction depends on the ionomer-carbon ratio and the ionomer coverage, as demonstrated in [16]. In [16], two ink compositions of I/C ratio of 0.5 and 0.2 respectively are considered. It is reported that 60% (76% respectively) of the ionomer present in the ink is available for the covering for the ink with I/C ratio of 0.5 (0.2 respectively). As can be seen from Fig. 8b, although the I/C ratio in our case is different (0.8), these values are in quite good agreement with our data for the average ionomer film thickness of 7.5 nm consistently with the average film thickness reported in [16].

Based on the results reported in Table 2 and in the literature, the fraction of the ionomer not in the coating and the filaments is therefore significant. This missing fraction in our reconstruction procedure, also referred to as the non-reconstructed ionomer, is further discussed now.

Fig. 13 HRTEM image of ionomer confined between carbon particles.

The structure of the CCL is classically described as resulting from the aggregation of carbon particle aggregates during the drying step [4]. In the ink, the ionomer is considered to coat the aggregates [62]. In [62], this fraction of the ionomer in the ink is referred to as the adsorbed ionomer whereas the complementary ionomer fraction in the solvent is referred as the free ionomer [18]. According to [18] or [62], the adsorbed ionomer in the ink evolves into the ionomer coating on the platinum loaded carbon particles whereas the free ionomers randomly deposit to form ionomer bridges (called deposited ionomers in [18] or [62]) during the drying step. This description is consistent with the presence of ionomer filaments, which would correspond to deposited ionomer according to the terminology used in [18], [62]. The deposited ionomer could also correspond indeed to bridges as illustrated in Fig. 13. These considerations and especially Fig. 13 provide hints as regards the non-reconstructed ionomer.

Consider Fig. 13. The vertical length of the ionomer region between the two dashed yellow lines is about 60 nm. This region is identified as solid phase in the segmentation procedure.

Since the maximum ionomer film thickness considered in our reconstruction procedure is 10 nm (Fig. 9a) only a fraction of the ionomer visible in Fig. 13 will be identified as ionomer.

The central region of vertical extension of at least 40 nm is considered as carbon phase in our reconstruction procedure. Of course, Fig. 13 is a two-dimensional image. This means that we have implicitly considered that the thickness of the ionomer region in the direction perpendicular to the image was at least greater than 20 nm. Assuming that the situation depicted in Fig. 13 exists also in three-dimensions, it can therefore be considered that a fraction of the non-reconstructed (or missing) ionomer is actually in the phase identified as the carbon phase in our reconstruction procedure. One may argue that some ionomer is in the nanoporosity of the HSA carbon particles. However, this option is generally discarded and it is considered that the ionomer cannot intrude into the very small pores of the HSA carbon particles [3,50]. As mentioned in the introduction, a fraction of the not detected ionomer could correspond to ionomer deposited on the surfaces onto which the CCL is in contact during the fabrication process [1]. Nevertheless, it is reasonable to consider that most of the non-detected ionomer correspond to ionomer in the CCL as illustrated in Fig. 13, i.e. in between the carbon particles. This is in line with the schematics reported in [17] where the ionomer distribution is split into a network corresponding to the coating ionomer and an ionomer network in between the coated carbon, which would therefore correspond to the nodetected ionomer. It is argued in [17] that the latter network does not form a percolating network in the case of Nafion. This would be an important feature since an objective of the ionomer reconstruction is to characterize the proton and the oxygen transfers from numerical simulations, e.g. [31,34] and references therein.

Since a fraction of the ionomer is not reconstructed in the images, the question arises indeed as to whether the non-reconstructed fraction plays a significant role or not in the proton or oxygen transfer. Thus, in the case of Nafion, based on the consideration reported in [17],

reconstructing only the ionomer coating could be sufficient to characterize the transport. On the other hand, the proton conductivity is generally underestimated in simulations based on reconstructed CCL or CCL sub-region [30]. Since, as discussed in the introduction, very often, only the coating ionomer is considered in the reconstruction procedure, this might be due to the impact of the transport within the ionomer localized elsewhere than in the coating, i.e. in particular the non-reconstructed ionomer in our approach. This question will be explored further in our next research by computing the effective proton conductivity from our reconstructed ionomer distribution

Conclusion

In this work, a real cathode catalyst layer structure was successfully imaged using FIB-SEM.

Data analysis using the macro-porosity determined from the segmented FIB-SEM image led to the conclusion that the ionomer is included in the segmented solid phase. However, the ionomer cannot be directly identified in the three-dimensional segmented image of the CCL microstructure obtained by FIB-SEM. To overcome this limitation, a computational method was presented for studying the ionomer distribution in the CCL and reconstructing the ionomer distribution in the three-dimensional image. The method is based on the consideration of two forms of ionomer structure in the image: the ionomer distributed under the form of ionomer filaments partially occupying the pore space and the ionomer thin films coating the platinum loaded carbon phase. The ionomer filaments were identified by using an upper bound filter in the size distribution in the solid phase identified from the Maximum Sphere Inscription (MSI) algorithm whereas the reconstruction of the coating ionomer film was performed from the computation of the segmented solid phase surface curvature and the consideration that the concave regions of the surface are preferential sites at the surface for ionomer accumulation.

The method allows reconstructing the ionomer coating and the ionomer filaments over the full size of the CCL FIB-SEM image, i.e. over representative images of the CCL microstructure. It was also shown how platinum particles, which are not directly visible in the segmented FIB-SEM images, can be added computationally at the surface of the carbon phase. Thus, a comprehensive representation of the full catalyst structure is obtained, which includes ionomer, carbon, and platinum particles. The reconstructed coating ionomer is found to form a well-connected percolating network over the CCL thickness whereas the ionomer filaments are likely to have an impact on the gas transport as they can occlude, at least partially, pores, as well as on the proton conductivity. However, the study indicates that not all the ionomer is reconstructed. The missing fraction is considered to primarily correspond to ionomer trapped in between the carbon particles, which could correspond to the ionomer bridges mentioned in the literature [18].

In spite of this limitation, due in part to the considered image resolution, our results therefore pave the way towards the computations of the CCL effective transport properties, such as the oxygen effective diffusion tensor or the CCL proton conductivity, from numerical simulations of the considered transport phenomena on three-dimensional digital images obtained via high resolution imagery techniques.

Chapter 4

Pore-scale simulations of oxygen effective diffusion and proton effective conductivity tensors from FIB-SEM three-dimensional images of the cathode catalyst layer of proton exchange membrane fuel cells

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 threedimensional 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 the 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 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

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 in order to reduce development cycles via systematic studies of the design and optimization of fuel cell systems [4][5][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], finite volume [12] or finite element [13,14] methods, or image-based methods, such as the full-morphology (FM) method [14,15,16], are performed over threedimensional digital images of component microstructures. Pore network models have also been frequently used, e.g. [14,16,17,18] 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 performed using macro-homogeneous models [4,5,19]. This type of models involves effective parameters, such as the permeability tensor or the effective diffusion tensor, e.g. [20], 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,14,[21][22][23][24][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,21,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 images numerically according to a numerical reconstruction method [27]. Illustrations of this approach for the CCL can be found in [3], [21][22][23][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) [20,26,[30][31][32][33][34], see [35] 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%, a 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 pore space, this is of course essential in order to characterize the proton transport across the CCL. For this reason, as detailed in [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 [36], since it combines an imagery technique (FIB-SEM) with a numerical approach to reconstruct in the image the ionomer, 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 [37] to solve the oxygen transport equation in the gas phase and the proton transport equation in the ionomer.

The paper is organized as follows. Images from [26] are presented in Section 2 together with a summary of the method developed in [26] to reconstruct the various solid components, namely the ionomer, the platinum particles 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 consists of the main conclusions of the study.

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 3 . The segmented image is illustrated in Fig. 1. The thickness of the CCL sample shown in Fig. 1 is comparable to the CCL thickness (~ 6 µm). Details on the image processing used to obtain the segmented image can be found in [26]. As mentioned above, the ionomer cannot be distinguished from the carbon phase in the segmented image. Also, the platinum particles 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 particles are reconstructed in the segmented image using numerical techniques. The method for distributed the Pt particles 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]. Two types of ionomer structures, namely the ionomer filaments and the ionomer coating, 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 [38]. 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 [39][40][41][42].

The obtained ionomer distribution is illustrated in Fig. 2 over a sub-region of the CCL of 250

x 250 x 250 nm 3 . Since the voxel size in the original image is of 5 x 5 x 5 nm 3 , these voxels are subdivided into voxels of size 2.5 x 2.5 x 2.5 nm 3 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 consistently with experimental studies showing that the coverage can vary depending on the I/C ratio, e.g. [39] and references therein. Consistently with the experimental observations, [39][40][41], 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 ionomer reconstruction procedure does not lead to identify all the volume (or mass) of ionomer expected from the I/C ratio (0.8 in the case of the considered CCL). Some ionomer is missing. It is concluded in [26] that the fraction of the ionomer possibly destroyed during the milling steps of the FIB-SEM procedure is negligible. Futhermore, the CCL porosity of 0.55 obtained from the segmented FIB-SEM image is consistent with the consideration that the ionomer is present in the segmented solid phase [26]. As reported in [26], the ionomer filament identification is sensitive to the thresholding size used to separate the ionomer filaments from the remaining segmented solid phase. Thus, it is possible that the filament ionomer fraction is a bit underestimated. A fraction of the missing ionomer could also correspond to ionomer deposited on the membrane or the decal film [43],

i.e. the substrate onto which the CCL is formed. However, it is concluded in [26] that most of the missing ionomer is present within the segmented solid phase under the form of bridges in between the carbon agglomerates. This missing ionomer fraction is therefore considered as included in the solid phase identified as the carbon phase in our reconstruction procedure. In this respect, the results on the proton effective conductivity computation presented in Section 4 will also provide insights on the missing ionomer. For this reason, the segmented solid phase after reconstruction of the ionomer coating and the filaments can be decomposed into the reconstructed ionomer and the complementary solid phase (CSP). For the latter, although mostly composed of carbon, this terminology is preferred to "carbon phase" since some ionomer is supposedly present in the complementary solid phase.

It has to be noted that the Ketjen Black carbon particles in the CCL are nanoporous particles with a porosity of about 0.25 [44] and micropores less than a few nanometers [44]. This 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 [45],

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 [38] 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.

Furthermore, it is often considered that the primary pores are filled with water under usual PEMFC operating conditions [45]. 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.

Computation of oxygen and proton effective transport properties.

Oxygen and proton transport pore scale governing equations.

The method to compute the oxygen transport in the gas phase is the same as that presented in [20] for the case of a microporous layer (MPL). The oxygen transport by diffusion at the pore scale in the pore space is described by a steady-state Fickian diffusion equation

∇. (𝐷 𝑂 2 (𝒙)∇𝑐 𝑂 2 ) = 0 (1)
where 𝑐 𝑂 2 is the oxygen concentration, 𝒙 denotes a position vector and 𝐷 𝑂 2 is the local oxygen diffusion coefficient. Given the pore sizes in the CCL (of an average value of 60 nm as mentioned in Section 2), Knudsen diffusion [46] must be considered in conjunction with Fick's diffusion. As in several previous works, e.g. [9,47] and references therein, our approach is based on the Bosanquet relation [48] combining the molecule-wall and moleculemolecule interactions according to the formula,

𝐷 𝑂 2 (𝑥) = 1 𝐷 𝑂 2 -𝑏𝑢𝑙𝑘 -1 +𝐷 𝑂 2 -𝐾𝑛 -1 (2)
where 𝐷 𝑂 2 -𝑏𝑢𝑙𝑘 is the bulk diffusion coefficient and 𝐷 𝑂 2 -𝐾𝑛 is the Knudsen diffusion coefficient. The latter is expressed as

𝐷 𝑂 2 -𝐾𝑛 = 𝑑 𝑝 3 √ 8𝑅𝑇 𝜋𝑀 (3)
where R, M, T, d p are the ideal gas constant, molar mass, temperature and pore diameter, respectively. The latter is determined using a maximum sphere inscription algorithm [20].

As discussed for instance in [1], [44] or [49], the 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. [34], [49]. However as pointed out in [34], 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. As the result, the 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 focusses on the characterization of O 2 transport in gas phase.

The governing equation for the proton transport in the ionomer is also expressed via a diffusion equation [9], as

∇. (𝜎 𝑃 (𝒙)∇𝜑 𝐼 ) = 0 (4)
where 𝜎 𝑃 (𝑥) is the proton conductivity in the ionomer and 𝜑 𝐼 is the potential in the ionomer [22]. As reported for instance in [50], 𝜎 𝑃 strongly varies with the relative humidity, i.e the ionomer water uptake λ. Also, as pointed out in [1], specifying the conductivity in Eq.( 4) can be an issue since it has been shown that the conductivity of ionomer thin films differs from that of the membranes [51,52]. 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 [53,54]. As in several previous works, e.g. [9] and references therein, a simplified approach is considered, neglecting the impact of liquid water and 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. Nevertheless, the question of the possible impact of liquid water on the proton transport will be considered again in Section 4.

Upscaled transport governing equations

According to the volume averaging method [55], the upscaled oxygen diffusion transport governing equation under steady-state conditions reads,

∇. (𝜀𝐃 𝑂 2 𝑒𝑓𝑓 . ∇𝐶 𝑂 2 ) = 0 ( 5
)
Where 𝜀 is the CCL macro-porosity (it is recalled that the porosity inside the carbon particles is not considered), 𝐃 𝑂 2 𝑒𝑓𝑓 is the effective diffusion tensor and 𝐶 𝑂 2 is the intrinsic average concentration over the averaging volume, see [55] for more details.

Similarly, the upscaled proton transport governing equation is expressed as,

∇. (𝜀 𝐼 𝛔 𝑃𝑒𝑓𝑓 . ∇𝜙 𝐼 ) = 0 (6)
where 𝜀 𝐼 is the volume fraction of ionomer in the computational domain, 𝛔 𝑃𝑒𝑓𝑓 is the proton effective conductivity tensor and 𝜙 𝐼 is the intrinsic average potential.

On the definition of the effective transport coefficients

Although the formulation of Eqs.( 5) and ( 6) is fully consistent with the rigorous framework of the volume averaging method [55], it has to be noted that a different formulation is frequently used, notably in the PEMFC literature. The difference lies in the volume fraction term in Eq.( 5) and ( 6), which is often implicitly included in the effective tenor. In other words, Eqs.( 5) for instance (and similar equation for the proton transport) is often expressed in the literature as,

∇. (𝐃 𝑂 2𝑒𝑓𝑓 * . ∇𝐶 𝑂 2 ) = 0 (7)
Thus,

𝐃 𝑂 2𝑒𝑓𝑓 * = 𝜀𝐃 𝑂 2 𝑒𝑓𝑓 (8)
And

𝛔 𝑃𝑒𝑓𝑓 * = 𝜀 𝐼 𝛔 𝑃𝑒𝑓𝑓 (9)
This is important to notice for the comparison with data from literature.

Also, it is customary to discuss the components of these tensor in dimensionless form using the concept of tortuosity. For instance, the dimensionless through-plane component of oxygen effective diffusion tensor is expressed as

𝐷 𝑂 2𝑒𝑓𝑓𝑦𝑦 * 𝐷 𝑂 2 -𝑏𝑢𝑙𝑘 = 𝜀𝐷 𝑂 2𝑒𝑓𝑓𝑦𝑦 𝐷 𝑂 2 -𝑏𝑢𝑙𝑘 = 𝜀 𝜏 ( 10
)
where 𝐷 𝑂 2 -𝑏𝑢𝑙𝑘 is the molecular diffusion coefficient of oxygen for the considered pressure and temperature and 𝜏 is the tortuosity (written without the subscript yy for simplicity).

Similarly,

𝜎 𝑂 2𝑒𝑓𝑓𝑦𝑦 * 𝜎 𝑃 = 𝜀 𝐼 𝜎 𝑂 2𝑒𝑓𝑓𝑦𝑦 𝜎 𝑃 = 𝜀 𝐼 𝜏 (11)

Numerical simulation methodology (sequential approach)

As exemplified in [20] in the case of a MPL, the method for determining the effective transport properties consists in solving the pore-scale transport equations, Eqs. ( 1) and ( 4), and identifying the effective transport coefficients by expressing that the macroscopic fluxes obtained by averaging the pore-scale fluxes should be equal to their macroscopic mathematical expressions (this is detailed below). However, solving Eqs. ( 1) and ( 4) 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 [20], it is more efficient to perform the resolution of the transport problems at the CCL scale by a grid-block technique. The grid-block technique consists in partitioning the CCL matrix digital image in adjacent cubic blocks of equal size L 1 . The pore-scale governing equations, Eqs. ( 1) and ( 4), are solved for each grid block and effective transport properties are identified for each block.

Then, the macroscopic equations ( 5) and ( 6) 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 .

Effective transport properties over a block (scale L1).

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 [56] for lateral conditions. See [20] for the details on the Bamberger conditions. Zero flux conditions are imposed on the internal limiting boundaries via a penalization technique (see again [20]).

For the oxygen transport, Eq.( 1) is solved together with the above-discussed boundary conditions using the LaplacianFoam solver of the open source CFD software OpenFOAM [20].

Once the pore scale concentration field 𝑐 𝑂 2 (𝑥) is obtained, the average diffusive flux over the computational domain is computed as

〈𝒋 𝑂 2 〉 = 1 𝑉 ∫ (-𝐷 𝑂 2 (𝒙)∇𝑐 𝑂 2 ) 𝑑𝑉 𝑉 𝑓 ( 12
)
where 𝑉 is the volume of the computational domain and 𝑉 𝑓 is the volume occupied by the pore space in V.

Then, it is expressed that the average flux given by Eq. (7-12) must be equal to the macroscopic diffusion flux, namely

𝑱 𝑂 2 = -𝜀𝐃 𝑂 2 𝑒𝑓𝑓 . ∇𝐶 𝑂 2 = 〈𝒋 𝑂 2 〉 (13)
Consider a cubic grid block of side-length L 1 with a cartesian coordinate system and Dirichlet conditions 𝑐 𝑂 2 -1 and 𝑐 𝑂 2 -2 imposed in the x-direction, i.e. at the inlet and outlet faces of the black in the x -direction. Expressing the macroscopic gradient ∇𝐶 𝑂 2 as

∇𝐶 𝑂 2 = 𝑐 𝑂 2 -1 -𝑐 𝑂 2 -2 𝐿 1
𝒆 𝑥 , where 𝒆 𝑥 is the unit vector in the direction of x and H x is the dimension of the computational domain in the x direction and using Eq.(8-13) leads to the following relationships

𝐷 𝑂 2 𝑒𝑓𝑓𝑥𝑥 = - 𝐿 1 𝜀(𝑐 𝑂 2 -1 -𝑐 𝑂 2 -2 ) 〈𝑗 𝑂 2 𝑥 〉 (14)
𝐷 𝑂 2 𝑒𝑓𝑓𝑦𝑥 = - 𝐿 1 𝜀(𝑐 𝑂 2 -1 -𝑐 𝑂 2 -2 ) 〈𝑗 𝑂 2 𝑦 〉 (15)
𝐷 𝑂 2 𝑒𝑓𝑓𝑧𝑥 = - 𝐿 1 𝜀(𝑐 𝑂 2 -1 -𝑐 𝑂 2 -2 ) 〈𝑗 𝑂 2 𝑧 〉 (16)
Applying the same procedure in the y and z directions allows computing all the components of tensor 𝐃 𝑂 2 𝑒𝑓𝑓 ,

( 𝐷 𝑂 2 𝑒𝑓𝑓𝑥𝑥 𝐷 𝑂 2 𝑒𝑓𝑓𝑥𝑦 𝐷 𝑂 2 𝑒𝑓𝑓𝑥𝑧 𝐷 𝑂 2 𝑒𝑓𝑓𝑦𝑥 𝐷 𝑂 2 𝑒𝑓𝑓𝑦𝑦 𝐷 𝑂 2 𝑒𝑓𝑓𝑦𝑧 𝐷 𝑂 2 𝑒𝑓𝑓𝑧𝑥 𝐷 𝑂 2 𝑒𝑓𝑓𝑧𝑦 𝐷 𝑂 2 𝑒𝑓𝑓𝑧𝑧 ) 𝐿 1 = - 𝐿 1 𝜀 𝐿 1 (𝑐 𝑂 2 -1 -𝑐 𝑂 2 -2 ) ( 〈𝑗 𝑂 2 𝑥 〉 1 〈𝑗 𝑂 2 𝑥 〉 2 〈𝑗 𝑂 2 𝑥 〉 3 〈𝑗 𝑂 2 𝑦 〉 1 〈𝑗 𝑂 2 𝑦 〉 2 〈𝑗 𝑂 2 𝑦 〉 3 〈𝑗 𝑂 2 𝑧 〉 1 〈𝑗 𝑂 2 𝑧 〉 2 〈𝑗 𝑂 2 𝑧 〉 3 ) (17)
where indices 1, 2 and 3 correspond to the application of the inlet-outlet concentration difference (𝑐 𝑂 2 -1 -𝑐 𝑂 2 -2 ) in the x, y and z directions, respectively.

Effective transport properties at the larger scale (scale L2)

Once the effective tensor 𝐃 𝑒𝑓𝑓-𝐿 1 has been computed for each block of size L 1 according to the procedure described in Section 3.5, the diffusion problem at the larger scale L 2 is then formulated as

∇. (𝜀 𝐿 1 (𝒙)𝐃 𝑂 2 𝑒𝑓𝑓-𝐿 1 (𝒙). ∇𝐶 𝑂 2 ) = 0. (18)
with appropriate boundary conditions. In this approach, the effective diffusion properties determined for one block and is 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.5. Eq. ( 18) is solved numerically using the same boundary conditions as presented in Section 3.5 considering successively the directions x, y and z and the effective tensor components at scale L 2 are determined using relationships similar to Eq.( 17), namely

( 𝐷 𝑂 2 𝑒𝑓𝑓𝑥𝑥 𝐷 𝑂 2 𝑒𝑓𝑓𝑥𝑦 𝐷 𝑂 2 𝑒𝑓𝑓𝑥𝑧 𝐷 𝑂 2 𝑒𝑓𝑓𝑦𝑥 𝐷 𝑂 2 𝑒𝑓𝑓𝑦𝑦 𝐷 𝑂 2 𝑒𝑓𝑓𝑦𝑧 𝐷 𝑂 2 𝑒𝑓𝑓𝑧𝑥 𝐷 𝑂 2 𝑒𝑓𝑓𝑧𝑦 𝐷 𝑂 2 𝑒𝑓𝑓𝑧𝑧 ) 𝐿 2 = - 1 𝜀 𝑳 𝟐 (𝑐 𝑂 2 -1 -𝑐 𝑂 2 -2 ) ( 𝐻 𝑥 〈𝑗 𝑂 2 𝑥 〉 1 𝐻 𝑦 〈𝑗 𝑂 2 𝑥 〉 2 𝐻 𝑧 〈𝑗 𝑂 2 𝑥 〉 3 𝐻 𝑥 〈𝑗 𝑂 2 𝑦 〉 1 𝐻 𝑦 〈𝑗 𝑂 2 𝑦 〉 2 𝐻 𝑧 〈𝑗 𝑂 2 𝑦 〉 3 𝐻 𝑥 〈𝑗 𝑂 2 𝑧 〉 1 𝐻 𝑦 〈𝑗 𝑂 2 𝑧 〉 2 𝐻 𝑧 〈𝑗 𝑂 2 𝑧 〉 3 ) (19)
where 𝐻 𝑥 = 10.5 µm, 𝐻 𝑦 =12 µm, 𝐻 𝑧 =4.5 µm (Table 1) and

〈𝒋 𝑂 2 〉 = 1 𝑉 ∫ (-𝜀 𝐿 1 (𝒙)𝐃 𝑂 2 𝑒𝑓𝑓-𝐿 1 (𝒙). ∇𝐶 𝑂 2 )𝑑𝑉 𝑉 (20)
As for the first step (Section 3.5), Eq. ( 18) is solved using the LaplacianFoam solver in OpenFOAM.

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

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

O 2 1500 nm × 1500 nm × 1500 nm 7 L 1 × 3 L 1 × 8 L 1 (10.5 µm × 4.5 µm × 12 µm) 5 nm (L 1 ) L 1 /10 (L 2 ) H+ 250 nm × 250 nm × 250 nm 6 L 1 × 6 L 1 × 6 L 1 (1.5 µm × 1.5 µm × 1.5 µm) 0.833 nm (L 1 ) L 1 /10 (L 2 )
4.4 Oxygen and proton effective transport properties: results of simulations.

Oxygen effective diffusion tensor

Simulations are performed considering a total gas pressure of 2 bars and a temperature of 80°C. Fig. 3 presents the results obtained at scale L 1 (L 1 = 1500 nm, i.e. cubic grid blocks of size 1500 × 1500 × 1500 nm 3 ). Each tensor component is expressed in dimensionless form in term of tortuosity according to classical relationships of the form

𝐷 𝑂 2 𝑒𝑓𝑓𝑥𝑥-𝐿 1 𝐷 𝑂 2 -𝑏𝑢𝑙𝑘 = 1 𝜏
(as introduced in Section 3.3). We have plotted in Fig. 3 the components of the tensor 𝐃 𝑂 2𝑒𝑓𝑓 * = 𝜀𝐃 𝑂 2 𝑒𝑓𝑓 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 𝜀 𝜏 in Fig. 3.

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. Consistently with the results on the porosity discussed in [26], blocks of the considered size can therefore be considered as representative elementary volume (REV). .

The through-plane direction corresponds to the y axis. Somewhat surprisingly, the throughplane component 𝐷 𝑂 2 𝑒𝑓𝑓𝑦𝑦-𝐿 1 is in between the two components 𝐷 𝑂 2 𝑒𝑓𝑓𝑥𝑥-𝐿 1 and 𝐷 𝑂 2 𝑒𝑓𝑓𝑧𝑧-𝐿 1 in the in plane-directions. In this respect, it can be noted than the milling direction in the FIB-SEM procedure is the z-direction. Thus, 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 [20] 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. Fig. 3 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,22,57]. Fig. 4 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 𝜏 is often estimated using Bruggeman's relationship, 𝜏 = 𝜀 -0.5 [58]. As shown in Fig. 5, Bruggeman like expressions can be used to represent our simulation results. However, consistently with [58], 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 the Knudsen diffusion. As pointe out in [20], the exponent actually depends on the considered pressure, temperature and gas when the 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.

𝐷 𝑒𝑓𝑓𝑦𝑦 or 𝐷 𝑒𝑓𝑓𝑦𝑦 * with our notations and considered coordinate system. The comparison with previous works presented in Fig. 6 is therefore performed considering only 𝐷 𝑒𝑓𝑓𝑦𝑦 or more precisely 𝐷 𝑒𝑓𝑓𝑦𝑦 * .

The results of our simulations in Fig. 6 are for the scale L 1 (cubic blocks with 1500 nm sides).

As can be seen our results when the Knudsen diffusion is considered agree with those reported in He et al. [9] and Lange et al. [22] based on simulations on synthetic images generated according to the same generation technique but not with the results reported in Chen et al. [61], 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. [61] is much better when Knudsen diffusion is ignored whereas the agreement with the results reported in Lange et al. [22] become much less good under this condition. It is likely that this results from differences in the geometry and topology between our FIB-SEM images and the synthetic images used by Lange et al. [22]. 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 the Knudsen diffusion. The comparison with experimental results in Fig. 6 is not obvious since the agreement between two series of experimental measurements by the same group, i.e. Yu et al. [59,60] 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 [59] is different from the one used in [60]).

The fact that our results are similar to those reported in He et al. [9] and Lange et al. [22] when the 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 the O2 transport since the ionomer filaments are not considered in [9] or [22]. This is consistent with the quite low volume fraction of the ionomer filaments (≈ 1%). As in [9], [22] or [61], our simulations are performed considering a dry pore space. As noted in [9] and consistently 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 [59].

Tables 2 and3 summarizes the result obtained at the scale of the whole CCL FIB-SEM image (L 2 scale). Consistently 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, 𝐷 𝑒𝑓𝑓𝑥𝑥-𝐿 2 < 𝐷 𝑒𝑓𝑓𝑦𝑦-𝐿 2 < 𝐷 𝑒𝑓𝑓𝑧𝑧-𝐿 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.

Proton effective conductivity

As mentioned in section 3.7, a finer computational mesh is needed for computing the effective proton conductivity tensor. 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 0.833 nm of side for the proton transport computation. Due to the finer discretization, the block size L 1 is 250 × 250 × 250 nm 3 , thus smaller than for the oxygen transport computation. This corresponds to grid-block of 300 × 300 × 300 computational voxels. These blocks are partitioning a CCL sub-image of 1500 × 1500 × 1500 nm 3 . The considered sub-image of 1500 × 1500 × 1500 nm 3 was selected based on the results presented in Fig. 6 and Table 3 so as to be representative of the values reported in Table 3 as regards the oxygen effective diffusion tensor. × 250 × 250 nm 3 ). 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 ionomer volume fraction (𝜀 𝐼 ). As for the oxygen transport, the through plane diagonal component (σ 𝑃𝑒𝑓𝑓,𝑦𝑦 *

) tends to lie in between the diagonal components in the in-plane directions with the same trend as regards these The coverage, i.e. the fraction of the interfacial surface between the pore and the solid phase corresponding to ionomer, varies from one block to the other as indicated in Fig. 8a 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. 8 as a function of the I/C ratio (8b) and the block porosity (Fig. 8c). The I/C ratio considered here is 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 ionomer voxels and complementary solid phase voxels. Thus,

𝑉 𝐹𝐼𝐵-𝑆𝐸𝑀-𝑆𝑆𝑃 = 𝑉 𝐼 + 𝑉 𝐶𝑃 (21)
where 𝑉 𝐹𝐼𝐵-𝑆𝐸𝑀-𝑆𝑆𝑃 is the volume of the FIB-SEM segmented solid phase in the considered image, 𝑉 𝐼 is the volume of reconstructed ionomer in the image and 𝑉 𝐶𝑃 is the volume of the solid phase not identified as ionomer, i.e. the complementary solid phase. Then the apparent I/C ratio is defined as 𝐼/𝐶 = 𝜌 𝐼 𝑉 𝐼

𝜌 𝑐 𝑉 𝐶𝑃

, where 𝜌 𝐼 is the ionomer density (𝜌 𝐼 = 2 g/cm 3 ) and 𝜌 𝑐 is the carbon density (𝜌 𝑐 = 1.8 g/cm 3 ). This ratio is defined as apparent since, as explained previously, some ionomer, located elsewhere in the segmented solid phase than in the coating or the filaments, is missing in the reconstruction procedure. ) is illustrated in Fig. 7d. 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 × 1500 × 1500 nm 3 corresponding to the block assembly depicted in Fig. 8d), are summarized in Tables 4 and5.

Chapter 4

123 Table 4: Dimensionless components of CCL proton effective conductivity tensor 𝛔 𝑃𝑒𝑓𝑓 , i.e.

( The results presented in Tables 4 and5 confirm the trends discussed above for the scale L 1 , i.e. 𝜎 𝑒𝑓𝑓𝑥𝑥-𝐿 2 < 𝜎 𝑒𝑓𝑓𝑦𝑦-𝐿 2 < 𝜎 𝑒𝑓𝑓𝑧𝑧-𝐿 2 ,the off-diagonal components 𝜎 𝑒𝑓𝑓𝑧𝑦-𝐿 2 and 𝜎 𝑒𝑓𝑓𝑦𝑧-𝐿 2 are not negligible compared to the diagonal components contrary to the other offdiagonal components. From Table 4, the proton effective conductivity tensor can be considered as symmetrical as expected from the up-scaling theory [55].

A comparison with results from the literature [62][63][64][65][66][67][68] is proposed in Fig. 9 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. [64] and are close to the experimental data of Suzuki and al. [65] and Modestov et al. [67]. They are not too far from the experimental measurements of Du et al. [66]. However, the deviation is noticeable with the measurements of Yakolev et al. [63], Liu et al. [68] and Boyer et al. [62]. 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 (𝜎 𝑃 in Eq.( 4)) which can be different in thin ionomer films from the value for thicker films generally considered to obtain the factor ( 𝜀 𝐼

𝜏

) 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. 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 𝜀 𝐼 (as indicated by Figs. 7a and8a). Thus, a large variability of the coverage is not expected for a given 𝜀 𝐼 . 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. 9. As noted in [9], they are representative of the results obtained by simulations on synthetic images in other works [22,57]. 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 morphology 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 is not considered). However, as mentioned before and detailed in [26], a fraction of the ionomer present in the ink is missing with our reconstruction procedure. As reported [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 missing (for the whole image). It is considered in [26] that the missing ionomer is mostly with the carbon and the platinum in the segmented solid phase complementary to the reconstructed ionomer and filaments.

The sensitivity of the results to the presence of ionomer inside the complementary solid phase and of water (inside the primary pores) is explored as follows.

Sensitivity to proton transport in complementary solid phase.

Since a fraction of the ionomer is missing in the reconstruction process, it can be hypothesized, as discussed above, that at least a fraction of the missing ionomer is in the solid phase complementary to the reconstructed ionomer. With the considered assumption, this phase is actually a nanoporous medium made of carbon, platinum, ionomer and voids (possibly containing water). The proton transport governing equation in the complementary solid phase is expressed as

∇. (𝜀 𝐼𝐶 𝛔 𝑃𝑒𝑓𝑓-𝐶𝑃 . ∇𝜑 𝐼 ) = 0 (22)
Where 𝜀 𝐼𝐶 is the volume fraction of ionomer in the complementary solid phase and 𝛔 𝑃𝑒𝑓𝑓-𝐶𝑃 is the proton effective conductivity in the complementary solid phase. The latter is expressed

as 𝛔 𝑃𝑒𝑓𝑓-𝐶𝑃 = 𝜎 𝑃 𝜏 𝐶𝑃
, where 𝜏 𝐶𝑃 is the tortuosity in the complementary solid phase. For simplicity a homogeneous and isotropic conductivity tensor is assumed in the complementary solid phase. Eq.( 22) is solved together with Eq.( 4) in the reconstructed ionomer, namely

∇. (𝜎 𝑃 ∇𝜑 𝐼 ) = 0 (23)
The procedure consists in solving Eq.( 24) over the computational domain illustrated in The volume fraction 𝜀 𝐼𝐶 in Eq.( 22) is determined using Eq.(A.21) (derived in the Appendix).

For the considered block, 𝜀 𝐼𝐶 = 0.1. The tortuosity 𝜏 𝐶𝑃 in the complementary solid phase is varied and the effective through plane proton conductivity is obtained as described in Section As can be seen from Fig. 11 (square red symbols), the impact of the missing 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. 11 is by a factor of 1.75. To obtain the results shown in Fig. 11, we have implicitly assumed that the ionomer network in the complementary solid phase is percolating. However, it is argued in [69] that this network does not form a percolating network in the case of Nafion. On the other hand, it is considered in [45] 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.( 24) replacing 𝜀 𝐼𝐶 in Eq.( 22) by 𝜀 𝐼𝑊𝐶 = 𝜀 𝐼𝐶 + 𝜀 𝑊𝐶 where 𝜀 𝑊𝐶 is the volume fraction of primary pores in the complementary solid phase (𝜀 𝑊𝐶 = 0.25 as considered in [26]). Varying again the tortuosity leads to the results shown with blue symbols in Fig. 11. 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 considered 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 [58], namely 𝜏 𝐶𝑃 = 𝜀 𝐼𝑊𝐶 -0.5 , with 𝜀 𝐼𝑊𝐶 = 𝜀 𝐼𝐶 when only the missing ionomer is considered and 𝜀 𝐼𝑊𝐶 = 𝜀 𝐼𝐶 + 𝜀 𝑊𝐶 when water is considered with the missing ionomer. For the considered representative block, this leads to the results reported in Fig. 9.

In this figure, "correction 1" refers to the case where only the missing ionomer is considered whereas "correction 2" refers to the case where both the missing ionomer and water are considered. As can be seen, the consideration of the missing 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 [70,71]. Although the results shown in Fig. 9 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].

Impact of oxygen transport in ionomer

From the results obtained for the proton conductivity, 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 [34], namely,

∇. (𝐃 𝐼𝑂 2𝑒𝑓𝑓 * . ∇𝐶 𝐼𝑂 2 ) = 0 (25)
Where 𝐃 𝐼𝑂 2𝑒𝑓𝑓 * is the oxygen effective diffusion tensor in the ionomer and 𝐶 𝐼𝑂 2 is the oxygen concentration in the ionomer. The latter is related to the oxygen concentration in the gas phase by Henry's law [34,49], 𝐶 𝐼𝑂 2 = 𝐶 𝑂 2 𝑅𝑇𝐻 𝑒 , where 𝐻 𝑒 is the Henry's constant for oxygen dissolved in the ionomer, T is the temperature and R is the gas constant. 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 ionomer leads to define the global oxygen effective diffusion coefficient as

Eq.( 26) can be also expressed as

𝐷 𝐺-𝑂 2𝑒𝑓𝑓𝑦𝑦 * 𝐷 𝑂 2𝑒𝑓𝑓𝑦𝑦 * = 1 + 𝑅𝑇𝐻 𝑒 𝐷 𝐼𝑂 2-𝑏𝑢𝑙𝑘 𝐷 𝑂 2-𝑏𝑢𝑙𝑘 𝜏 𝑂 2 𝑦𝑦 𝜏 𝑃𝑦𝑦 𝜀 𝐼 𝜀 (27)
For representative values of the literature (𝐻 𝑒 = 5 ×10 -5 mol Pa -1 m 3 [72]), [34,72]) and the data on the ratios , is on the order of 10 -7 .

𝐷 𝐼𝑂 2-𝑏𝑢𝑙𝑘 𝐷 𝑂 2-𝑏𝑢𝑙𝑘 ≈ 10 -5
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.

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 proton conductivity and the effective oxygen diffusion in the secondary pores of a CCL 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.

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 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 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 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 Chapter 4

131 ionomer network discussed in [69]. 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 ionomomer 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.

Introduction

Polymer electrolyte membrane fuel cells (PEMFC) are the subject of a considerable interest in the context of the transition towards a decarbonized energy. Yet, improvements in terms of performance, durability and cost reduction are needed [1,2]. In this context, the cathode catalyst layer (CCL) is of central importance as critical component of PEMFC where the oxygen reduction reaction (ORR) takes place and costly platinum particles are used. This is also the place where water forms as the results of the ORR reaction and thus a key component as regards the crucial water management problem in PEMFCs [3,4]. The water management in the CCL and the impact of liquid water formation in the CCL have been the subject of many studies, e.g. [5][6][7][8][9][10][11] and references therein. These issues are intimately linked to the wettability properties of the CCL [12,13], and thus to the composition of the CCL. The latter is a nano-porous material obtained by drying a colloidal suspension made of platinum loaded carbon particles, ionomer and solvents [14]. The wettability of each constituent is quite different. Platinum is strongly hydrophilic whereas carbon is moderately hydrophilic with contact angles in the range [60 -90°] [13]. The ionomer can be hydrophilic or hydrophobic depending on its interaction with the platinum-carbon substrate [13]. As a result, as pointed out in [12], models for water management must take into consideration the heterogeneity of the wettability properties in the CCL. As discussed in [12] or [13], the wettability properties notably have an impact on the water retention curve (WRC), which is a key element of the macro-homogeneous models used for the simulations of two-phase flows in CCL, e.g. [15][16][17][18][19][20][21]. As exemplified in [13], the WRC can be obtained by combining the pore size distribution and the contact angle distribution via statistical approaches. However, a shortcoming of the approach presented in [13] is the oversimplification of the pore space geometry as a bundle of capillary tubes of circular cross-section. An alternative to obtain the WRC is to rely on numerical simulations on digital images of the microstructure of the porous medium, [22][23][24]. This can be achieved by direct simulations at the pore scale using numerical simulations techniques such as LBM [24] or volume of fluid methods (V.O.F) [26], or imagebased techniques such as the full morphology approach (FM) [22,23,27,28] but also via semi-integrated methods, such as pore network models [25,27,28,29]. An advantage of these methods is a much better consideration of the geometry and the topology of the pore space. However, using this type of approach implies to specify the contact angle distribution within the considered microstructure. This should be done locally in the case of the direct methods such as VOF or LBM, i.e. at a scale consistent with the numerical lattice or grid elements used with these methods or at the pore scale, i.e. on average in each pore, in the case of the PNM simulations. In this context the objective of the present article is to propose an approach to specify the contact angle distribution in a CCL microstructure at various scales of interest, typically at a scale smaller than the pore scale in the prospect of LBM or VOF simulations, or other direct numerical simulation techniques, e.g. [30], and at the pore scale, in the prospect of the PNM simulations.

As discussed in [31], three-dimensional digital images of CCL microstructures can be obtained via numerical generation techniques, e.g. [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46], or via experimental techniques such as nano-X-ray tomography [47][48][49] or focused ion beam-scanning electron microscopy (FIB-SEM) [30,[50][51][52][53][54][55][56][57][58]. In this article, our approach is illustrated using a 3D digital image of a CCL obtained by FIB-SEM but can be applied to microstructures obtained via other techniques. The input is in fact a voxelized three-dimensional digital image in which the various CCL constituents, namely carbon, platinum, ionomer and voids, are segmented. A motivation for considering an image obtained via an experimental technique, i.e. FIB-SEM, is that this type of image is presumably more representative of the real microstructure geometry and topology than images generated numerically. However, the FIB-SEM procedure only size in the solid phase obtained via a maximum sphere inscription algorithm available in the image analysis toolbox Porespy [59]. The coating ionomer corresponds to the fraction of the ionomer coating the carbon phase. The coating ionomer 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 [31].

The obtained ionomer distribution is illustrated in Fig. 1 over a sub-region of the CCL of 250

x 250 x 250 nm 3 . Since the voxel size in the original image is of 5 x 5 x 5 nm 3 , these voxels are subdivided into voxels of size 2.5 x 2.5 x 2.5 nm 3 in order to obtain a coating ionomer minimum film thickness of 2.5 nm. As illustrated in Fig. 1, the coating by the ionomer is not perfect and some regions of the carbon particle surface are not covered consistently with experimental studies showing that the coverage can vary depending on the I/C ratio, e.g. [60] and references therein. A coverage of 0.8 is considered in this work. Consistently with the experimental observations, [60][61][62], the convex zones of the carbon agglomerate surface are less likely to be covered than the concave zones. As discussed for instance in [13], two categories of pores are usually distinguished in a CCL: the primary pores and the secondary pores. The size of the primary pores is typically in the range (1-10 nm) in radius whereas the secondary pores are larger, 10-100 nm in radius.

Based on the voxel size of 5 nm of the FIB-SEM image, it is considered that the pore space identified in the FIB-SEM image corresponds to the secondary pores, also referred to as the macro-porosity. The primary pores are not visible and therefore not considered in the present article, which therefore focuses on the wettability properties in the pore space formed by the secondary pores.

Computation of contact angle distribution in CCL

Two type of methods are considered for computing the distribution of contact angles in the CCL image after the reconstruction of ionomer in the image. The first one is a fully deterministic approach considering the reconstructed distribution of the ionomer with the Pt particles. Its objective is the reconstruction of the local contact angle distribution, i.e.

considering a scale where the platinum particles can be considered individually. The second one is a hybrid method since it combines a statistical approach regarding the Pt particle distribution and the reconstructed ionomer in the image. The objective of the hybrid method is the reconstruction of effective contact angle distributions. The effective contact angle, which is a descriptor of the wettability properties on an average sense, is considered at different scales: at the scale of the voxel size considered for the ionomer coating reconstruction (2.5 x 2.5 x 2.5 nm 3 ), at the scale of sub-blocks of 250 x 250 x 250 nm 3 partitioning the CCL FIB-SEM image and at the pore scale after identification of pores in the pore space. In both approaches, i.e. deterministic or hybrid, the method to determine the contact angle as a function of ionomer coating thickness and the nature of the substrate, i.e. carbon or platinum, closely follows the approach presented in [13].

Local contact angle distribution

The aim is to determine the local contact angle distribution at a scale where the platinum particles can be identified individually in the image. The contact angle considered in this section is referred to as "local" contact angle since it corresponds to the smallest voxel size considered in this work, namely 0.5 x 0.5 x 0.5 nm 3 , and varies from one interfacial voxel to the other in the considered digital images of the CCL.

As explained in [31], the distribution of platinum particles illustrated in Fig. 2a for instance is obtained considering an artificial resolution of 0.5 nm. This means that the image of 250 x 250 x 250 nm 3 corresponding to the block illustrated in Fig. 2 obtained from FIB-SEM with voxels of 5 x 5 x 5 nm 3 is sub-divided into smaller voxels of size 0.5 x 0.5 x 0.5 nm 3 . Since this scale is consistent with the expected size of the platinum particles (Fig. 3), the contact angle is defined at this scale in this section. As illustrated in Fig. 2, the interfacial surface between the pore space and the solid phase is heterogeneous. The pore space is in contact over this interface with carbon free of ionomer (under the considered approximation that the impact of the missing ionomer is negligible, see Section 2 and Section 5) and platinum (shown in grey in Fig. 2b), platinum particles free of ionomer (in red in Fig. 2b), ionomer filaments (in green) and ionomer coating of variable thickness (Fig. 2c). The objective is to determine the contact angle for each type of interfacial voxel (an interfacial voxel is a solid voxel in contact with the pore space).

First, let us consider the region of the solid phase interfacial surface free of ionomer. The contact angle for the interfacial voxels corresponding to carbon is taken equal to 80° as in [63] or reported on carbon fiber in [64] whereas the contact angle is zero for the interfacial voxels corresponding to platinum, which is known to be strongly hydrophilic [12]. Regarding the ionomer coating interfacial voxels, we follow the approach presented in [13]. According to [13], the contact angle on the ionomer coating depends on the substrate and the thickness of the ionomer film. This dependence is expressed as follows. The cosine of local contact angle, 𝜃 𝑖𝑐 , at the coating ionomer surface is expressed as cos(𝜃 𝑖𝑐 ) = ∆ 𝐷𝐴 (cos(𝜃 𝐻𝐼 ) -cos(𝜃 𝐻𝑂 )) + cos (𝜃 𝐻𝑂 )

where 𝜃 𝐻𝐼 is the contact angle when the ionomer side chains are exposed to the pores (𝜃 𝐻𝐼 =30°) and 𝜃 𝐻𝑂 is the contact angle when the ionomer backbone is exposed to the pores (𝜃 𝐻𝑂 =120°). ∆ 𝐷𝐴 in Eq.( 1) is the degree of alignment (of the molecular structure in the ionomer), which depends on the substrate and the ionomer film thickness. The latter is expressed as

∆ 𝐷𝐴 = ∆ 𝐷𝐴,∞ + (∆ 𝐷𝐴,0 -∆ 𝐷𝐴,∞ )exp (- (𝑒-𝑒 0 ) 𝑒 0 ) ( 2
)
where 𝑒 is the coating ionomer film thickness (which in our case can take the discrete value indicated in Fig. 2c), 𝑒 0 is the thickness of the first layer of ionomer when it starts forming the coating (in the ink). The value of 𝑒 0 considered in [13] is 6.5 nm, which is clearly not consistent, viz. too large, with the film thickness distribution depicted in Fig. 2c. Thus, a smaller value, namely 2.5 nm, consistent with the value reported in [65], is used in what follows. In Eq. ( 2), ∆ 𝐷𝐴,∞ is the expected value of the degree of alignment for sufficiently thick ionomer film where a random mixture of side chains and backbone is expected at the ionomer surface. The value of ∆ 𝐷𝐴,∞ considered in [13] is 0.37. ∆ 𝐷𝐴,0 is the degree of alignment for a sufficiently thin film, for which the interactions with the substrate imposes the contact angle. In [13], ∆ 𝐷𝐴,0 = 0.25 for the carbon substrate (∆ 𝐷𝐴,𝐶 = 0.25) whereas ∆ 𝐷𝐴,0 = 0.

for platinum (∆ 𝐷𝐴,𝑃𝑡 = 0.). For platinum as a substrate and 𝑒 0 =2.5 nm, Eq.( 2) leads to the evolution of the degree of alignment shown in Fig. 4. With our choice of 𝑒 0 =2.5 nm, the degree of alignment goes from zero for the thinner ionomer film considered (of 2.5 nm in thickness) to almost 1 for the thickest film (10 nm). In other words, the impact of the degree of alignment is maximized (since it varies almost over the full range [0, 1]).

Application of the above formula for the various cases in our images leads to the contact angles values reported in Table 1. At the surface of the ionomer filaments, a "thick" ionomer value is considered, namely 90° (Table 1). As can be seen from Table 1, the contact angle varies over the interfacial surface area between zero (on platinum particles free of ionomer) to 120° (ionomer coating film of 2.5 nm in thickness on a platinum particle). The contact angle is maximum at the surface of the thinnest ionomer coating film on platinum as the result of the inversion of the substrate wetting behavior described in [13]. For thicker ionomer films, the difference in contact angle between the two substrates gradually vanishes with an increasing ionomer coating film thickness.

However, since the surface fraction of the pore space -solid phase interfacial surface corresponding to platinum particles is quite small ( ≈ 5 % ), the contact angle mostly varies in the range [90° -100°].

To distribute the contact angle in the image by using the contact values reported in Table 1, we have to determine the substrate associated with each ionomer coating voxel in contact with a pore (the ionomer film thickness is determined in the ionomer coating reconstruction procedure [31] as illustrated in Fig. 2c). To this end, two distances are determined: 𝑑 𝐶 and 𝑑 𝑃𝑡 .

Distance 𝑑 𝐶 is the distance between the center of the considered interfacial voxel and the closest carbon phase voxel whereas distance 𝑑 𝑃𝑡 is the distance between the center of the considered voxel and the closest platinum particle. Distance 𝑑 𝐶 (𝑑 𝑃𝑡 respectively) is determined from a distance map established considering only the coating ionomer and the carbon phase (the platinum particles respectively). Then the substrate considered to compute the contact angle for the considered coating ionomer interfacial voxel is the one corresponding to min (𝑑 𝐶 , 𝑑 𝑃𝑡 ). The result of this procedure is illustrated in Fig. 5. In Fig. 5b, the ionomer coating interfacial voxels are shown in orange where they are closer to platinum than to carbon and in blue in the opposite case. In the latter case, we consider that the contact angle is the one corresponding to the carbon substrate in Table 1 whereas the value corresponding to the platinum substrate is considered when the interfacial voxel is closer to platinum. This is an approximation since often it is likely that the contact angle, i.e. the structuration of the ionomer in the ionomer coating film, depends both on the closer substrate (platinum for instance) and the less close substrate (then carbon). Fig. 5c (shown also in Fig. 6 for better visibility) illustrates the local contact angle distribution so obtained. The largest value of the contact angle, i.e. 120° correspond to a thin ionomer coating film on platinum and the smallest value, i.e. 0° on platinum free of ionomer, correspond to rather small regions, consistently with the sparse distribution of the platinum particles. As a result, the contact angle varies mostly between its value on carbon free of ionomer and platinum, i.e. 80°, and the value on 2.5 nm thick ionomer coating film on carbon, i.e. 99°.

As illustrated in [13], the contact angle distribution is of particular interest in relation with the characterization of CCL global properties such as the water retention curve for instance, which, as mentioned in the introduction is an important input data for macroscopic models aiming at modelling water transfer in the CCL. The water retention curve can be obtained from direct simulation methods such as the Lattice Boltzmann method (LBM) [27] or the volume of fluid (VOF) methods [26] directly on digital images of the CCL such as the one illustrated in Fig. 6. However, this type of approach is computationally expensive [27] and simulating the very low capillary number flows needed to compute the capillary pressure nm in contact with the pores (and similar definitions for the three other ionomer coating thicknesses).

Note that contrary to section 3.1, the platinum particles are not explicitly considered in the decomposition (3). The platinum particles are considered according to a statistical approach.

The platinum particles are considered as spherical with the diameter distribution illustrated in Fig. 3. This is the distribution considered in [30]. This distribution is approximated by a truncated Gaussian law, with an average diameter of 4 nm and a standard deviation of 3nm.

As in [31], the platinum particles are considered as distributed randomly over the microstructure. It is also considered that 25 % of the platinum particles are on the carbon phase external surface (the complementary fraction of 75% is considered to be inside the carbon phase [30,31]). Let 𝑚 𝑝𝑡 denotes the platinum loading (𝑚 𝑝𝑡 =0.

where we have considered for simplicity that A without-ionomer ≈ A, since the coating follows the shape of the carbon phase. A without-ionomers represents the interfacial surface area of the solid phase in the absence of ionomers (i.e. when the ionomer voxels identified in [31] are discarded).

Let 𝑓 be the probability density of the platinum particle diameter, denoted by 𝑑 𝑝 . Then the mass of platinum particles over the considered surface can be expressed as,

𝑀 𝑝𝑡-𝑆-𝑗 = 𝑁 𝜌 𝑝𝑡 ∫ 𝑓(𝑑 𝑝 ) +∞ -∞ 𝜋 6 𝑑 𝑝 3 𝑑(𝑑 𝑝 ) (5)
where 𝜌 𝑝𝑡 is the platinum density and N is the number of platinum particles over the considered surface. Combining Eqs. ( 4) and (5)

which can be also expressed as,

𝑎 𝑝𝑡 = 3 8 𝑚 𝑝𝑡 𝜂 𝜌 𝑝𝑡 𝑙 𝐶𝐶𝐿 ∫ 𝑓(𝑑 𝑝 ) +∞ -∞ 𝑑 𝑝 2 𝑑(𝑑 𝑝 ) ∫ 𝑓(𝑑 𝑝 ) +∞ -∞ 𝑑 𝑝 3 𝑑(𝑑 𝑝 ) (9)
where 𝜂 = 𝐴 𝐿 3 is the CCL specific surface area, i.e. the interfacial area between the pore and the segmented solid phase per unit volume. Fig. 7 shows the variation of the "voxel" specific area as a function of L computed over cubic blocks of size 250 x 250 x 250 nm 3 paving a CCL image of size 1500 x 1500 x 1500 nm 3 , which corresponds to 216 blocks.

The voxel "specific surface area" is the dimensionless ratio defined by , where 𝐿 𝑣𝑜𝑥𝑒𝑙 is the voxel size (𝐿 𝑣𝑜𝑥𝑒𝑙 = 2.5 nm for the computations leading to the results shown in Fig. 7). The average specific surface area so obtained is equal to 2.84 × 10 7 m -1 .

It can be seen that the specific surface area varies from one block to the other. On the other hand, it is likely that the number of platinum particles over the external surface of the carbon phase in each block also varies, which could compensate the variation of the interfacial area.

In other words, it is hypothesized that the specific density of particles (number of particles per unit interfacial surface area) does not vary much from one block to the other. For the considered platinum particle diameter Gaussian distribution of average diameter 𝑑 ̅ 𝑝 (4 nm) and standard deviation 𝜎 𝑑 𝑝 (3 nm), Eq.( 9) can be expressed as, The complementary fraction is the surface fraction of carbon 𝑎 𝐶 = 1 -𝑎 𝑝𝑡 . As can been from Eq.( 8) or Eq.( 10), 𝑎 𝑝𝑡 does not depend on the particular surface considered, i.e. 𝐴 𝐶 , 𝐴 𝑖𝑐-2.5 , etc. This is because no spatial correlation between the presence of platinum particles and the ionomer film thickness is assumed, or more precisely between the presence density of platinum particles and the local curvature of the carbon phase surface.

𝑎 𝑝𝑡 =
Then, as in [13], an effective value for ∆ 𝐷𝐴,0 is obtained by linear interpolation as ∆ 𝐷𝐴,0 = ∆ 𝐷𝐴,𝐶 (1 -𝑎 𝑝𝑡 ) + ∆ 𝐷𝐴,𝑃𝑡 𝑎 𝑝𝑡 (11) Using this estimate in Eq.( 2) and then using Eq.( 1) allows determining the contact angle at the surface of the considered surface when this surface is coated with ionomer.

Fig.8. Variation of effective contact angle as a function of ionomer coating thickness

Fig. 8 illustrates the variation of the contact angle over the ionomer coating as a function of its thickness so determined. As can be seen, the contact angle variation over the ionomer coating is of about 10°, between 90° and 100°, as the result of the impact of the ionomer coating thickness on the contact angle.

On the platinum loaded carbon phase surface free of ionomer, the contact angle is estimated as cos (𝜃 𝐶 ) = cos (𝜃 𝐶𝑎 )(1 -𝑎 𝑝𝑡 ) + cos (𝜃 𝑃𝑡 )𝑎 𝑝𝑡 (12) where 𝜃 𝐶𝑎 is the contact angle on carbon free of platinum (𝜃 𝐶𝑎 = 80°) and 𝜃 𝑃𝑡 is the contact angle on platinum (𝜃 𝑃𝑡 =0°).

The filaments are considered as sufficiently thick for the contact angle to correspond to a thick film, thus ∆ 𝐷𝐴-𝑖𝑓 = ∆ 𝐷𝐴,∞ .

For the discrete distribution of ionomer coating film thickness illustrated in Fig. 2c, this leads to the values of contact angle summarized in As indicated in Table 2 and illustrated in Fig. 9, the local effective contact angle varies in the range [78° -100°] according to our computations. The pore sizes determined in [31] are mostly in the range [10 nm -120 nm] with a peak at 70 nm. Fig. 9 thus indicates that the local local effective contact angle fluctuates spatially within a pore. The range [78° -100°] is significantly broader than the one considered in [13] for the secondary pores (90.5° ± 1.5°). It could be considered that this is mainly because some carbon phase is free of ionomer in our case and because the minimum ionomer film thickness considered in [13] is significantly greater (6.5 nm instead of 2.5 nm in our work). However, the difference in the range of variations can be in fact considered as to be essentially a question of scale since, as mentioned before, this is the pore scale, i.e. domains of the size of the pores, which is considered in [13],

and not the smaller scale of the wettability heterogeneity illustrated in Fig. 9. To get further insights into this scale issue, we consider blocks of size L = 250 nm in the next section. Each color corresponds to a value of the contact angle. Carbon phase free of ionomer in gray, ionomer filament in green, ionomer coating in cyan (2.5 nm), yellow (5 nm), orange (7.5 nm), red (10 nm).

Effective contact angle at the block scale

An effective contact angle can be determined at the scale of a block of size L from the relationship, cos(𝜃 𝐿 ) = cos(𝜃 𝐶 ) 𝛼 𝐶 + cos(𝜃 𝑖𝑓 ) 𝛼 𝑖𝑓 + cos(𝜃 𝑖𝑐2.5 ) 𝛼 𝑖𝑐2.5 + cos(𝜃 𝑖𝑐5 ) 𝛼 𝑖𝑐5 + cos(𝜃 𝑖𝑐7.5 ) 𝛼 𝑖𝑐7.5 + cos(𝜃 𝑖𝑐10 ) 𝛼 𝑖𝑐10 (13) where 𝛼 𝐶 is the fraction of the interfacial surface between the solid phase and the pores in the block where the carbon phase is free of ionomer (thus 𝛼 𝐶 = 𝐴 𝐶 /𝐴 ). Similarly, 𝛼 𝑖𝑓 = 𝐴 𝑖𝑓 /𝐴, 𝛼 𝑖𝑐2.5 = 𝐴 𝑖𝑐-2.5 /𝐴, etc. The variations of the various surface fractions over cubic blocks of size L = 250 nn is shown in Fig. 10. Then, applications of Eq.( 13) leads to the results shown in Fig. 11. As can be seen, the averaging procedure represented by Eq.( 13) leads to significantly dampen the local fluctuations of the local effective contact angle. The average effective contact angle at the scale of the blocks of size 250 nn is 91.73°, thus close to the value 90.5° considered in [13] for the pores. The standard deviation of ± 0.5° is a bit less than the variability of ± 1.5° reported in [13]. The average value is 91.73° and the standard deviation is 0.5° over the 216 blocs considered.

However, the blocks of size L = 250 nn are larger than the average pore size of 60 nm reported in [31], whereas the effective contact angle considered in [13] is supposed to be defined at the pore scale. This question is investigated in Section 3.4, where pores are identified in the image.

Effective contact angle at the pore scale.

As discussed in Section 3.1 and exemplified in [13], the pore scale is of particular interest for characterizing the wettability properties of the CCL via the water retention curve (WRC) [13].

An effective contact angle is associated with each pore in [13] to determine the WRC from a representation of the pore space as a bundle of capillary tubes. As illustrated in Fig. 12, the pore space can be partitioned into pores using the watershed algorithm [28,66,67] available in PoreSpy [59]. For each pore so identified, a pore size can be defined as the radius of the sphere of equivalent volume. For the considered CCL image of 1500 nm ×1500 nm × 1500 nm, this leads to the pore size distribution depicted in Fig. 13. The average radius value over the 4698 pores identified in the image is 74 nm, which is about twice the average pore size reported in [31] using the maximum sphere inscription algorithm in PoreSpy [59]. As illustrated in Fig. 12b, the pores identified from the watershed algorithm are not spherical. Based on the pore shapes shown in Fig. 12, it is consistent that the pore volumes, and thus the pore equivalent radii, are greater than the volumes of the maximum spheres that can be inscribed in the pores. However, the shape of the histogram shows in Fig. 13 is similar to the one presented in [31] obtained with the maximum sphere inscription algorithm, but with the factor of about 2 in the pore size. Then, as for the blocks of size L, the various surface fractions 𝛼 𝐶 , 𝛼 𝑖𝑓 , 𝛼 𝑖𝑐2.5 , 𝛼 𝑖𝑐5 , 𝛼 𝑖𝑐7.5 , 𝛼 𝑖𝑐10 over the interfacial area between the solid phase and the pore space are computed for each pore in the 3D image. The effective contact angle 𝜃 𝑝𝑜 at the scale of the pores can then be determined as for the blocks according to the relationship, cos(𝜃 𝑝𝑜 ) = cos(𝜃 𝐶 ) 𝛼 𝐶 + cos(𝜃 𝑖𝑓 ) 𝛼 𝑖𝑓 + cos(𝜃 𝑖𝑐2.5 ) 𝛼 𝑖𝑐2.5 + cos(𝜃 𝑖𝑐5 ) 𝛼 𝑖𝑐5 + cos(𝜃 𝑖𝑐7.5 ) 𝛼 𝑖𝑐7.5 + cos(𝜃 𝑖𝑐10 ) 𝛼 𝑖𝑐10 (14) where the various contact angles are specified using the values reported in Table 2.

The variations of the various surface fractions 𝛼 𝐶 , 𝛼 𝑖𝑓 , 𝛼 𝑖𝑐2.5 , 𝛼 𝑖𝑐5 , 𝛼 𝑖𝑐7.5 , 𝛼 𝑖𝑐10 at the pore scale are shown in Fig. 14.

Compared to the variations of the surface fractions at the scale of blocks of size 250 nm (Fig. 10), it can be seen that the range of variations are greater in Fig. 5.13, which is a consequence of the consideration of many smaller volumes than the block size of 250 nm. On average, the hierarchy is similar with comparable values for 𝛼 𝑖𝑐2.5 and 𝛼 𝑖𝑐5 , and (𝛼 𝑖𝑐2.5 , 𝛼 𝑖𝑐5 ) > 𝛼 𝑖𝑐7.5 ,> 𝛼 𝑖𝑐10 , and 𝛼 𝐶 > 𝛼 𝑖𝑓 . The values of the effective contact angle at the pore scale so determined are plotted in Fig. 15.

The effective contact angle varies between about 85° and 95° (with even a few values outside this range) which is a significantly broader range than at the block scale of 250 nm (Fig. 11). This is consistent with the consideration of many volumes smaller than 250 x 250 x 250 nm 3 .

The average value is 91.7° and the standard deviation is 1.3°, which is quite comparable to the average value of 90.5° and the standard deviation of 1.5° considered in [13] for the secondary pores.

This is further illustrated in Fig. 15b, which shows the pore effective contact angle as a function of the pore radius. The variability of the contact angle is clearly greater for pores smaller than the average pore (of radius 74 nm) and clearly less for the larger pores. As can be seen from the curve obtained by running average (in red in Fig. 15b), the effective contact angle tends to increase with the pore size over the range of pore sizes below the average pore size (74 nm) and then stabilizes on average for the larger pores. From Fig. 14, it can be seen that the ionomer free carbon surface fraction tends be greater in the smaller pores as well as the surface fraction corresponding to the thickest coating ionomer film (10 nm). Both types of interface correspond to smaller values of the contact angle (78.5°

and 90° respectively from Table 2). As explained in [31], the thickest ionomer coating films are identified in concave zones of the solid phase interfacial surface. These zones are likely to be identified as small pores by the watershed algorithm. Thus, the fact that the contact angle is comparatively smaller in the smaller pores is consistent with our ionomer coating reconstruction procedure based on the local curvature of the interfacial surface between the pore space and the solid phase.

Water retention curve

A water retention curve can be determined analytically from the data shown in Fig. 15b as follows. Since the pores tend to be slightly hydrophobic, liquid water is considered as the nonwetting fluid and gas as the wetting fluid. Hence, the capillary pressure (at the scale of the considered CCL volume) is defined as 𝑃 𝑐 = 𝑃 𝑤 -𝑃 𝑔 , where 𝑃 𝑤 is the pressure in the liquid water and 𝑃 𝑔 is the pressure in the gas phase. Then the capillary pressure invasion threshold of a pore is defined as 𝑝 𝑐 (𝜃 𝑝𝑜 , 𝑟 𝑝𝑜 ) = -2𝛾cos (𝜃 𝑝𝑜 ) 𝑟 𝑝𝑜 (15) Thus 𝑝 𝑐 > 0 for a hydrophobic pore (𝜃 𝑝𝑜 > 90°) and 𝑝 𝑐 < 0 for a hydrophilic pore (𝜃 𝑝𝑜 < 90°). The pore radii 𝑟 𝑝𝑜𝑗 are ranked from smallest to largest. The index j varies from 1 to 4698, 4698 being the number 𝑁 𝑝𝑜 of pores identified in the considered CCL image. The volume of pore j is given by 𝑉 𝑝𝑜𝑗 =

The water retention curve, or capillary pressure curve, is obtained by plotting 𝑝 𝑐 (𝜃 𝑝𝑜𝑗 , 𝑟 𝑝𝑜𝑗 ) as a function of 𝑆 𝑗 . The obtained curve is shown in Fig. 16.

The curve is compared with the curve reported in [13] for the secondary pores. Although a purely statistical approach is developed in [13], the method used in [13] to obtain the curve is also based on the representation of the pore space as a bundle of capillaries. As can be seen, the range of pressure variation is greater in our case due to broader distributions of both the contact angle and the pore size. This points out the importance of the combined distributions of effective contact angle and pore size. It should be clear however that the approach used in [13] or in this paper to obtain the WRC is just indicative since it is based on an oversimplified representation of the pore space as a bundle of capillaries. Potentially important physical phenomena during the invasion of the pore space by water cannot be captured with this type of model. For instance, it is known [71] that the local invasion rules are different in the range of contact angles expected in the CCL (Fig. 15) than the classical invasion percolation rules often used to obtain the capillary pressure curve from pore network models, i.e. [72]. In particular, cooperative pore filling mechanisms between adjacent menisci [63,71,[START_REF] Tranter | Capillary Hysteresis in neutrally Wettable Fibrous Media: A Pore Network Study of a Fuel Cell Electrode[END_REF] are likely when the contact angle is close to 90°. In this respect, the CCL 3D images with the effective contact angle evaluated in each pore obtained in the present work are an interesting input for advanced computations of the capillary pressure curve based on the images.

Discussion

As explained in [31], our ionomer reconstruction method allows reconstructing the ionomer coating for a given value of the coverage and the ionomer filaments. As discussed in [31], a significant fraction of the ionomer is not identified and is considered to be mainly located between carbon particles inside the solid phase remaining after identification of the ionomer coating and ionomer filaments. However, this missing ionomer fraction can be expected to play a minor role as regards the wettability properties of the secondary pores because the missing ionomer is expected to be mostly localized away from the secondary pore spacesolid phase interfacial area considered in our approach. On the other hand, since the contact angle depends on the nature of the substrate for the sufficiently thin ionomer coating film and the surface fractions free of coating ionomer, a more refined approach where the missing ionomer would be identified is perhaps desirable.

The use of an effective contact angle also raises questions. For the effective contact angle concept to make sense, the scale of the contact angle heterogeneities must be sufficiently separated from the scale at which the effective contact is computed. Based on the size of the smaller pores identified with the watershed algorithm (as low as 10 -20 nm in Fig. 15), and the relatively sparse distribution of the platinum particles in Fig. 2a (over a domain of size 250 nm), one can wonder whether the effective contact distribution at the pore scale is well captured by our approach due to the lack of length scale separation. Consider the distribution of local contact angle obtained with the fully deterministic approach illustrated in Fig. 6. An effective contact angle can be computed at the scale of the image shown in Fig. 6

where 𝑁 𝑖𝑣 is the number of interfacial voxels and 𝜃 𝑗 is the contact angle assigned to voxel j according to the values given in Table 1.

For the block of size 250 nm illustrated in Fig. 6, a value of 91.6° is obtained from Eq. ( 17) whereas a value of 91.26° is obtained for this particular block with the hybrid approach presented in Section 3.3. Although the comparison is only for a particular block, this suggests that the hybrid approach gives results comparable to the direct approach in terms of effective contact angle at this scale.

The results also depend on the values of the intrinsic contact angles, i.e. the contact angle on the carbon free of ionomer and platinum particles or the contact angle on the platinum free of ionomer. As indicated before we used a value of 80° for carbon free of ionomers and platinum, which is an average value based on the values proposed in [12]. As discussed in [12], the value can be lower due to the presence of edges, defects and amorphous sections involved and can decrease to 60°-70°. A value of 62° is for instance used in [13]. Although the surface fraction corresponding to carbon free of ionomer is not very large (0.2 on average), considering a lower value for the intrinsic contact angle on carbon will increase a bit the range of the effective contact angle values at the pore scale. Regarding the platinum, a value of 10° is considered in [13]. This is slightly greater than the value 0 we have used but the impact of this variation should be minor since the values of the corresponding cosines are very close and the surface fractions corresponding to platinum quite low.

A more serious shortcoming than the uncertainty on the intrinsic contact angle values lies in the fact that only a dry ionomer is considered in our work. As pointed out in [13], the interaction of the ionomer with the various substrates can be more complex than considered in the present work (and in [13]) under humid conditions. Future works are therefore necessary to consider the impact of the humidity conditions.

Conclusion

The wettability properties of the porous network in fuel cell catalyst layers are crucial for the water management of PEMFC and a key input data for the models. This work provides an approach to determine the contact angle distribution in three dimensional digital images of the CCL microstructure. The approach was illustrated from 3D images of a CCL obtained by FIB-SEM after reconstruction of the ionomer coating and ionomer filaments, but it can be applied as well on synthetic images generated numerically. Our results are an illustration of the key role of the ionomer on the wettability properties. The latter highly depend on the ionomer distribution within the CCL microstructure and in particular on the ionomer coating the carbon particle agglomerates. The impact of ionomer on the wettability properties was modelled using the approach presented in [13] considering the molecular alignment of ionomer sidechains and backbones to be decisive for the resulting wetting properties and taking into account the inversion phenomenon of the wetting properties of the Pt/C support for sufficiently thin ionomer coating films.

Various scales of interest were considered. The reconstruction of the contact angle distribution at the smallest scale, i.e. a sufficiently small voxel size of 0.5 nm for considering the platinum particles individually in the image, provides 3D digital images that can be used as input data for the detailed simulations of two-phase flow via direct numerical simulation (DNS) techniques. The consideration of the larger scale corresponding to voxelized images obtained after the ionomer coating reconstruction leads to the definition of an effective contact angle for each interfacial voxel. This also leads to 3D images well suited for twophase flow direct numerical simulations and the image-based computations of important properties for macrohomogeneous models such as the water retention curve or the relative permeabilities. The concept of effective contact angle was also considered at the pore scale after segmentation of the pore space into pores. This leads to associate an effective contact with each pore. The 3D images so obtained are well suited for simulations of two-phase flows and computations of the water retention curve and relative permeabilities from PNM simulations or image-based techniques such as the full morphology approach. The effective contact angle was also characterized at the scale of CCL sub-blocks of size 250 x 250 x 250 nm 3 . The results indicate that the effective wettability does not vary significantly at this scale within the CCL.

Our results allow to investigate the correlation between the pore size and effective contact angle. The trend is a lower effective contact angle on average in the smaller pores with a greater variability in the smaller pores than in the larger pores.

As illustrated from a simple model, the combined information on the pore size and the effective contact angle for each pore can be used to obtain the water retention curve. This, however, cannot replace the characterization of the water retention curve from DNS or PNM simulations as evoked above due to the oversimplification of the pore space geometry as a bundle of tubes in the simple model. This is especially important in the case of the CCL because the effective contact angle at the pore scale in the CCL is close to 90°. It is known that the local rules controlling the pore filling mechanisms and thus the water distribution in a porous medium are particularly sensitive to the contact angle distribution when the latter is close to 90° [71].

In summary, this work contributes to bridging the detailed representation of the microstructure and its wettability properties to macroscopic catalyst layer properties. catalyst surface. Kongkanand et al. have highlighted that this resistance can represent a significant fraction of the overall performance limitation, equivalent to the global oxygen pore diffusion resistance considered here.

Another significant performance limitation may arise from proton conduction or water management, which we explored in Chapters 3 to 5. However, quantifying these transport phenomena requires a comprehensive understanding of the catalyst structure, including ionomer distribution and platinum content. Such details cannot be provided by FIB-SEM images alone. This led us to develop a reconstruction model in Chapter 3 to obtain local ionomer coating and platinum distribution. The agreement between our findings and existing literature data lends credibility to our hypothesis that ionomer primarily coats the cavities of carbon agglomerates. The advantage of this ionomer coverage model is that it combines the catalyst layer microstructure from FIB-SEM with a modeling approach, minimizing biases that may be present in the fully randomly reconstructed synthetic structures often used in the literature.

To validate and further investigate our ionomer reconstruction approach, we computed the effective transport properties of the catalyst layer and compared them to experimental data. In Chapter 4, we evaluated the proton conduction within the ionomer film through direct numerical simulations. Our fitted tortuosity law exhibited very good agreement with experimental results from the literature, addressing the underestimation typically observed in numerical simulations based on synthetic 3D images compared to experiments. We also demonstrated that this gap can be corrected by considering proton transport in the carbon particle microporosity, which is generally considered to be liquid water-filled under PEMFC usual operation conditions.

Water management can significantly affect the effective transport properties, either for oxygen or protons. The presence of liquid water in the secondary pores can introduce additional oxygen diffusion resistance, while liquid water promotes proton transport. Chapter 5 quantified the catalyst layer wetting properties at different scales, indicating that the effective contact angle at the pore scale vary over a narrow range around 90° within the CCL, which is consistent with recent studies in the literature. We also evaluated the capillary pressure curve using a simple model.

In term of future work, the fully reconstructed CCL sample can be used for direct numerical simulations, considering gas and proton transport coupled with electrochemical reactions on the platinum catalyst's reconstructed surface. This is planned as part of the Further-FC project by partner DLR. Additionally, we suggest continuing the investigation of the catalyst structure by considering micropores within the carbon phase. An option could be to generate numerically micropores into the previously reconstructed carbon, accounting for platinum within these micropores, though challenges regarding sample size remain.

As illustrated in this work, in spite of the encouraging results obtained with our ionomer reconstruction approach in the FIB-SEM images, the reconstruction of the morphology of the ionomer and its distribution in the CCL microstructure remains a challenge. Within the framework of our approach, the "missing" ionomer for instance is an issue. The fact that our approach leads to better results on the proton conductivity than previous works based on synthetic images would deserve to be studied in more details so as to characterize the morphological differences between the synthetic images from the literature and our images.

Overall, predicting the proton conductivity remains a challenge.

The question of the impact of liquid water in the secondary pores in the CCL has been only very superficially addressed in the present work with the computation of the capillary pressure curve (CPC) using an oversimplified analytical approach. Simulations of the liquid water within the secondary pores and the CPC should be undertaken by simulations on our reconstructed images. This can be presumably done using various techniques (PNM or full morphology for instance). However, this is not as straightforward as one could consider based on the works on this topic in the literature. The fact that the contact angle is close to 90° is a difficulty due to the occurrence of local cooperation mechanisms between adjacent menisci controlling the pore body invasion. This difficulty is overlooked in most works in the PEMFC literature. These mechanisms need to be considered in the simulation algorithms.

A somewhat natural outcome for the results presented in the present thesis lies in the macrohomogeneous (continuum) models coupling the various layers (GDM-MPL-CCL) considered in this study, and taking into account gas diffusion, proton conduction, liquid water formation and transport and ORR kinetics. This model should incorporate the effective transport properties computed in this study. This type of model facilitates the quantification of different performance limitations and is better suited than the direct simulations for establishing recommendations in order to optimize the design of PEMFCs. This work is

Contents 1 . 1

 11 General context .......................................................................................................... 1.2 Applications: .............................................................................................................. 1.3 PEM fuel cell structure: ............................................................................................. 1.4 Cell performance limitations ...................................................................................... 1.5 PEM fuel cell modeling: ............................................................................................ 1.6 Further-FC project ...................................................................................................... 1.7 Outline of PhD thesis manuscript .............................................................................. 1.8 References ..................................................................................................................

Fig. 1 .

 1 Fig.1. Atsushi Tsutsumi. Ragone plot comparing the performance of various electrochemistry systems (energy density vs. output density) [3].

Fig. 3 .

 3 Fig.3. a) ZeroAvia Successfully Completes Initial Dornier 228 Flight Test Campaign [6] . b) Energy observer 2 | KADEG BOUCHER/EPRON DESIGN

Fig. 6 .

 6 Fig.6. Detail structure of GDL (SGL 25BC), MPL, and CL[9]

Fig. 8 .

 8 Fig.8. Analysis of performance limitations based on recent study literature presented in [11].

Fig. 9 .

 9 Fig.9. Comparison between different 3D tomography techniques for volume of analyzed material and voxel resolution [24].

Fig. 1 .

 1 Fig.1. (a) 3D segmented FIB-SEM image of the MPL (pore space in light grey, solid phase in black), (b) MPL porosity as a function of averaging volume (see text), (c) MPL pore size distribution.

Fig. 2 .

 2 Fig.2. a) Cross-section image of the GDM-MPL assembly extracted from XCT analysis, b) segmented 2D image of the GDM -MPL assembly with the MPL in red and the GDM in green, c) 3D rendering of the GDM-MPL assembly, d) 3D rendering of the GDM, e) side view showing traversing cracks, f) GDM-MPL assembly macro-porosity variation in the through-plane direction (i.e. along the z direction shown in Fig. 2e). Insets: crack width distribution, distribution of in-plane distance between cracks, image of MPL cracks, GDM pore size distribution and image of the GDM fibrous medium.

Fig. 3 .

 3 Fig. 3. (a) Local pore size radius in a MPL sub-sample from MSI method (solid phase in light grey, inscribed spheres in pore space in color), (b) Knudsen number distribution in MPL

Fig. 4 .

 4 Fig.4. a) Sketch of direct numerical simulations with inlet-outlet concentration difference applied in the x-direction. b) fictitious continuum medium equivalent block, (c) and (d)

Fig. 5 .

 5 Fig.5. Distribution of the components of MPL effective diffusion tensor, (a) with block size L 1 = 500 nm and without Knudsen diffusion, (b) with block size L 1 = 500 nm and with Knudsen diffusion, (c) with block size L 1 =1500nm and without Knudsen diffusion, (d) with block size L 1 =1500nm and with Knudsen diffusion. (e) Probability density function of the relative diffusivity tensor component 𝐷 𝑒𝑓𝑓𝑥𝑥-𝐿 1 , for both block sizes, with and without considering Knudsen diffusion

Fig. 6 .

 6 Fig.6. Comparison of the evolution of 𝜀𝐷 𝑒𝑓𝑓𝑥𝑥-𝐿 1 𝐷 𝑏𝑢𝑙𝑘 , i.e. 𝜀 𝜏 , as a function of porosity with

  ℓ 1 is the vertical distance between first fiber (bottom) and average position of MPL bottom surface. The thickness of the GDM-MPL overlap region is denoted by ℓ 2 : ℓ 2 is the vertical distance between the most advanced fiber and the average position of MPL bottom surface. The thickness of the MPL region free of fibers is ℓ 3 : ℓ 3 is the vertical distance between the most advanced fiber and the average position of MPL top surface.A horizontal reference plane is defined at some distance from the MPL top surface to facilitate the implementation of the top boundary condition. The thickness of this additional layer is ℓ 4 : ℓ 4 is the average vertical distance between the MPL top surface average position and this plane. To compute the assembly through-plane diffusion coefficient, the top BC (C = c 2 ) is applied over this top plane surface. The region in blue is a region where a very large value of the diffusion coefficient (1 m 2 /s, thus 10 5 greater than the diffusion coefficient in the pores) is imposed so that the diffusion transport resistance of this region is negligible. For this reason, the value of the thickness of the additional layer, i.e. ℓ 4 , has no impact on the result. The bottom boundary condition (C = c 1 ) is imposed over the plane adjacent to the first fiber voxel in the vertical direction in the image.

Fig. 7 .

 7 Fig. 7. a) Illustration of multiscale approach combining FIB-SEM (MPL matrix) and X-ray tomography images (GDM-MPL assembly). The MPL 3D cracks network and the penetration of MPL (in red) in GDM (in green) are visible in the X-ray tomography images. b) Sketch of the various layers and equivalent resistance model. Expressions of the various resistances are given in Table3.

  Determining the ionomer three-dimensional distribution in the microstructure of the cathode catalyst layer (CCL) of a polymer electrolyte membrane fuel cell (PEMFC) is pivotal in order to perform numerical simulations of the electro-chemical and transport processes directly at the scale of the microstructure and develop optimized design of CCL. In this work, the ionomer distribution is reconstructed from a 3D image of the CCL obtained by focused ion beam-scanning electron microscopy (FIB-SEM). The SEM image contrast of the solid phase does not allow to separate the carbon, the ionomer and the catalyst. In order to identify and distribute the ionomer in the segmented solid phase, numerical techniques are used, based on physical considerations. Two forms of ionomer are distinguished: the ionomer filaments and the coating ionomer. The ionomer filaments are identified by applying a filtering procedure to the size distribution in the segmented solid phase obtained via a maximum sphere inscription algorithm. The coating ionomer is identified via a filtering procedure of the curvature distribution computed over the segmented solid surface based on the consideration that concave regions 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 agreement with the available experimental data. In a last step, catalyst particles, also not visible in the FIB-SEM image, are distributed over the carbon phase surface.Keywords: Ionomer, Pt/C, Coverage, I/C, CCL, PEMFC

Fig. 1 .

 1 Fig.1. a) Schematic of segmentation procedure, b) CCL segmented microstructure from FIB-SEM image with segmented solid phase in light gray and pore space in black.

Fig. 2 .

 2 Fig.2. HRTEM images of different CCL areas where the ionomer is clearly visible together with carbon phase; (a) example of ionomer filament, (b) coating ionomer. In (b) the detectable limit between the ionomer and the carbon phase surface is indicated by a yellow line.

Porespy [ 54 ]Fig. 3 (

 543 Fig.3 (a) Pore Size Distribution (PSD) obtained from using the maximum inscribed sphere algorithm in the void phase, (b) Size Distribution (SSD) obtained from using the maximum inscribed sphere algorithm in the solid phase. (c) Local porosity distribution over the full segmented image as a function of the averaging volume size.

Fig. 4 .

 4 Fig.4. (a) Cross section of membrane with the CCL (on top) and anode catalyst layer; (b) Computed CCL thickness as a function of macro-porosity for various values of the carbon particle micro-porosity; (c) Computed solid phase fractions as a function of carbon particle micro-porosity.

Fig. 5 .

 5 Fig.5. (a) 3D visualization of the Maximum Inscribed Sphere size distribution in a sample of 250x250x250 nm 3 . (b)-(c) detected filament distribution in image (a) for a upper thresholding size of 18 nm (b) and 25 nm (c) corresponding to relative ionomer volume fraction of 7%, (b) and 25 % (c) respectively. (d-f) Identification of ionomer filament (in green) in a sample of 250x750x750 nm 3 : (d) upper thresholding size of 18 nm, (e) upper thresholding size of 25 nm, (f) upper thresholding size of 30 nm.

(Fig. 6 .

 6 Fig.6. (a) Illustration of curvature computation procedure on the voxel structure after exclusion of ionomer filaments; (b) 3D rendering of curvature distribution at the surface of the solid phase in a sub-region of the CCL of 250 x 250 x 250 nm 3 . The curvature in a concave zone is negative due to the chosen unit normal vector orientation.

Fig. 7 .

 7 Fig.7. Examples of ionomer coverage determined from thresholding the local solid surface curvature distribution in a sub-region of the CCL of 250 x 250 x 250 nm 3 . The coverage increases with an increasing curvature threshold.

Fig. 8 .

 8 Fig.8. (a) Ionomer coverage as a function of curvature threshold (see text); (b) Computed ionomer volume fraction (with and without consideration of the ionomer filaments) as a function of ionomer coverage for two values of coating ionomer film thickness (assumed spatially uniform) with comparison with some studies from the literature.

Fig. 9 (

 9 Fig.9 (a) Coating ionomer film thickness histogram based on the curvature obtained considering a 1500 x 1500 x 1500 nm 3 image, (b) 3D rendering of coating ionomer film thickness distribution obtained for the coverage of 82% determined from the considered correlation (Fig.10a) between the local surface curvature and the film thickness in a sub-

Fig. 10 .

 10 Fig.10. Computed coating ionomer distribution in a CCL slice of 75 nm in width and 1500 nm in lateral extension for a coverage of 0.82. In the left column, only the coating (in blue)

Fig. 11 .

 11 Fig.11. Illustration of the coating ionomer percolating (in blue) and non-percolating (in red) coating ionomer clusters in a CCL sample of 1500 x 1500 x 1500 nm 3 for four different values of the coverage. The right column shows only the percolating cluster.

Fig. 12 .

 12 Fig.12. Reconstructed CCL sample of 250 x 250 x 250 nm 3 with platinum particles after sub-

Fig. 1 .

 1 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 3 , y is the through-plane direction, z is the FIB-SEM milling direction.

Fig. 2 .

 2 Fig.2. a) 3D rendering of the segmented solid phase. b) reconstructed 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 3 .

Fig. 3 .

 3 Fig.3. Distribution of the dimensionless components of CCL oxygen effective tensor 𝐃 𝑂 2 𝑒𝑓𝑓 : (a) for block size L 1 =1500 nm without Knudsen diffusion, (b) for block size L 1 =1500 nm with Knudsen diffusion, (c) same as in (a) but for components of tensor 𝐃 𝑂 2𝑒𝑓𝑓 *

Fig. 4 .

 4 Fig.4. Illustration of the through plane oxygen transport conductivity in the various blocks of size L 1 .

Fig. 5 .

 5 Fig.5. Fitting of the evolution of the dimensionless oxygen effective diffusion tensor 𝐃 𝑂 2 𝑒𝑓𝑓 diagonal components as a function of porosity with the analytical model of Bruggeman.

Fig. 6 .

 6 Fig.6. Comparison of through-plane effective diffusivity obtained in the present work with other pore -scale simulation results and experimental data. 𝐷 𝑏𝑢𝑙𝑘 in this figure is a shorthand for 𝐷 𝑂 2 -𝑏𝑢𝑙𝑘 .

Fig. 7 .

 7 Fig.7. Distribution of the dimensionless components of CCL proton effective conductivity tensor 𝛔 𝑃𝑒𝑓𝑓 *

Fig. 8 .

 8 Fig.8. Distribution of the diagonal dimensionless components of CCL proton effective conductivity tensor 𝛔 𝑃𝑒𝑓𝑓 *

Fig.8c shows

  Fig.8cshows that the proton conductivity decreases with the increasing block porosity, consistently with Fig.7a, 8a and 8b since the ionomer volume fraction, the coverage and the I/C ratio increases with a decreasing block porosity. The spatial heterogeneity of the through plane diagonal component (σ 𝑃𝑒𝑓𝑓,𝑦𝑦 *

Fig. 9 .

 9 Fig.9. Comparison of the evolution of the dimensionless through plane effective proton

Fig. 10 ,

 10 Fig.10,

Fig. 10 .

 10 Fig.10. Local proton conductivity distribution in the representative block (𝜀 = 0.52, 𝜀 𝐼 = 0.126, coverage = 0.82). Coating ionomer and filament in red. Pore space in darker blue (within the framework of the used penalization technique a very low conductivity value is assigned to pore voxels), complementary solid phase in lighter blue.

3Fig. 11 .

 11 Fig.11. Variation of dimensionless through-plane proton effective conductivity 𝜎 𝑒𝑓𝑓𝑦𝑦-𝐿 2 *

Fig. 1 .

 1 Fig.1. a) 3D rendering of solid phase (in green). b) Example of coating ionomer distribution with coating ionomer in blue and ionomer filaments in green. The platinum loaded carbon phase free of ionomer is in grey. c) Example of coating ionomer film thickness distribution. Ionomer filament in green. Platinum loaded carbon phase free of ionomer in grey. All images are for a sub-region of the CCL of 250 x 250 x 250 nm 3 with a voxel size of 2.5 x 2.5 x 2.5 nm 3 .

Fig. 2 .

 2 Fig.2. a) CCL sample of 250 x 250 x 250 nm 3 with reconstructed platinum particle distribution after subdivision of the sample in 0.5 x 0.5 x 0.5 nm 3 voxels (carbon phase in grey, platinum particles in red, the reconstructed ionomer is removed in this image), b) same sample as in (a) with the reconstructed ionomer (ionomer coating in blue, ionomer filaments

Fig. 3 .

 3 Fig.3. Normalized platinum particle size histogram from [30].

Fig. 4 .

 4 Fig.4. Plot of degree of alignment as a function of ionomer film thickness for the platinum substrate.

Fig. 5 .

 5 Fig. 5. a) Reconstructed ionomer distribution in a block of size of 250 x 250 x 250 nm 3 with voxel size of 0.5 x 0.5 x 0.5 nm 3 , b) same image as in a) with indication of closer substrate below the ionomer coating interfacial voxels (ionomer coating on platinum in orange, ionomer coating on carbon in blue), c) same image as in a) (but slightly tilted) showing the local contact angle distribution (expressed in degree in the color bar).

Fig. 6 .

 6 Fig. 6. Same image as in Fig.5c but enlarged for better visibility: local contact angle distribution obtained for the coverage of 82% in a sub-region of the CCL of 250 x 250 x 250 nm 3 with voxel size of 0.5 x 0.5 x 0.5 nm 3 . The color bar gives the contact angle in degree.

  2 mg/cm 2 for the considered CCL). Then consider a CCL cubic image of volume 𝐿 × 𝐿 × 𝐿. The mass 𝑀 𝑝𝑡 of platinum in this image is on average equal to 𝑀 𝑝𝑡 = (𝑚 𝑝𝑡 × 𝐿 × 𝐿) × ( 𝐿 𝑙 𝐶𝐶𝐿 ) where 𝑙 𝐶𝐶𝐿 is the CCL thickness (𝑙 𝐶𝐶𝐿 = 6 µm for the considered CCL). The mass of platinum on the carbon phase external surface in the considered image of size L is therefore 𝑀 𝑝𝑡-𝑆 = 0.25𝑀 𝑝𝑡 = 0.25(𝑚 𝑝𝑡 × 𝐿 × 𝐿) × ( Consider one of the surfaces in Eq.(3) (except the ionomer filament surface). Let denote such a surface by 𝐴 𝑗 . The mass of platinum over this surface is 𝑀 𝑝𝑡-𝑆-𝑗 = 𝑀 𝑝𝑡-

  𝜂 𝑣𝑜𝑥𝑒𝑙 = 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑣𝑜𝑥𝑒𝑙𝑠 𝑖𝑛 𝑡ℎ𝑒 𝑏𝑙𝑜𝑐𝑘 𝑓𝑜𝑟𝑚𝑖𝑛𝑔 𝑡ℎ𝑒 𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 𝑛𝑢𝑚𝑏𝑒𝑟 𝑜𝑓 𝑣𝑜𝑥𝑒𝑙𝑠 𝑖𝑛 𝑡ℎ𝑒 𝑏𝑙𝑜𝑐𝑘 . To obtain 𝜂 from 𝜂 𝑣𝑜𝑥𝑒𝑙 , we apply the relationship 𝜂 ≈ 𝜂 𝑣𝑜𝑥𝑒𝑙 𝐿 𝑣𝑜𝑥𝑒𝑙

Fig. 7 .

 7 Fig.7. Variation of voxel specific surface area 𝜂 𝑣𝑜𝑥𝑒𝑙 over 216 blocks of size 250 x 250 x 250 nm 3 paving a CCL image of 1500 x 1500 x 1500 nm 3 .

Fig. 9 .

 9 Fig.9. 3D rendering of local effective contact angle distribution obtained for the coverage of 82% in a sub-region of the CCL of 250 x 250 x 250 nm 3 with voxel size 2.5 x 2.5 x 2.5 nm 3 .

Fig. 10 .

 10 Fig.10. Surface fractions of the various components over the interfacial surface area between the solid phase and the pores in blocks of size 250 x 250 x 250 nm 3 with a voxel size of 2.5 × 2.5 × 2.5 nm 3 : a) carbon phase free of ionomer and ionomer filaments, b) ionomer coating for the four thicknesses considered.

Fig. 11 .

 11 Fig.11. Variation of effective contact angle at the scale of blocks of size 250 x 250 x 250 nm 3 .

Fig. 12 .

 12 Fig.12. a) 1500 nm × 1500 nm 2D image of the various reconstructed components in the CCL: pore space in dark blue, carbon phase in light blue, ionomer filaments in yellow, ionomer coating in green. b) partition of the pore space in pores using the watershed algorithm in PoreSpy [59]. The dark blue color corresponds to the solid phase (carbon and ionomer) whereas each colored region in the pore space corresponds to a pore. c) Effective contact angle distribution in the pores.

Fig. 13 .

 13 Fig.13. Pore size histogram obtained after application of the watershed algorithm [59].

Fig. 14 .

 14 Fig.14. Surface fractions of the various components over the interfacial surface area between the solid phase and the pore space in the pores identified in a 3D CCL image of size 1500 x 1500 x 1500 nm3 with a voxel size of 2.5 x 2.5 x 2.5 nm3. Carbon phase free of ionomer and ionomer filaments in the main plot. Ionomer coating for the four thicknesses considered in the inset.

Fig. 15 .

 15 Fig.15. a) Variation of effective contact angle at the pore scale; b) Variation of effective contact angle in pores as the function of pore size. The red curve corresponds to a running average over 100 successive values.

  of the pore equivalent radius. The liquid water saturation corresponding to a given capillary pressure threshold 𝑝 𝑐 (𝜃 𝑝𝑜𝑗 , 𝑟 𝑝𝑜𝑗 ) is then computed as,

Fig. 16 .

 16 Fig.16. Water retention curve.

  from the local contact angle distribution shown in Fig. 6 using the relationship, cos(𝜃 𝐿-𝑑𝑖𝑟𝑒𝑐𝑡 )

Table 1 :

 1 MPL effective diffusion tensor relative components (

	1 𝜏 ) for both block sizes

Table 2 .

 2 GDM and GDM-MPL assembly through plane diffusion properties.

				GDM without	GDL	with	GDL	with
				MPL		cracked MPL	crack free MPL
	Thickness ℓ (µm)				220	236	236
	Porosity				0.88	0.81	0.8
	Through plane effective relative	0.8	0.66	0.6
	diffusion coefficient	1 𝜏
	Through plane diffusive resistance	1	1.4	1.6
	𝑅/𝑅 𝐺𝐷𝑀-𝑤𝑜𝑀𝑃𝐿 (𝑅 =	ℓ𝜏 ℰ𝐷 𝑏𝑢𝑙𝑘	)

  Contents 3.1 Introduction ................................................................................................................ SEM........................................................................................................... 3.2.2 Transmission Electron Microscopy (TEM) ...................................................... 3.3 Size distributions in void and solid phases ................................................................ 3.4 Impact of focused ion beam milling on ionomer .......................................................

	3.2.1	FIB-

3.2 Characterization of the CCL by FIB-SEM, Transmission Electron Microscopy (TEM) ...................................................................................................................................

3.5 Ionomer filaments and coating ionomer film ............................................................. 3.5.1 Ionomer filaments ............................................................................................. 3.5.2 Coating ionomer film. Computation of the solid phase surface local curvature .. .......................................................................................................................... 3.5.3 Coating ionomer film. Coverage modeling ......................................................

3.5.4 Coating ionomer film. Ionomer film thickness distribution ............................ 3.6 Platinum particles ....................................................................................................... 3.7 Discussion .................................................................................................................. 3.8 Conclusion.................................................................................................................. 3.9 References ..................................................................................................................

Table 1 .

 1 Filament relative ionomer fraction (fraction of total ionomer volume detected as

	Upper thresholding size	18 nm	25 nm	30 nm
	Filament relative	7%	25 %	48%
	Ionomer Fraction
	Relative Solid Fraction	2.3%	8.1%	16%

ionomer filaments) and filament ionomer relative solid fraction (fraction of total volume of solid phase) obtained for various values of the upper thresholding size (see text).

Table 2 .

 2 Filament ionomer, coating ionomer and not detected ionomer volume fractions relative to the total volume of ionomer in the CCL for various values of the coverage for the I/C ratio of 0.8, an average ionomer coating film thickness of 5 nm and the filament upper filtering size of 18 nm.

	(with ionomer	0.28	0.47 0.67	0.84	1.0
	filaments)
	Filament ionomer relative	7%	7 %	7%	7%	7%
	volume Fraction
	Coating ionomer relative	16% 28% 43% 54% 66%
	volume Fraction
	Not detected ionomer	77% 65% 50% 39% 27%
	relative volume fraction

Table 1 .

 1 Computational domain sizes and cubic finite volume sizes considered in the

	simulations
	L 1	L 2	Computational
			finite volume size

Table 2 :

 2 Dimensionless components of CCL oxygen effective diffusion tensor 𝐃 𝑂 2𝑒𝑓𝑓 , i.e. at the CCL scale (L 2 ). The through plane coefficient is 𝐷 𝑒𝑓𝑓𝑦𝑦-𝐿 2 . 𝐷 𝑏𝑢𝑙𝑘 is a shorthand

	( 1 𝜏 ), for 𝐷 𝑂 2 -𝑏𝑢𝑙𝑘
		𝐷 𝑒𝑓𝑓𝑥𝑥-𝐿 2	𝐷 𝑒𝑓𝑓𝑥𝑦-𝐿 2	𝐷 𝑒𝑓𝑓𝑥𝑧-𝐿 2	𝐷 𝑒𝑓𝑓𝑦𝑥-𝐿 2	𝐷 𝑒𝑓𝑓𝑦𝑦-𝐿 2	𝐷 𝑒𝑓𝑓𝑦𝑧-𝐿 2	𝐷 𝑒𝑓𝑓𝑧𝑥-𝐿 2	𝐷 𝑒𝑓𝑓𝑧𝑦-𝐿 2	𝐷 𝑒𝑓𝑓𝑧𝑧-𝐿 2
		𝐷 𝑏𝑢𝑙𝑘	𝐷 𝑏𝑢𝑙𝑘	𝐷 𝑏𝑢𝑙𝑘	𝐷 𝑏𝑢𝑙𝑘	𝐷 𝑏𝑢𝑙𝑘	𝐷 𝑏𝑢𝑙𝑘	𝐷 𝑏𝑢𝑙𝑘	𝐷 𝑏𝑢𝑙𝑘	𝐷 𝑏𝑢𝑙𝑘
	With	0.14	0.0001	0.002	-0.001	0.17	0.05	0.0001	0.048	0.24
	Knudsen
	diffusion
	Without	0.363	0.004	0.003	0.0026	0.457	0.11	0.0017	0.10	0.578
	Knudsen
	diffusion

Table 3 :

 3 Dimensionless components of CCL oxygen effective diffusion tensor 𝐃 𝑂 2𝑒𝑓𝑓 * at the CCL scale (L 2 ). The through plane coefficient is 𝐷 𝑒𝑓𝑓𝑦𝑦-𝐿 2 . 𝐷 𝑏𝑢𝑙𝑘 is a shorthand for 𝐷 𝑂 2 -𝑏𝑢𝑙𝑘 .

	, i.e.

  at scale L 2 . The through plane coefficient is 𝜎 𝑃𝑒𝑓𝑓𝑦𝑦-𝐿 2 .

	𝜎 𝑒𝑓𝑓𝑥𝑥-𝐿 2	𝜎 𝑒𝑓𝑓𝑥𝑦-𝐿 2	𝜎 𝑒𝑓𝑓𝑥𝑧-𝐿 2	𝜎 𝑒𝑓𝑓𝑦𝑥-𝐿 2	𝜎 𝑒𝑓𝑓𝑦𝑦-𝐿 2	𝜎 𝑒𝑓𝑓𝑦𝑧-𝐿 2	𝜎 𝑒𝑓𝑓𝑧𝑥-𝐿 2	𝜎 𝑒𝑓𝑓𝑧𝑦-𝐿 2	𝜎 𝑒𝑓𝑓𝑧𝑧-𝐿 2
	𝜎 𝑝	𝜎 𝑝	𝜎 𝑝	𝜎 𝑝	𝜎 𝑝	𝜎 𝑝	𝜎 𝑝	𝜎 𝑝	𝜎 𝑝
	0.14	0.001	-0.004	-0.002	0.2	0.112	-0.004	0.116	0.23

𝜏

)

Table 5 :

 5 Dimensionless components of CCL proton effective conductivity tensor 𝛔 𝑃𝑒𝑓𝑓 *

	, i.e.

Table 1

 1 Contact angle over the pore / solid phase interfacial voxels

		Material	Contact angle (°)
	𝜃 𝐶	Carbon (free of ionomer and	80°
		platinum)
	𝜃 𝑃𝑡	Platinum (free of ionomer)	0°

  allows to express N as,

	𝑎 𝑝𝑡 =	𝐴 𝑗-𝑝𝑡 𝐴 𝑗	=	𝑁 ∫ -∞ +∞	𝑓(𝑑 𝑝 ) 𝐴 𝑗	𝜋 4 𝑑 𝑝 2 𝑑(𝑑 𝑝 )	(7)
	Substituting the expression of 𝑁 into Eq.(7) leads to
			𝑎 𝑝𝑡 =	𝐴 𝑗-𝑝𝑡 𝐴 𝑗	=	3 2	(	0.25(𝑚 𝑝𝑡 ×𝐿×𝐿)×( 𝑙 𝐶𝐶𝐿 𝐿 𝐴 𝜌 𝑝𝑡	) )	∫ -∞ +∞ ∫ +∞ -∞	𝑓(𝑑 𝑝 ) 𝑓(𝑑 𝑝 ) 𝑑 𝑝 𝑑 𝑝 2 𝑑(𝑑 𝑝 ) 3 𝑑(𝑑 𝑝 )
	𝑁 =	0.25(𝑚 𝑝𝑡 ×𝐿×𝐿)×( 𝜌 𝑝𝑡 ∫ 𝑓(𝑑 𝑝 ) +∞ -∞ 6 𝑑 𝑝 3 𝑑(𝑑 𝑝 ) 𝐿 𝜋 𝑙 𝐶𝐶𝐿 )	𝐴 𝑗

𝐴

(6)

The next step is to determine the surface fraction of platinum in 𝐴 𝑗 . The latter is given by

Table 2 Table 2 .

 22 Local effective contact angle over the pore / solid phase interfacial area

		Material	Contact angle (°)
	𝜃 𝐶	Carbon (free of ionomer)	78.5
	𝜃 𝑖𝑐2.5 Ionomer coating (2.5 nm)	100°
	𝜃 𝑖𝑐5 Ionomer coating (5 nm)	93°
	𝜃 𝑖𝑐7.5 Ionomer coating (7.5 nm)	91°
	𝜃 𝑖𝑐10 Ionomer coating (10 nm)	90°
		𝜃 𝑖𝑓 Ionomer filament	90°

Remerciements

Cette Joint Undertaking reçoit le soutien du programme de recherche 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 and Innovation program, Hydrogen Europe and Hydrogen Europe Research.

Abbreviations

This missing volume fraction is denoted by 𝜀 𝐼-𝐶𝑃 . The complementary segmented solid phase (subscript CP) 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 × 250 × 250 nm 3 ). Then 𝑉 = 𝑉 𝑣𝑜𝑖𝑑 + 𝑉 𝑃𝑡 + 𝑉 𝐶 + 𝑉 𝐼 (A-1)

Where 𝑉 𝑣𝑜𝑖𝑑 is the volume of pores, 𝑉 𝑃𝑡 is the volume of platinum, 𝑉 𝐶 is the volume of carbon and 𝑉 𝐼 is the volume of ionomer.

We also have

𝑉 𝑣𝑜𝑖𝑑-𝑚𝑎𝑐𝑟𝑜 is the volume of macropores (the one identified in the FIB-SEM image, also referred to as the secondary pores) and 𝑉 𝑣𝑜𝑖𝑑-𝑚𝑖𝑐𝑟𝑜 is the volume of micro-meso pores (the pores in the Pt-carbon particles, also referred to as the primary pores).

𝑉 𝐼 = 𝑉 𝐼𝑅 + 𝑉 𝐼𝑀 (A-3)

Where 𝑉 𝐼𝑅 is the volume of reconstructed ionomer (the identified coating and filament ionomer) and 𝑉 𝐼𝑀 is the volume of missing ionomer.

We have also the following relationships,

Where 𝜀 is the porosity identified from the FIB-SEM image (𝜀 ≈ 0.55 for the whole image).

Where 𝑉 𝑆𝑆𝑃 is the volume of the segmented solid phase in the FIB-SEM image. For the considered block of size 250 × 250 × 250 nm 3 and porosity 0.55, 𝑉 𝑆𝑆𝑃 and 𝑉 𝐼𝑅 are known.

We also known the expected I/C ratio (0.8), denoted by

where 𝜌 𝐼 is the ionomer density (2 g/cm 3 ) and 𝜌 𝐶 is the carbon density (1.8 g/cm 3 ) and the platinum carbon ratio (0.462),

Where 𝜌 𝑃𝑡 is the platinum density (21.45 g/cm 3 The next step is to express 𝑉 𝑣𝑜𝑖𝑑-𝑚𝑖𝑐𝑟𝑜 as a function of 𝑉 𝐶 Let 𝜀 𝑚𝑖𝑐𝑟𝑜 denotes the porosity of the carbon particles (𝜀 𝑚𝑖𝑐𝑟𝑜 ≈ 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

Where 𝑉 𝑃𝑡-𝑐𝑎𝑟𝑏𝑜𝑛-𝑝𝑎𝑟𝑡𝑖𝑐𝑙𝑒𝑠 is the volume of the Pt-carbon-particles. Combining Eqs. (A-14)

and (A-15) yields,

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

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

The sole unknown in Eq.(A-18) is 𝑉 𝐶 . To obtain 𝑉 𝐶 we start from Eq. Which leads to express 𝑉 𝐶 as

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

Chapter 5

Reconstruction of contact angle distribution in three-dimensional digital images of catalyst layer microstructure of PEM fuel cells

Abstract:

The wettability properties of the porous components in a proton exchange membrane fuel cell (PEMFC) are a key factor in relation with the water management issue and a key information in order to predict the water retention curve, a crucial input data for the macrohomogeneous models. In this work, we propose an approach to construct the contact angle distribution in three-dimensional digital images of the cathode catalyst layer. Various reconstruction scales are considered, from the platinum particle scale to the pore scale and the scale of sub-blocks in the CCL. This leads to obtain effective contact angles distribution depending on the scale of interest. The approach also allows to investigate the link between the pore size and the effective contact angle at the pore scale. The water retention curve is obtained from the combined information of the pore size and the effective contact angle for each pore.

Keywords: Ionomer distribution, contact angle, PEM fuel cell, Catalyst layer

provides a filtered image since the voxel size (5 nm in our case) is greater than the size of the smallest elements (which can be as small as 1 nm in radius [13], platinum particles ~ 4 nm on average [30]) in the CCL. Also, as discussed in detail in [31] and recalled in the next section, the used FIB-SEM procedure does not allow to distinguish the ionomer from the carbon in the obtained segmented solid phase. The ionomer, and in particular the ionomer coating, is reconstructed in the FIB-SEM segmented 3D image via numerical procedures as detailed in [31]. In spite of these shortcomings, we believe that valuable information on the wettability properties of the considered CCL is obtained. Also, as mentioned before, the main objective is to present a methodology which can be applied to any 3D image of CCL microstructure.

The paper is organized as follows. Information on the CCL and the CCL 3D image obtained via FIB-SEM and the ionomer reconstruction procedure are recalled in Section 2. Section 3 presents the method to determine the contact angle distribution at various scales of interest as well as illustrative results. The water retention curve is determined in Section 4 via a simplified approach and compared to results from literature. A discussion is proposed in Section 5. Section 6 consists of the main conclusions of the study.

Cathode catalyst layer microstructure

The CCL considered in the present study is a commercially available catalyst layer (TEC10E50E; Tanaka Kikinzoku Kogyo K.K., Japan) with a Pt content of 46.2 wt%, I/C ratio of 0.8 and Ketjen Black high surface area (HSA) carbon. As described in [31], a 10.9 µm ×13.2 µm ×5.2 µm digital image of the CCL microstructure is obtained by FIB-SEM with a voxel size of 5 × 5 × 5 nm 3 . Details on the image processing used to obtain the segmented image can be found in [31]. As discussed in [31], the ionomer cannot be distinguished from the carbon phase in the segmented image. Also, the platinum particles cannot be distinguished since their size (around 4 nm and below) is smaller than the voxel size (5 nm). As a result, both the ionomer and the platinum particles are reconstructed in the segmented image using numerical techniques. The method for distributing the Pt particles is recalled in Section 3. In this section, the main features of the ionomer reconstruction are recalled. The details can be found in [31].

Two forms of ionomer, namely the ionomer filaments and the coating ionomer, are distinguished. The ionomer filaments are ionomer protrusions partially clogging the pores.

The ionomer filaments are identified by applying a filtering procedure to the distribution of curve in systems with surfaces of highly heterogeneous wettability is still a considerable challenge. One way of avoiding these difficulties is to introduce the concept of effective contact angle. The effective contact angle [66,67], also referred to as the apparent contact angle [67], describes the wettability property on average, i.e. over sub-regions of the pore space of size larger that the voxel size in Fig. 6. The use of the effective contact angle concept is deemed to be sufficient for determining the water retention curve as considered in [13].

In [13], the considered contact angle is the effective contact angle at the pore scale, and not at the scale of each voxel as illustrated in Fig. 6. The average pore size in our CCL determined in [31] by application of the maximum inscribed sphere algorithm available in PoreSpy [59] is 60 nm, thus much greater than the voxel size of 0.5 nm in Fig. 6. The next step is therefore to determine the effective contact angle distributions. To get insights into the scale dependence of the effective contact angle, we consider, as mentioned before, three different scales. The first scale is still a voxel-scale but with a voxel size greater than the size of 0.5 nm illustrated in Fig. 6. The voxel size of 2.5 × 2.5 × 2.5 nm 3 used in the ionomer coating reconstruction [31] is considered. The effective contact angle at the scale of the 2.5 × 2.5 × 2.5 nm 3 voxels is referred to as the "local" effective contact angle. The second scale one is a cubic block scale of 250 nm, which is the first scale considered for the characterization of the CCL macroporosity in [31] or CCL effective transport properties in [69]. The third scale is the pore scale where the pores are obtained by partitioning the pore space into pores according to a watershed algorithm [59]. As we will see, this latter scale is in fact generally smaller than the block scale of 250 nm

where 𝐴 𝐶 is the surface where the carbon phase is in direct contact with the pores, 𝐴 𝑖𝑓 is the surface of the ionomer filaments, 𝐴 𝑖𝑐-2.5 is the surface of the coating ionomer of thickness 2.5

Chapter 6

Conclusions and prospects

In this study, we have extensively investigated two crucial layers within the PEM fuel cell, employing a combination of imaging data and numerical simulations. This approach has enabled us to quantify various structural and physical properties of these layers.

We began by computing the oxygen diffusion coefficient in the microporous layer and integrating it into the macrostructure of the gas diffusion layer. In Chapter 2, we evaluated the overall resistance within this layer, and a similar computation was performed for the catalyst layer in Chapter 4, focusing solely on oxygen diffusion within the pore space. By summarizing the oxygen diffusion resistance across these distinct layers, we are able to characterize the overall oxygen diffusion resistance as a composite effect of these layers. This presented graphic illustrates the contribution of each layer to this resistance:

Our findings suggest that the primary source of oxygen transport resistance within the pores on the cathode side can be primarily attributed to the gas diffusion medium. These results