Modelling of solid oxide cell oxygen electrodes

Numerical models are versatile tools to study and predict efficiently the performance of solid oxide cells (SOCs) according to their microstructure and composition. As the main contribution to the cell polarisation is due to the oxygen electrode, a large part of the proposed models has been focused on this electrode. Electrode modelling aims to improve the SOCs performance by serving as a guide for the microstructural optimisation, and helps to better understand the electrochemical reaction mechanisms. For studying the electrode microstructure, three categories of models can be distinguished: homogenised models, simplified geometry based models, and reconstructed microstructure based models. Most models are based on continuum physics, while elementary kinetic models have been developed more recently. This article presents a review of the existing SOCs models for the oxygen electrode. As a perspective, the current challenges of electrode modelling are discussed in views of a better prediction of the performance and durability, and more specifically for the case of thin-film SOCs.


Introduction
For the ecological transition of our society, it becomes nowadays essential to develop and commercialize cost-effective and environmental-friendly electricity. Among the variety of energy conversion devices, solid oxide fuel cells (SOFCs) convert efficiently the chemical energy fuels (hydrogen or hydrocarbons) into electrical energy through electrochemical reactions. In the reversible operation mode (electrolysis mode), solid oxide electrolysis cells (SOECs) consume electricity to produce hydrogen by reducing H 2 O. With the development of the hydrogen-based mobility, SOECs are becoming an efficient way to produce 'green hydrogen' when connected to renewable electricity sources, such as solar panels or wind turbines [1]. Reversible SOCs are very convenient for electricity storage applications as they allow smoothing the intermittence of renewable energy sources. Thus, hydrogen is produced when the energy supply is higher than the electricity demand and inversely [2].
Commercial SOFCs operate generally at high temperatures (800 • C-900 • C), as the transport properties of charged species (O 2-) in the electrode and electrolyte materials are thermally activated. However working at these high temperatures results in high fabrication costs, and long-term degradation in the cell performance. Thus, the main challenge is to decrease the operating temperature of SOFCs to gain in durability and cost while maintaining relatively good performance [3,4]. The second-generation of SOFCs (intermediate temperature-SOFCs or IT-SOFCs) are characterized by a thick anode that provides mechanical support and a thinner electrolyte (∼10 µm), which allows operating at lower temperatures (∼650 • C-800 • C). Nowadays, efforts are made to reduce the operating temperature below 650 • C, while maintaining acceptable performance (above 1 W cm −2 ) [5]. To reduce the electrodes and electrolyte polarization resistances at low temperatures, two strategies are generally followed. On the one hand, new electrode and electrolyte materials with better catalytic and transport properties are studied. On the other hand, the electrode microstructure is optimised to reduce the total polarisation resistance.
Traditionally, in high temperature SOFC, the electrolyte is made of yttria-stabilized zirconia (YSZ) while a composite of lanthanum strontium manganite and YSZ (LSM-YSZ) is used as oxygen electrode. In this electrode composite, the charge transfer reactions are limited to the line where the gas phase, the electronic conductor (EC) LSM phase and the ionic conductor (IC) YSZ phase are in contact (known as the triple phase boundary, TPBs line). In order to extend the electrochemically active region to the surface of the electrode particles, mixed ionic-electronic conducting (MIEC) materials are now commonly used at intermediate working temperatures. Among MIEC electrodes, mainly two types of mechanisms exist for the oxygen transport by solid-state diffusion. For simple perovskites, such as lanthanum strontium cobalt (LaSrCoO 3-δ ), the charge carriers for the ionic conductivity is ensured by oxygen vacancies [6], whereas for Ruddlesden-Popper oxides, the ionic conductivity is provided by oxygen interstitials (e.g. La 2 NiO 4+δ ) [7,8], through an interstitially diffusion mechanism. Even if the use of MIEC materials in SOCs allows increasing the electrode performance, their chemical instability and reactivity with the electrolyte, along with the thermal expansion coefficient mismatch with YSZ, remain an issue [9]. To overcome chemical and mechanical electrode degradation problems, a third category of electrodes based on composite formed by the mixture of a MIEC with a pure IC (e.g. LSCF/CGO) provide a promising solution. Indeed, it has been shown that the cathodic polarization resistance can be significantly reduced when the pure MIEC electrode is replaced by a MIEC/IC composite. In particular, several studies have shown that the addition of 30-60 wt.% CGO in lanthanum strontium cobalt ferrite (LSCF) electrode leads to a significant improvement (decrease) of the electrode polarisation resistance [10][11][12][13][14]. This improvement is explained by the role of the IC with higher oxygen diffusivity than the MIEC, which acts as a new pathway for oxygen ion transport in the electrode, leading to an increase of the effective electrode thickness [5]. On the other hand, the high density of TPBs in the composite electrode also participates to enhance the performance [15].
The performance of a given electrode is commonly measured by electrochemical impedance spectroscopy (EIS), and plotted in the form of the area specific resistance (ASR) against the inverse of temperature. These measurements are typically performed on symmetrical cells by measuring the impedance response under a small sinusoidal perturbation of the current or voltage. By varying the frequency over several decades, typically in the range from the mHz to the MHz, an electrochemical impedance spectrum can be obtained, which is usually presented as a Nyquist or Bode plot [16]. The analysis of these spectra under varying experimental conditions can help reveal the limiting electrochemical or physical processes taking place at the electrode. By using a reference electrode EIS can be acquired not only at the open circuit potential (OCP) but also under polarization. These characterizations can be complemented by cyclic voltammetry or polarization curve (potential vs current density plot) measurements.
Modelling is a relevant tool that allows to understand and interpret the physical and electrochemical phenomena observed experimentally. Beyond that, a predictive model should enable to optimize the studied system, with inexpensive cost compared to an experimentally based trial-error method.
For any modelling approach, the general process can be divided into three main parts: the choice of the approach, the implementation of a numerical method or the derivation of an analytical solution, and the validation. Firstly, the choice of the approach is usually imposed by the problematic of each specific study and the objectives: what is being modelled? At which scale? In addition, the set of hypotheses must be carefully chosen according to the geometry and the physics of the studied system: what is the electrode microstructure? Which are the relevant electrochemical processes? Relevant hypotheses must ensure a reasonable balance between the accuracy of the results and the model complexity. Furthermore, the physics and boundary conditions must be clearly identified. Secondly, the model is solved through either an analytical approach or numerical methods using either a commercial software (e.g. Matlab, COMSOL Multiphysics [21][22][23], ANSYS [24]) or within the framework of an in-house or open-source code (e.g. Parcell3D [25,26]). The geometry and the physical equations are defined and implemented according to the hypotheses previously made. After solving the governing equations, a validation step is crucial to ensure good confidence in the results. Figure 1 illustrates the different relevant modelling scales for the optimisation of a complete SOC system. First, the electrochemistry and transport mechanisms can be modelled at the scale of the electrode microstructure. Sometimes called 'micro-modelling' , as opposed to 'macro-modelling' for the larger system models, electrode modelling is focused on the materials performance, the influence of the microstructure and the reaction mechanisms. It is worth noting that a large part of the developed micro-models concerns the oxygen electrode, which is known to contribute significantly to the cell polarization. At the next level, the single repeat unit (SRU) models, which take into account a single cell sandwich between two interconnects, can predict the cell polarisation curve or temperature gradients arising in operation. These models can be used to optimise the cell geometry or the design of the gas channels. They can be extended to simulate the distribution of gas composition, partial pressures, currents or temperatures within a stack, which includes several SRUs. Lastly, the system level takes into account external components such as turbines, compressors, condensers, and heaters, to optimise the system design. The system, stack and Modelling of SOCs from the electrode to the system scale. Assumptions, physics and state variables are adapted to the scale of the study. Typical state variables are: pO 2 the oxygen gas phase partial pressure, cO the oxygen concentration in electrode material, xO 2 and xH 2 the oxygen and hydrogen molar fraction in gas channels, j the current density, Ve the electrical potential, T the cell temperature, P the pressure, U the gas velocity, η the system efficiency. Adapted from [17], Copyright (2020), with permission from Elsevier. Adapted from [18] SRU models have been extensively discussed in previous reviews (e.g. [27][28][29][30][31] in SOFC mode and [32,33] in SOEC mode). However, to the best of our knowledge, no recent review is dedicated to the models at the microstructure scale, models which take into account multi-component mass transport as well as the electrochemical behaviour of the electrode material.
The present article reviews the existing microscale models proposed for the SOCs oxygen electrode. It aims to highlight how electrode modelling can help to gain an in-depth understanding of two main issues: (i) the relation between microstructure and electrode polarization resistance and (ii) the oxygen reduction reaction (ORR) mechanisms occurring in MIEC electrodes.
This review revises the various modelling approaches, which have been implemented to study the influence of the oxygen electrode microstructure. The second part deals with the elementary kinetic models that aim to investigate the reaction mechanisms occurring under anodic and cathodic polarization. Finally, the last section examines the modelling of very thin electrodes for thin film or micro SOC (µSOC) applications.

SOC modelling at the scale of the electrodes microstructure
This section deals with the models that focus on the influence of the microstructural parameters on the electrode performance. SOFC electrodes have been extensively studied and modelled in order to optimize their microstructure. These models can be divided into three categories depending on how the microstructure is implemented. Their goal is generally to predict the impedance response of an electrolyte-electrode system.
The first category of models involves homogenized or analytical models. The microstructure is represented in terms of macrohomogenous parameters such as porosity, volume-specific surface area and tortuosity. It generally provides good results as long as the microstructural parameters are known. However, these models are inadequate if the size of the active region is in the range of the microstructural features (particle or pores). The homogenized models do not take into account local inhomogeneities, current constriction, graded structures or anisotropy in the structure.
The second category of models involves simplified microstructures. The microstructure is represented by regular geometrical objects such as spheres, cubes or rods. Contrary to the homogenized models, simplified microstructure models can capture some effects of 3D microstructures. However, they assume that the complex microstructure can be modelled by monodispersed simple geometrical objects homogeneously distributed in space. Some microstructural details of (e.g. rough interfacial contact areas) are thus neglected. To overcome this issue, more complex microstructural models have been recently proposed to be representative of the real electrode morphology [34][35][36]. For the third category, models are built considering the real microstructure as the geometry. This approach takes into account all the microstructural details and is able to describe precisely the shape and size of the particles, the current constriction gradients or anisotropy inside the microstructure. The microstructure can be reconstructed by focused ion beam-scanning electron microscopy (FIB-SEM) [37] or x-ray tomography techniques [38].
This review is focused on models based on the resolution of continuum-physics conservation equations (such as Fick laws of diffusion, conservation of the electrochemical gradient, current density conservation). The solutions can be analytical after some drastic simplification. Otherwise, the equations need to be solved via conventional numerical methods e.g. finite element method (FEM), finite volume method, and lattice Boltzmann method (LBM). Other classical modelling approaches are based on the use of a resistance network or equivalent circuit modelling (e.g. transmission line models [40][41][42], finite length Gerischer [43], finite length Warburg [44]). Equivalent circuit models generally do not take explicitly into account microstructural parameters. It is worth mentioning that some authors, starting from continuum physics and applying discretisation techniques, derived equivalent circuit models involving some microstructural parameters [45,46].

Homogenised electrodes
When it comes to study the ORR mechanisms in a MIEC electrode, the Butler-Volmer equation may be insufficient, as it describes only the charge-transfer occurring in the electrode, while other surface and bulk mechanisms are not taken into account. Adler et al constructed a model (figure 2) with the following macroscopically defined processes: O 2 diffusion in the gas phase, O 2 exchange through the MIEC surface, oxygen transport in the bulk by solid-state diffusion, and charge-transfer at the electrolyte/electrode interface [47]. In this contribution, only the main analytical equations of the Adler-Lane-Steele (ALS) model are summarized. The model is based on a porous electrode with a known microstructure. The electrode microstructural parameters needed for the model are the porosity fraction ε, the surface area density a and the solid-phase tortuosity factor τ .
In the case of a thick porous semi-infinite macro-homogeneous electrode and in the absence of limitations due to gas phase diffusion, Z chem is the convoluted impedance of the surface exchange and solid-state diffusion in the electrode material at OCP and described by the following equations: where R is the ideal gas constant, T the temperature, F Faraday's constant, a the specific surface area, D chem the chemical diffusion coefficient, k chem the chemical surface exchange coefficient, c O,eq the equilibrium oxygen concentration in the MIEC (i.e. in open circuit condition) and γ O the thermodynamic factor with p O2 the partial pressure of gaseous oxygen). Originally, Adler et al used a formalism based on oxygen vacancies; however, the equation is also valid for materials with oxygen interstitial point defects [48][49][50].
The hypothesis of a semi-infinite electrode is valid if the cathode thickness L is at least three times higher than the materials characteristic length l δ (or 'penetration length'), L > 3l δ with: As both ionic conduction and surface exchange impact the penetration length, Adler et al showed using their model that in the case of L ≫ 10µm, most of the electrode thickness is electrochemically inactive for typical MIEC materials considered. In their seminal work, the authors contributed to move from the traditional way of considering the solid state electrochemical polarization based on charge-transfer at the electrode/electrolyte interface, towards a macroscopic thermodynamic gradient of chemical species view, where the overall mechanism is limited by the chemical reaction and solid-state diffusion occurring several microns away from the charge-transfer interface [47,51].
Equation (1) is mathematically equivalent to a Gerischer impedance, so that the model fits successfully the impedance diagram of a porous MIEC electrode obtained by EIS. The ALS model-as it is a well-detailed model-has been the starting point of many studies, some of which deal with its limitations. As an example, Mortensen et al extended the ALS model to MIEC-IC composite cathodes [52,53] and used it to fit impedance spectra for a (La 0.6 Sr 0.4 ) 0.99 CoO 3−δ /Ce 0.9 Gd 0.1 O 1.95 composite (LSC40/CGO10) under various temperature and oxygen partial pressures.
In addition, Nielsen et al derived a 1D analytical model for pure MIECs and composite cathodes in the case of finite-length electrodes (i.e. when L < l δ ). In their work, the finite-length-Gerischer impedance is introduced, which is more general than the Gerischer impedance, as it is not restricted to thick electrodes [41]. Figure 3 illustrates the change of shape of the impedance spectra from a Gerischer-type spectrum (for L/λ o = 10) to a pure semi-circle curve (for L/λ o = 0.2) with decreasing electrode thickness.
When the electrode thickness is high compared to the characteristic length, the diffusion process in the bulk is rate determining, whereas for thin cathodes, the polarization resistance is determined by the surface exchange process only. Then, for thin cathodes (L < l δ ), the impedance takes the following expression: This expression has the form of the impedance of a RC equivalent circuit and results in a perfect semi-circle in a Nyquist plot.
To take into account gas diffusion and its oxygen partial pressure dependence in the impedance calculation, 1D numerical models have been developed by Deseure et al to model DC and AC impedance of a homogeneous porous LSM-YSZ composite electrode [54].
It can be also noted that the ALS model has been recently extended to express the semi-analytical solutions in terms of linear sweep voltammetry and cyclic voltammetry of porous MIEC electrodes [56].

Electrode models based on a simplified microstructure geometry
The second type of models for a more detailed description of the electrode microstructure are based on the construction of a simplified geometry using spheres, cubes or rods. The size and the percolation of particles, the active surface, and the space between particles are chosen so that they correspond to the effective microstructural parameters. As the current density is inhomogeneous in the particles, current constriction can increase the polarization resistance. The simplified geometry-based models are especially relevant to study the effect of porosity, volume fraction or particle size, especially for composite electrodes where percolation and connectivity can induce limitations on the mass transport. These models involve a three-dimensional geometry, so no analytical solution can be derived, and subsequently they need to be solved numerically. Fleig was one of the first to study current density inhomogeneities in electrode materials by using a 2D FEM model [57]. While the ALS model is based on conservation of species flux (vacancies or oxygen interstitials), Fleig's model is formulated in terms of electrochemical potential. Both approaches are equivalent as the electrochemical potential is linked to the species concentration. The 2D current distribution is simulated in the MIEC cathode and the case where the bulk diffusion is rate limiting is presented [57]. Fleig and Maier extended these calculations to cases in which the oxygen incorporation into the cathode significantly influences the overall reaction rate. They studied the relation between geometry and polarization resistance for different material parameters [58]. In their model, the first electrode geometry consist of pillars of cylindrical-shaped MIEC, whereas the second geometry consists of a stack of small and large (in diameter) cylinders to create a more realistic topology [58].
Several years after the publication of the ALS model, Adler's group addressed some of the deficiencies of the 1D ALS model with an extended model [55,59]. Firstly, the model includes surface diffusion. Secondly, the 3D diffusion within the particles is implemented by using MIEC cylinders, where the geometry is chosen to match with the porosity and the surface area of the microstructure. The choice of using this 'pseudo-particles' model, as shown figure 4 is based on a compromise between accuracy and simplicity of the model. The authors provided an analytical solution for a MIEC column, in which the prediction of their model is equivalent to the ALS one when the microstructure feature (i.e. the column) is small enough compared to the penetration length (l δ ). They acquired impedance data on porous lanthanum strontium cobalt (LaSrCoO 3-δ ) electrodes over a much wider range of conditions than previously reported, and drew a more complete comparison to theoretical predictions. Besides, the extended model takes into account surface diffusion, as well as 3D diffusion within particles near the interface when the characteristic utilisation length l δ is lower than the size of individual morphological features d, l δ < d [59].
A more complex 3D model was developed by Rüger et al by using a cubic-based geometry [23]. The 3D model is based on the description of the oxygen concentration in the MIEC. The processes of gas diffusion, surface exchange, bulk diffusion and charge transfer are implemented in a 3D FEM commercial software (COMSOL Multiphysics), as shown in figure 5. Pore and MIEC cubes are randomly distributed to obtain a homogeneous structure as depicted in figure 6. The goal of the model is to calculate the area specific resistance ASR cat of a MIEC cathode as a function of the material's parameters and its microstructure. The chemical diffusion coefficient (D chem ) and the chemical surface exchange coefficient (k chem ) are taken from the literature. The oxygen ion equilibrium concentration c O,eq was evaluated as a function of oxygen partial pressure for each temperature. The use of non-effective parameters allows implementing complex electrode microstructures [23].
The gas diffusion of oxygen gas molecules in the MIEC pore is modelled using the dusty-gas-model. This model is commonly used to describe the molecular diffusion of the binary mixture of oxygen and nitrogen and the Knudsen diffusion related to the interaction of the oxygen molecules with the pore wall: − → J O2 is the oxygen molecule flux density, p the total pressure, M O2 and M N2 are the molecular mass of oxygen and nitrogen respectively, x O2 is the oxygen molar fraction, D O2N2 and D K O2 are the binary and the Knudsen diffusion coefficients, respectively. The surface exchange at the gas/electrode interface is driven by . Due to the assumption of constant electric field in the MIEC, the driving force of the ionic bulk diffusion results only from the concentration gradient. Thus, the transport of oxygen ions is described by Fick's first and second laws: In the electrolyte, a constant chemical potential is assumed. Then ionic current − → j ion is determined only by the gradient in the electrical potential: where σ El is the ionic conductivity of the electrolyte and Φ El the electrical potential. The ionic current conservation law is then solved: The model allows studying the optimal electrode thickness according to the penetration depth ( figure 6). The impact of the particle size on the ASR was investigated while the volumetric composition of the MIEC and the porosity remain constant. It was shown that the ASR decreases by a factor of 2 when the particle size decreases from 1 µm to 250 nm.
As a continuation of Rüger's model, Häffelin et al implemented a time-dependent 3D impedance model for a MIEC SOFC cathode that considers the complex coupling of gas diffusion, surface exchange, ionic bulk diffusion and electrolyte conductivity. The goal was to carry out a detailed analysis of the formation of a Gerischer-type impedance. The surface exchange reactions dominate the low-frequency region, whereas the 45 • Gerischer ramp is related to the ionic diffusion. The model also manages to separate and quantify the contribution of the gas diffusion in the cathode pores. For instance, it was shown that, for a typical cathode with 30% of porosity at an oxygen partial pressure of 0.21 atm, the gas diffusion accounts for 2% of the total polarization, whereas a depletion of oxygen to 0.01 atm multiplies the gas diffusion impact by almost a factor 20 [25,62].
More recently, Celikbilek et al adapted Rüger's model to LSCF and LSCF-CGO composite SOFC cathodes with a multiscale microstructure [63]. The authors' goal was to compare the performance of both composition, while taking into account the real microstructure of the electrode layers. Both compositions show an architecture containing macropore channels and nanoporous pillars, even though the composite electrode shows lower nanoporosity. Macrochannels are designed in the model geometry. A source term is included in the homogeneous material volume to take into account the oxygen incorporation through the nanopores.
Some of these works include comparison with experimental data showing that even a simplified geometry model provides a good estimation of the electrode performance [59,64], and can be used to study the influence of the microstructure on the performance. Besides, no fitting procedures are used for the model input parameters as they are obtained by independent experimental characterization. In this case, it is crucial to characterize the bulk material properties, such as the D chem , k chem and the non-stoichiometric oxygen-as well as their variation with temperature and oxygen partial pressure.

Reconstructed microstructure electrodes
The last category of model consists to use for the simulated domain a reconstructed electrode. This method is the most accurate to study the microstructure impact on the performance, but it is also the most challenging to implement for various reasons, such as the high computational cost.
The 3D reconstruction of the real microstructure of the electrode became possible with the development of the dual-beam FIB-SEM technique. The method is briefly described in the following. The sample is milled and tilted so that a cross-sectional SEM image can be taken. Then, the area observed is slightly milled in the z direction by the FIB (around 10-50 nm are removed) so that a new x-y observation surface is exposed for SEM imaging (see figure 7). The process is repeated automatically to yield a series of 2D SEM images that will be stacked to reconstruct the final 3D image. The method was used to study, for instance, the correlation between microstructural and electrochemical characteristics during redox cycles for a Ni-YSZ SOFCs anode [65], or to compare two sintering methods used to manufacture anode supported SOFCs [66]. This technique is also widely used, especially in composite electrodes, to extract real microstructural parameters, such as porosity, tortuosity, active surface area, volume fraction, or TPB length [67][68][69]. As described in [68], the main advantage of the FIB-SEM technique is the very high quality of the collected volumetric data. However, the simulation of the electrochemical impedance on reconstructed geometries appears challenging, as it requires a large amount of computational time due to the need of a large volume combined with high resolution.
In order to study more accurately the interaction of microstructure and performance, Rüger et al showed for the first time the possibility of using as numerical model geometry the reconstructed microstructure [60]. The main concern was the size of the representative volume elements (RVE). The RVE has to be sufficiently large to consider the averaged-properties independent of its location. Also, the 3D FEM implemented in COMSOL Multiphysics showed, at that time, issues in the management of the geometry, due to limitation in the number of voxels, which define the solid material [60]. In the same group, this issue was later overcome by implementing a C++ program, called ParCell3D, which was developed specifically to solve the problem of the amount of available memory [26,70].
Two years later, Matsusaki's group managed to reconstruct and solve the oxygen chemical potential (µ O ) inside a MIEC electrode LSCF by using the LBM [61]. Figure 8 illustrates the oxygen chemical potential gradient (a) and both electronic and ionic current streamlines (b) through the reconstructed the MIEC electrode. The authors showed that the predicted polarisation curve in cathodic polarisation matches with the experimental results. In addition, they demonstrated that their model is an effective tool for investigating local oxygen potential field.
As alternative to the FIB-SEM technique, Lynch et al managed to reconstruct a 3D micro/nano-scale porous structure using x-ray nanotomography. The authors coupled COMSOL Multiphysics and iso2mesh software with custom (Matlab) code for automating recognition, meshing, refinement, equation assignment and solution. Compared to FIB-SEM, x-ray nanotomography has the advantages of being a non-destructive and distortion-free reconstruction method with good phase sensitivity [71]. In this study, the effect of the local nanostructured morphology on the global electrochemical response has been demonstrated using a LSM electrode.
More recently, Hsu et al used a real reconstructed microstructure for modelling a heterogeneous LSM-YSZ composite electrode [72]. In this particular case, four phases are meshed and divided into four subvolumes: the EC (LSM), the IC (YSZ), the pores, and the TPBs. The same model was used to study the effect of the LSM infiltration into a YSZ backbone [73]. As other examples, Yu et al [17] and Goel et al [74] studied the diffusion within a real 3D electrode microstructure reconstructed by FIB-SEM. Their model is based on a complex formalism and the governing equations are solved by the smoothed boundary method. As demonstrated by Cooper et al [75], these 3D models have the ability to reproduce experimentally observed departures from ideal impedance spectra shown in figure 3, such as semi-circle deformations and flattening. However, it is important to note that the electrode impedance spectra is also significantly influenced by the underlying reaction mechanisms, which are not yet fully understood and modelled.
Thus, 3D tomography techniques such as the FIB-SEM and x-ray computed tomography, allow studying the optimal electrode microstructure for the best SOFC performance, but also the electrode degradation that can occur at prolonged operating time. While these cutting-edge 3D characterization techniques provide very accurate 3D microstructures, they remain very costly and time-consuming. Alternatively, numerical simulations are also a good solution to investigate the influence of the 3D microstructure on the electrochemical performance of SOFC electrodes, without performing any experiments. As an example, Yan et al developed a numerical framework based on successive Discrete Element Method (DEM), kinetic Monte Carlo method, and LBM to predict SOFC performance from the raw powder [76]. Three steps are simulated: (i) the powder forming, (ii) the numerical sintering, and (iii) the electrochemical calculation.

Modelling the oxygen reaction mechanisms in SOC electrodes
One crucial step in the development of high performing SOCs is the understanding of reaction mechanism in the electrodes. The study of the oxygen reaction mechanisms is a complex task as many parameters are involved, such as the microstructure (thickness, specific surface area, porosity, and tortuosity), the transport properties of the materials (ionic and electronic conductivity), the external conditions (temperature, potential, and gas-phase composition) and the electrochemical processes.
In the case where only one rate-determining step is modelled, an analytical approach can be used. On the other hand, if several rate-determining steps are taken into account, a numerical solution of the set of partial differential equations is required. Several studies focus on modelling elementary kinetics to identify the rate determining steps that are linked to the current-voltage characteristics and dependencies of the polarization resistance with the partial pressure.
Adler et al suggested that the oxygen reduction in a MIEC follows two steps [47]. First, the oxygen is reduced chemically at the gas/MIEC interface. Then, it is transported through the MIEC by solid-state diffusion. This mechanism is called the 'bulk path' by opposition of the 'surface path' that occurs for electrodes made of poorly ionic conductivity materials (e.g. LSM). In this particular case, the charge transfer occurs only at the TPB.
Several authors combine the experimental impedance measurements with 1D equivalent circuit models to identify more precisely the oxygen reduction mechanisms. Using a LSCF/CGO composite electrode, Grunbaum et al suggested that the high frequencies contribution is due to the oxygen diffusion in the bulk of the MIEC material, whereas, the low frequencies contribution is due to the dissociative oxygen adsorption in the MIEC surface and/or the oxygen gas diffusion in the pores of the electrode [77]. The low frequencies region strongly depends on the oxygen partial pressure and gas carrier. The simulation and study of impedance spectra allows for a better insight into the reaction mechanisms.

Elementary kinetic models for SOFCs
The emergence of elementary kinetic models was initiated by the theoretical work of Svensson et al followed a few years later by that of Deseure et al on the mathematical modelling of oxygen exchange and transport in IC-EC composite electrodes [78,79].
Independently, Bessler et al identified two limitations of the common approach for modelling a SOFC based on the subtraction of overpotential from the Nernst equilibrium potential [80]. Most of the models apply the Butler-Volmer equations to describe current/overpotential relationships. However, these equations do not take into account complex electrochemical processes at the three-phase boundaries. The authors proposed a new modelling approach, where the electrochemistry is described by elementary-kinetic reactions, especially for coupled charge-transfer and surface chemistry. In their new framework for electrochemical simulation, the Nernst equation is not used. The charge-transfer reaction is given by an elementary-kinetic equation. Physically, the potential step at the electrode/electrolyte interface originates from electrical double layer, which is formed of a space-charge region of excess ionic and electronic charge [80]. For the mass transport, the gas flow (1-10 mm scale) is modelled using the Navier-Stokes conservation equations, a continuum approach is used for the charge transport and mass transport in the pores (1 mm-1 µm scale), and the surface diffusion towards the three-phase boundaries (1 nm-1 µm scale) is modelled using Fick's law [80]. The implementation allows the calculation in steady state and transient regime to simulate polarization curves and impedance spectra, respectively.
As an extension of Bessler's model, Yurkiv et al studied the impedance of the surface double layer at the LSCF surface of a LSCF/CGO composite electrode represented in a continuum approach using effective mass and charge transport coefficients. The two interfacial double layers are represented in figure 9. In contrast to the ALS model, ion diffusion is assumed not to be rate limiting and is described by a simple resistance [81]. The analysis of the simulated impedance spectra identified three distinct features depending on the frequency: (i) gas conversion at the oxygen electrode (low frequencies), (ii) electrochemical oxygen reduction on the MIEC surface (intermediate frequencies) and (iii) charge transfer of the double negatively charged oxygen at the two-phase boundary between the MIEC and IC (high frequencies). These observations agree with the experimental impedance measurements carried out by Grunbaum et al [77].
The limitations of these models are mostly related to the non-spatially resolved mechanisms. The transport and the oxygen ion concentration change are modelled by 1D conservation equations across the electrode thickness using effective properties and thus not resolved explicitly at the scale of the bulk material.

Elementary kinetic models for reversible-SOCs
As shown in many studies for porous MIEC materials, in fuel cell mode, the kinetics are controlled by the oxygen incorporation, followed by the ionic transfer at the MIEC/IC interface. On the other hand, in electrolysis mode, the direct electrochemical oxidation/reduction at the TPBs also need to be taken into account, especially for composite electrodes. The electrolysis mode has received much less interest compared to the fuel cell mode. Nevertheless, Laurencin et al have developed an elementary kinetic model taking into account both bulk and surface paths, in both cathodic and anodic polarization [82]. The aim of their study  was to analyse and compare the reaction mechanisms of reversible-SOCs composed of a single-phase MIEC and a MIEC/CGO composite. As depicted in figure 10, the reversible electrochemical kinetic reactions rely on (i) the charge transfer through the MIEC/CGO interface, (ii) the production (or consumption) of oxygen adsorbates at the gas/MIEC interface, followed (or preceded) by (iii) the associative desorption (or adsorption) of oxygen. The surface path is modelled by (iv) the direct oxidation (or reduction) at the MIEC/gas/electrolyte TPBs that produces (or consumes) adsorbed oxygen.
The model developed by Laurencin's group allows studying the polarisation impact on the reaction mechanisms in both modes [82]. It is shown that, in cathodic polarization (fuel cell mode), the surface diffusion does not affect the electrode performance, as expected, even though the surface diffusivity is increased. The adsorption and incorporation steps followed by bulk diffusion control the reaction mechanism. At lower oxygen partial pressure, oxygen adsorption becomes rate limiting. In anodic polarisation, the rate limiting step of the overall mechanism changes from an oxidation occurring on the gas/MIEC surface at low polarisation values to the direct oxidation at TPBs at higher polarisation.
In the same group, several successive studies apply this model for various purposes. Hubert et al performed microstructural sensitivity analysis using both LSCF single phase and LSCF/CGO composite electrodes reconstructed by x-ray nanotomography. The LSCF/CGO composite, which is fully controlled by the surface whatever the polarisation [15], was found to have a better performance than the LSCF single-phase electrode, especially in anodic polarisation (electrolysis mode) [84]. In addition, Monaco et al focused on the study of LSCF under anodic polarization to explain the transition detected at low polarization from the bulk to the surface path [85]. The analysis of the local quantities show that the change of mechanism is due to the depletion of oxygen vacancies. Furthermore, the temperature influence on the reactive mechanisms has been explained by the combination of the two following effects. When the temperature increases, the oxygen non-stoichiometry in the MIEC at equilibrium and the reaction kinetics are both enhanced. Then, at higher temperatures and at high anodic polarisation the bulk path is favoured, as the oxygen vacancies are less easily depleted.
More recently, Effori et al investigated the effect of the LSCF operation and decomposition on the electrode response and demonstrated that the model can accurately fit experimental impedance spectra [87]. In another study, Sdanghi et al investigated the rate-determining steps of a Ruddlesden-Popper type oxide with oxygen interstitials (La 2 NiO 4 ) in both fuel -cell and electrolysis mode [50]. It can be noted that the elementary model described in [87] was also extended to simulate the cyclic voltammetry response of an LSCF-type electrode [88]. It has been shown that the peaks of the voltammograms are related to the transient change of oxygen concentration in LSCF.

Description of thin film SOCs
Thin film solid oxide fuel cell (or micro SOFC, µ-SOFC) technology is a relatively new and recent field, where thick film (tens to hundreds of µm) and bulk materials give way to materials in thin film form (typically tens to hundreds of nm) [89]. The development of µSOFCs has become possible thanks to the progress in thin film technology (deposition and characterization techniques) and microfabrication of thin-film components relevant to SOFCs [90]. With the growing interest of Internet-of-Thing (IoT), µSOFCs are designed to complement and/or replace lithium batteries in electronic devices for health care, automated production lines or smart cities. The fuel flexibility, the high energy density, and the potential to provide continuous power are three attractive characteristics of µSOFCs. To reduce the operating temperature below 500 • C, the electrolyte is much thinner than that of conventional SOFCs (<1 µm). Micro-cells are designed and micro-fabricated on silicon-based substrates using microelectromechanical systems (MEMS) technology, as it allows large-scale and low-cost production of high quality devices. One promising approach is the freestanding electrolyte membrane ( figure 11). This method was used in Tarancón's team by Garbayo et al to fabricate large-area µSOFCs using YSZ electrolyte and LSC/CGO electrodes [86].

Modelling thin film SOC electrodes
Models for micro-SOFCs are scarce. Both Saranya et al and Navickas et al [24,91] independently developed a specific modelling approach to study the grain boundary effects on the oxygen reduction and diffusion kinetics. Sr-doped LaMnO 3 (LSM) was studied as thin film electrode material. By matching simple finite element model results with the measurements obtained by 18 O tracer exchange experiments and subsequently analysed by secondary ion mass spectrometry (SIMS) it was shown that grain boundaries facilitate both the oxygen diffusion and oxygen exchange [24,91]. Figure 12 shows that both diffusion coefficient and surface exchange coefficient must be highly enhanced in the grain boundary to fit the experimental results. A similar model was later extended by Tarancon's group to the study of other perovskite cathode compositions [92].
Using a similar approach, the influence of dislocations on the oxygen ion transport in epitaxial LSM thin films was also investigated by Fleig's group [93]. Depending on the perovskite substrate, the lattice mismatch between the film and the substrate is known to induce compressive or tensile in-plane strain in the LSM layer grown on top. The FEM simulations (figure 13) fit well the experimental data and show that the oxygen ion diffusion along dislocations is about 2-3 orders of magnitude faster than that in the bulk in LSM [93].
At the sub-micrometre scale, the electrode thickness is below the penetration length of the material-or at least close to it. Due to the very low thickness of the electrode, the diffusion of oxygen ions in the bulk material is sufficiently high to be non-limiting. The surface exchange processes become thus rate determining. Then, in the paradigm of very thin film electrodes, the specific surface area becomes a crucial parameter that needs to be enhanced to facilitate the oxygen incorporation. As illustrated for La 2 NiO 4 thin films grown on YSZ electrolyte, Stangl et al were able to improve the electrode performance by largely increasing the specific surface area from a dense layer to a nano-columnar structure [94]. The La 2 NiO 4 films were deposited by pulsed-injection metalorganic chemical vapour deposition and the structure was controlled by the deposition temperature and the film thickness.
Thus, at these low dimensions, microstructural features become important. Therefore, the modelling of µSOC faces the difficulty to implement representative geometries as the microstructural details (pores, grains boundaries) are in the scale of the resolution limit of the common tomography techniques (FIB-SEM or x-ray tomography) [95].   fitted with the model including bulk contribution (orange), the bulk and the dislocations (blue) as long as the addition of the contribution of the in-plane compressively strained interface region (red). Adapted [93]. CC BY 4.0.

Conclusions and future prospects
This review reports the large diversity of SOC models that exist at the scale of the electrode. Each model is designed and implemented to investigate specific aspects of the electrode performance. Regarding the oxygen electrode, modelling has proven to be an efficient tool to understand the oxygen transport mechanisms in the electrode and to study the process-microstructure-performance relationship.
Analytical models provide good approximations when comparing with experimental measurements on homogeneous porous electrodes. Yet, these models are based on very strong hypotheses that should be carefully examined and validated to ensure that the model is fully relevant. Reconstructed models are well adapted to study the effect of the complex electrode microstructure on the oxygen transport. However, these models require heavy computational power and nano-tomography techniques. As an alternative modelling strategy, simplified microstructure models represent a good compromise between accuracy and complexity.
Elementary kinetics models are interesting tools to investigate the oxygen reactive pathway according to the polarization bias, the oxygen partial pressure and the temperature. Nevertheless, for a better understanding of the oxygen transport along with the reaction mechanisms occurring in the electrode material, it would highly relevant to connect these elementary models with computations at the atomic scale. Besides, the simulation of the material decomposition in operation and its impact on the cell performance must also rely on models developed at the atomic scale. In this context, the developments of models coupling elementary kinetic simulations with computations based on the density functional theory should bring new insights in the mechanisms of electrode operation and degradation.
With the recent development of thin film SOCs, the existing models become inadequate at the nanometric scale. As the film thickness becomes smaller than the penetration length, the surface exchange reactions limit the oxygen transport. In addition, microstructural features ignored at higher dimensions such as grains boundaries have proved to play a key role, especially for materials with low oxygen diffusion coefficients.

Data availability statement
No new data were created or analysed in this study.