Modeling the impedance response and steady state behaviour of porous CGO-based MIEC anodes †

Mixed ionic and electronic conducting (MIEC) materials recently gained much interest for use as anodes in solid oxide fuel cell (SOFC) applications. However, many processes in MIEC-based porous anodes are still poorly understood and the appropriate interpretation of corresponding electrochemical impedance spectroscopy (EIS) data is challenging. Therefore, a model which is capable to capture all relevant physico-chemical processes is a crucial prerequisite for systematic materials optimization. In this contribution we present a comprehensive model for MIEC-based anodes providing both the DC-behaviour and the EIS-spectra. The model enables one to distinguish between the impact of the chemical capacitance, the reaction resistance, the gas impedance and the charge transport resistance on the EIS-spectrum and therewith allows its appropriate interpretation for button cell conditions. Typical MIEC-features are studied with the model applied to gadolinium doped ceria (CGO) anodes with diﬀerent microstructures. The results obtained for CGO anodes reveal the spatial distribution of the reaction zone and associated transport distances for the charge carriers and gas species. Moreover, parameter spaces for transport limited and surface reaction limited situations are depicted. By linking bulk material properties, microstructure eﬀects and the cell design with the cell performance, we present a way towards a systematic materials optimization for MIEC-based anodes.


Introduction
Solid oxide fuel cell (SOFC) technology is a promising solution for the efficient use of renewable fuels or natural gas for decentral heat and power generation.The most commonly used anode material in SOFCs is Ni-YSZ.Unfortunately, this material shows various types of degradation including Ni coarsening, carbon coking, sulfur poisoning and mechanical damage caused by redox cycling. 1 Moreover, the electrochemical reaction is bound to the triple phase boundaries, which induces a specific microstructure limitation towards the electrochemical activity in Ni-YSZ cermet anodes.As an alternative anode concept, mixed ionic and electronic conductive (MIEC) materials are drawing much attention. 2An obvious advantage of MIEC materials is the fact that the fuel oxidation reaction can take place on the complete MIEC/pore interface.Furthermore, since it is a mixed conductor, only one contiguous solid phase is needed to enable transport of both charge carriers.In advanced materials concepts, MIECS are also often used as a basis for composite electrodes, whereby the MIEC is combined with an additional contiguous phase of superior electric conductivity.This combination leads to a reduction of the transport resistance.An additional variation of the anode concept is the infiltration of porous MIEC backbone with catalytically active metal particles, which leads to an enhancement of the electrochemical activity as e.g.shown by Primdahl et al. 3,4 and Price et al. 5 The optimization of these different anode concepts is the focus of ongoing research in the fuel cell community.Thereby, the type and chemistry of the MIEC material itself plays an important role.The most relevant MIECs for SOFC anodes are either perovskite-type materials with a wide range of elemental compositions or MIECs with fluorite structure.Ceria is the most prominent fluorite example. 2 The importance of those materials is not restricted to fuel cells.They are also used as selective membranes to produce high purity oxygen, 6 as electrode material in batteries, as sensors e.g. to measure the oxygen partial pressure based on the change in stoichiometry to mention but a few.A more complete overview of MIEC materials and their applications a Zurich University of Applied Sciences, Institute of Computational Physics, Winterthur, Switzerland.E-mail: mame@zhaw.ch;Tel: +41 (0)58 934 70 80 b Hexis AG, Winterthur, Switzerland c Department of Physics, University of Fribourg, Fribourg, Switzerland can be found in the Handbook of Solid State Electrochemistry of Bouwmeester et al. 7 In this study we present a comprehensive numerical model for MIEC-based SOFC anodes on the button cell level, using gadolinium doped ceria (CGO) represented by Ce 0.9 Gd 0.1 O 1.95Àd .CGO is chosen as a reference system because it is widely used in commercial SOFC applications e.g. by Hexis AG 8 and Sunfire GmbH 9 due to its high performance and stability.However, similar or possibly even higher conductivities can be achieved by samarium doped ceria according to Koettgen et al. 10 The importance of reliable models for the optimization of CGObased anodes is emphasized here, because there are still many open questions regarding the underlying reaction mechanism and the rate limiting processes.A detailed physical model for MIEC anodes, such as the one presented in this study enables quantifying the resistive effects from different processes that are otherwise indistinguishably overlapping in EIS measurements.Furthermore, a model-based approach allows performing parametric studies for a wide range of parameters including microstructure and material properties, operating conditions and anode architectures.Such extensive parametric investigations would exceed by far the capacities and capabilities of conventional experimental approaches.
So far, the ceria-based anodes have been studied intensively and in great detail mainly in thin-film symmetric cell setups with transmission line (TL) and finite element (FE) models.TL-Models including corresponding measurements were presented by Lai et al. 11 for SDC15 and by Nenning et al. 2 for SrTi 0.7 Fe 0.3 O 3Àd .FE-models for SDC15 thin film electrodes were presented by Ciucci et al. 12 as a DC and by Chen et al. 13 as an AC (EIS) model (these references are discussed in more detail in Section 2.2).Modeling studies for the case of porous highperformance electrodes with ceria as presented by Kishimoto et al. 14 are less frequent.However, ceria-based anodes are now at the threshold of commercialization. 8,9The requirements of the market are improvements of performance and efficiency, longer lifetime and lower system costs.To reach these goals, a detailed understanding is required of the physico-chemical processes involved in ceria-based porous anodes.This is true for button cell, large area cell and stack levels.In this context it must be emphasized that only a suitable MIEC model -in combination with experiments -is capable to provide a fundamental understanding of the underlying mechanism, which is the basis for purposeful, knowledge-based optimization of different cell concepts and material systems.
In SOFC research, electrochemical impedance spectroscopy (EIS) is an essential characterization tool, which serves as a basis for materials optimization on the electrode, cell and stack levels.However, for CGO-based electrodes, there is no consensus how to interpret the impedance spectra yet.With EIS, the different processes in CGO-based anodes are often not distinguishable from each other because their frequency ranges overlap and because of their interdependencies.Especially for the low frequency arc(s) in impedance spectra of ceria-/CGObased anodes, the resistive contributions from surface reaction resistance/chemical capacitance, gas diffusion impedance and charge transport resistance are often unclear.In addition, the impact of these processes may strongly depend on the cell architecture and type of electrode (e.g.dense thin-film electrode vs. porous anodes in button cells vs. large area cells).Hence, for the same anode material (i.e.CGO) different processes might become rate limiting in dependence of the cell architecture and electrode type, which explains the controversial interpretations of the corresponding impedance spectra in literature.A detailed discussion of this controversy is provided by Marmet et al., 15 who compared literature data of EIS-spectra for dense CGO-based anodes in thin-film electrodes, with the spectra of porous CGO-anodes in button cells and large area cells.The outcome of this comparison can be summarized as follows: In principle, it is well accepted (e.g.ref. 11, 16 and 17) that ceriabased electrodes exhibit a chemical capacitance, showing the ability of ceria to store energy by a change in its local stoichiometry.Hence, this process potentially controls the observed low-frequency arc.In fact, for ceria-based dense thin-film electrodes the low-frequency arc is identified as a surface reaction resistance/chemical capacitance process by Chueh et al. 17 They employed conditions under which other effects like the gas impedance and transport resistance of the charge carriers are negligible.In contrast, for large area cells, there is strong evidence that the low frequency arc must be identified as a gas impedance process.Hence, large area cells exhibit a large gas conversion impedance due to the fuel consumption along the cell.This has been reported for a Hexis short stack with CGO-based anodes by Linder et al. 18 For this case it must be noted that the (low frequency) contribution from the surface reaction resistance/chemical capacitance process is relatively small and therefore hidden within the dominant gas impedance arc.The third case of CGO-based anodes is represented by a button cell setup with an excess fuel flow rate.In this setup the gas impedance is much smaller compared to large area cells, which was reported by Riegraf et al. 19 Hence, whether the low frequency arcs in button cells are dominated by the gas impedance, surface reaction resistance/chemical capacitance process or even by a different process like the transport resistance of the charge carriers, is usually not obvious in the experimental results.Consequently, different interpretations are found in the literature.In order to unravel this controversy, a model for MIEC/ceria-based anodes is presented in this study, which enables simulating the EIS-spectra for button cells, but also for other cell setups (e.g.thin films).Thus, the aim of this model (among others) is to reveal the different low frequency contributions on a quantitative level.
Another point of interest for SOFC developers is the contribution of the transport resistance of the charge carriers to the total area specific resistance (ASR).This contribution is closely linked to the extension of the reaction zone.Therefore, a quantitative description of the reaction zone and associated distribution of charge carriers are considered as important information for a controlled and knowledge-based optimization.Commonly it is suggested in the literature (e.g. by Primdahl et al. 3 and Nakamura et al. 20 ) that for MIEC anodes with a higher electronic than ionic conductivity, the reactive zone of This journal is © the Owner Societies 2021 the electrode will be localized at the electrolyte/anode interface.In some circumstances, however, the electrochemical conversion may spread over the entire anode layer or it can even take place preferentially in a zone adjacent to the current collector.A more detailed discussion of this topic can be found in Section 4.2.With our model operated in DC-mode, we also provide a quantitative description of the location of the reaction zone that gives new insight on the reaction mechanism of MIEC anodes.A reliable quantitative description of the extension of the reaction zone has, to the best of our knowledge, not been reported for porous CGO-based anodes so far.
From the discussion above it follows that reliable interpretations of impedance spectra and systematic optimizations of MIEC-based SOFC anodes require a detailed model, which captures (a) the relevant physico-chemical processes, (b) the associated material laws, (c) the influence of the microstructure and (d) the dependencies on varying operating conditions.Moreover, the possibility to run the model in both, AC and DC modes has great advantages.However, for porous MIEC/ CGO electrodes, such combined AC/DC models are hardly available in the literature.This motivated us to develop a comprehensive model for a CGO-based anode on the button cell level enabling an appropriate interpretation of EIS-spectra and providing a thorough understanding of the complex physico-chemical processes involved.
This paper is structured as follows: In Section 2 we present a thorough description of a 1D MIEC-anode model, which can be operated in DC as well as in AC modes.In a first stage, the model formulation focuses on general aspects of MIEC anodes (i.e. it describes processes that are common to all kinds of different MIEC materials).It also provides a quantitative description of the relationships between effective properties with the corresponding 3D microstructure characteristics.In this way, the 1D model is capable to capture in a reliable way the relevant microstructure effects based on input from 3D analysis.Furthermore, this first stage of model formulation also contains thorough descriptions of the following important processes and phenomena: charge transport in MIEC materials, hydrogen oxidation surface reaction, gas transport and gas concentration impedance.In a second stage, the model is then extended towards a specific MIEC material, which is Gd doped ceria (CGO10).Hence, CGO-specific material laws and phenomena are described, such as crystal chemistry and equilibrium density of charge carriers, nonstoichiometry, electric and ionic conductivities, surface reaction overpotential, exchange current density and chemical capacitance.
In Section 3 the results from DC-and AC-simulations with CGO-based anodes are presented.The DC-simulations provide profiles across the anode layer, which show for example local concentrations of charge carriers, and also the local reaction rates on the CGO surface.The AC-model enables to distinguish specific impedance contributions from different processes such as charge transport, gas transport and hydrogen oxidation surface reaction.In a first step, the simulations are performed for a specific set of parameters describing a standard button cell architecture, standard operating conditions, standard properties of CGO10 and a standard microstructure.In a second step, the DC-and AC-models are then used for parametric sweeps that illustrate e.g. the impact of anode layer thickness and microstructure variations on the cell performance.
The modeling results represent the basis for a detailed understanding of the overall reaction mechanism in CGO anodes, which is discussed in Section 4. The distribution of the reaction zone and associated current densities are discussed in the context of a complex interplay of different processes and associated resistances.The underlying mechanism of coupled processes and rate limiting effects are discussed specifically for charge transport, distribution of the reaction zone, gas impedance, surface reaction and chemical capacitance.Furthermore, the effect of parameter variations describing the anode architecture (e.g.layer thickness) and microstructure properties (e.g.porosity and surface area) are also discussed.The model thus allows us to distinguish and even quantify the overlapping impedance contributions from different processes.
Section 5 concludes with the lessons learned from combined DC-and AC-modeling of CGO anodes, with a special emphasis on the identification of rate limiting processes and their overlapping contributions to impedance spectra for different parameter constellations.The paper extends these aspects by providing an extensive illustration of simulated impedance spectra for different anode/cell architectures (i.e.model electrode, button cell, large are cell), different layer thicknesses and different microstructures (porous-dense-massive, coarse-fine).Together with the model-based identification of rate limiting processes, this overview helps solve the literature controversy about how to interpret the low frequency arcs in impedance spectra for MIEC/CGO-based anodes.

General properties of MIECs
MIECs exhibit both ionic and electronic conductivity.Ion transport occurs normally via interstitial sites or by vacancy motion, both caused by ionic defects.These ionic defects can be formed by thermal excitation, by change of stoichiometry or by doping. 7Electronic conductivity occurs via delocalized states in the conduction or valence bands or by localized states by a thermally assisted hopping mechanism. 7The electronic conductivity is generated by defects as well, but the mechanisms differ from those of the ionic conductivity.The conductivity depends on the mobility m and the number density n of mobile charge carriers.As discussed for the case of CGO in the Section 2.4.1, immobile charge carriers are important to be considered for the electro-neutrality condition but do not contribute to the conductivity.The mobility generally depends on the crystal structure and the temperature T. The number of mobile charge carriers depends on the doping y and/or on the oxygennonstoichiometry d, which is directly linked with oxygen vacancies.Especially when generated by reduction, the oxygennonstoichiometry depends strongly on the environmental conditions (oxygen partial pressure p O 2 ) and temperature.Moreover, the transport of the charges in a MIEC is not simply governed by the drift in the electric field, but as well by diffusion due to concentration gradients of the charge carriers. 11,16Thereby, the mobility and the diffusivity of the charge carriers are related to each other according to the Nernst-Einstein relation eqn (7).
As already mentioned, the simultaneous electronic and ionic conduction enables the electrochemical reaction at the MIEC/ pore interface and is therefore not restricted to the three phase boundaries as in classical cermets.However, this does not mean that the whole MIEC/pore interface is electrochemically active to the same extent.Already Adler et al. 21concluded with their fundamental model for MIEC cathodes that the extension of the reaction zone is limited by the ionic conduction.Nakamura et al. 20 and Primdahl et al. 3 concluded from their experimental studies on MIEC anodes with a higher electronic than ionic conductivity that the active site of the electrode is nearby the electrolyte.
3][24] Because the charges can be stored in the bulk, and not only as surface charges as in a dielectric capacitor, the chemical capacitance is a volumetric effect and can therefore be orders of magnitude higher than dielectric capacitance effects.
In the following, a MIEC anode model will be presented including the physical effects mentioned above.After a literature overview of MIEC anode models in Section 2.2, the model will at first be described in a general way, applicable for many MIEC materials in Section 2.3.In a second step, it will be specified and parametrized for Ce 0.9 Gd 0.1 O 1.95Àd in Section 2.4.

Literature overview of MIEC anode models
For dense thin-film electrodes, a number of sophisticated models were reported.Ciucci et al. 12 presented a 2D stationary simulation study with patterned metal current collectors to study the transport pathways of the charge carriers for symmetric thin-film electrodes.They concluded that the surface reaction rather than the migration of the charge carriers is the overall rate-limiting step for this particular cell setup.Chen et al. 13 performed a linear perturbation study in the frequency domain to model the impedance response of samarium doped ceria (SDC) thin-film electrodes.They showed that depending on the spacing of the contacts on the SDC thin-film, either the surface reaction resistance or the transport resistance of the charge carriers become dominant.Based on these results they suggested the existence of a maximum electronic penetration depth.
Jamnik and Maier 16 developed a general transmission line model for thin-film MIEC electrodes and discussed the decisive materials parameters.The latter determine whether or not the impedance spectrum is of Warburg type i.e. dominated by transport limitations of the charge carriers.Lai and Haile 11 presented a rigorous derivation of the impedance behaviour of mixed conducting materials using transmission lines and performed a case study for ceria.Nenning et al. 2 presented a two-mode measurement setup for thin-film MIEC electrodes, which enables the determination of the ionic and electronic conductivities as well as the electrochemical reaction resistance by fitting the EIS-spectra to corresponding transmission line models.Their approach was tested for SrTi 0.7 Fe 0.3 O 3Àd .
For porous MIEC high performance anodes, only few studies are available.Kishimoto et al. 14 developed a one-dimensional numerical DC model of a nickel-infiltrated GDC anode.He especially studied the microstructure-performance relationship, where the microstructure information was obtained from focused ion beam tomography.However, a detailed analysis of the transport mechanisms of charge carriers and their interdependency with the surface reaction and the influence of the gas impedance was out of their scope.

2.3
Model formulation for a general MIEC anode on the button cell level 2.3.1 Summary of the model features.This section is intended as a guide for the following MIEC anode model formulation but also as a comprehensive summary to prepare for the results and discussion sections.The model intends to represent the conditions of a button cell setup as illustrated in Fig. 1.
1D simulation setup.For this study, a 1D computational model is developed and implemented in the commercial software package COMSOL Multiphysics. 25The model consists of two computational domains as visualized schematically in Fig. 2. The domain x = [0,L] is the porous MIEC anode functional layer, where the hydrogen oxidation reaction (HOR), the transport of the charge carriers in the MIEC and the transport of the gas species in the porous media take place.As the pores and the MIEC phase are not spatially resolved, the effects of the porous microstructure are taken into account by using effective transport properties.They can be related to bulk material properties with a so called microstructure factor M for each phase as described in Section 2.3.2.The M-factor characterizes the considered microstructure by accounting for its phase volume fraction, tortuosity and constrictivity.The domain x = [ÀL Stag ,0] represents the stagnant gas layer at the interface between the gas flow channel and the anode functional layer.There, gas species diffusion corresponding to a button cell setup as depicted in Fig. 1 is taken into account.
Hydrogen oxidation reaction at the MIEC surface.In this contribution, humidified hydrogen is used as fuel.Due to the HOR, hydrogen is consumed and water, electrons and oxygen ion vacancies are produced at the MIEC/pore interface.The reaction kinetics depends on the actual material system and is formulated for CGO in Section 2.4.4.For CGO, the reaction kinetics especially depends on the concentration of electrons c eon and vacancies c vac .The surface reaction resistance in combination with chemical capacitance provoke the impedance contribution Z SR in the EIS-spectra.Note also that the production of an oxygen ion vacancy is equivalent to the consumption of an oxygen ion.Accordingly, the transport of an oxygen ion from the electrolyte to the MIEC/pore interface is equivalent to the transport of an oxygen ion vacancy in the opposite direction.In our model, we consider vacancies instead of the oxygen ions.
Charge transport.Electrons and vacancies are transported from the reaction site at the MIEC/pore interface throughout the MIEC to the CCL and to the electrolyte, respectively.The transport of the charge carriers is governed by drift and diffusion and is implemented based on the full Nernst Planck Poisson equation as described in Section 2.3.3.The charge transport losses in the MIEC provoke the impedance contribution Z chrgtpt in the EIS-spectra.
Gas transport.The gas species are transported from and to the reaction sites by convection and diffusion.The concentration changes along the gas supply cause a gas concentration impedance based on the Nernst equation (eqn (23)).Concentration gradients along the channel are caused by the fuel utilization (gas conversion impedance).This effect can be suppressed by an excess fuel supply, which is common practice for button cell experiments, as depicted in Fig. 1.Therefore, this effect is neglected in the current model.The gas transport in the gas channel is governed by convection.However, there exists a stagnant gas layer above the electrode as illustrated in Fig. 1, where the gas transport is governed by diffusion and not by convection of the excess fuel supply.In the domain x = [ÀL Stag ,0], the gas species diffusion in a stagnant gas layer is thus modeled.The thickness of this stagnant gas layer L Stag depends on the actual button cell setup as discussed in further detail in Section 2.3.6.A Stefan-Maxwell formulation is used for the gas transport in the stagnant gas layer as described in Section 2.3.8.
The gas transport within the pores of the MIEC-electrode is governed by diffusion.However, for fine porous microstructures with pores in the range of sub-microns, the transport can no longer solely be described by ordinary bulk diffusion but is additionally affected by Knudsen diffusion.The latter describes the scattering events of the gas molecules with the pore walls.Fine porous microstructures are relevant because of their potentially high electrochemically activity associated with their high specific surface area.To capture the Knudsen effects, the well known dusty-gas model (DGM) is used.It has proven able to describe the combined transport mechanism by bulk and Knudsen diffusion appropriately, see Section 2.3.7 and Sections B, C of the ESI.† The concentration gradients of the gas species in the stagnant gas layer and in the pores of the MIEC-electrode cause a gas diffusion impedance Z gas , as further discussed in Section 2.3.6.
By considering the gas transport and source terms from the reaction, the gas species concentrations are calculated.Moreover the oxygen partial pressure is calculated according to the dissociation of water.Therewith, oxygen partial pressure dependent material laws for the concentration and conductivity of electrons and vacancies and for the reaction kinetics can be used.
Model limitations.The described model entails some limitations.The electrolyte at x = L and the current collector layer (CCL) at x = 0 are modeled as boundary conditions and hence are not spatially resolved.The transport of electrons across the MIEC/CCL interface and the transport of vacancies across the MIEC/electrolyte interface are assumed to be lossless.In addition, the electrolyte and the CCL are modeled as ideal conductors.This means that the resistances of the CCL and the electrolyte are neglected and therefore do not induce any inhomogeneous in-plane current densities.In reality, losses from unfavourable in-plane conduction to the ribs of the interconnector might lead to electrochemically less active regions.According to the assumptions necessary for the 1D setup, the in-plane transport in y-and z-directions for the poreand the MIEC-phase of the porous electrode is assumed to be ideal.However, this simplification is quite appropriate.On the one hand, there is no macroscopic in-plane charge carrier transport because of the assumption of an ideal current collector and an ideal electrolyte.In the same manner for the porephase, the assumption of an excess fuel supply leads to the absence of a concentration gradient in y-and z-directions.There are thus no macroscopic in-plane diffusion fluxes.On the other hand, for the pore-scale level of the MIEC-phase, the in-plane transport of the charge carriers happens only from the surface to the center of a sintered MIEC-particle.These distances are on the sub-micrometer scale.In a similar way, no steep in-plane concentration gradients are expected within the single pore bodies.Moreover, these in-plane pathways on porescale are not tortuous and no constrictions due to bottlenecks are present within one single grain or one local pore body (bulge).As a result, the effective transport properties are assumed to be very close to the intrinsic ones.Therefore, the in-plane transport resistance can be regarded as negligible for gas-transport in the pores as well as for charge-transport in the MIEC-phase.However, for transports in the through-plane direction, limitations from microstructure due to tortuosity and constrictivity (bottlenecks) become effective and must be considered in the model.
2.3.2Effective transport properties from microstructure analysis.As already mentioned, the pores and the MIECphase are not spatially resolved.Therefore, the effects of the porous microstructure is incorporated using effective transport properties.They can be related to the bulk material properties with a microstructure factor M. 26 The relations for e.g. the effective diffusivity then reads: The M-factor can be determined based on 3D voxel data in two ways: (a) either using direct numerical simulation on the voxel grid, or (b) via morphological analysis.In (a) the M-factor is interpreted as a relative property (e.g.relative diffusivity M = D eff /D bulk ), which is indirectly derived from the knowledge of effective and intrinsic properties.In this work, the effective properties are determined by means of numerical simulations on the 3D geometry of the microstructure represented by a voxel grid using GeoDict 27 software.The diffusion equation in the solid MIEC phase is solved to determine the effective diffusivity and therewith the M-factors for the diffusion of the charge carriers.The derived M-factor can be used for the drift and diffusion of the charge carriers, as the equations for transport by diffusion and drift (i.e.electric conduction) are mathematically equivalent, i.e.M sim MIEC = D eff /D 0 = s eff /s 0 .The diffusion equation in the pore space is solved to determine the effective diffusivity and therewith the M-factors M sim pore = D eff /D 0 for the bulk diffusion of the gas species.Additionally, the gas flow permeability k flow for the parametrization of the Darcy law is determined by a flow simulation on the voxel grid.Moreover, an M-factor for the Knudsen diffusion is calculated.These M-factors for the pores space are needed for the parametrization of the dusty-gas model presented in Section 2.3.7.The used GeoDict solvers and their parametrization are documented in the results Section 3.6.
In (b) the M-factor is interpreted as a correction factor, which describes the limiting effects of microstructure and is directly obtained from microstructure characteristics.The M-factors for gas diffusion, as well as for conduction and diffusion of the charge carriers are additionally determined via such morphological analysis of the 3D microstructure data.Therewith, additional insight about the contribution of different microstructure effects on the transport properties can be achieved.Such methodologies have been developed and applied in previous publications by Holzer et al., [28][29][30] Gaiselmann et al., 31 Pecho et al., 32 Neumann et al. 33,34 and Stenzel et al. 35 Thereby, the M-factor includes the following microstructure effects: the phase fraction (e.g.porosity) f ph , the tortuosity t and the constrictivity b, relevant for drift and diffusion.A virtual materials testing (VMT) approach was applied in order to determine the M-factor empirically. 26,31,35hereby a large number (48000) of virtual 3D microstructures were created with a stochastic model.These structures cover a large range of values for f ph , b and t.For each 3D structure, the effective diffusivity D eff was determined by numerical simulation with GeoDict 27 software.The following expression was then found by statistical error minimization: This equation for M pred used in (b) was validated for different porous materials by comparing it with either results from simulation (M sim used in a) or from experiments (M exp ).Gaiselmann 31 and Stenzel 35 showed that the equation M pred gives good results for all three phases in SOFC cermet electrodes.Furthermore, validations for this equation were also performed successfully by Holzer et al. for very different types of porous microstructures, such as sintered ceramic membranes, 30 fibrous GDL in PEM fuel cell 36 and even open cellular materials (Holzer, unpublished data).These successful validations confirm that the equation for M pred has a global meaning, in the sense that it is capable to predict effective transport properties for all kinds of microstructures based on only three parameters (f ph , b, t).
For the solid MIEC-phase, the M-factor for the drift (i.e.electric conduction) and diffusion of the charge carriers within the MIEC yields: where f MIEC is the solid volume fraction, b MIEC the constrictivity and t MIEC the tortuosity of the MIEC-phase.
Likewise, for the bulk-diffusion of the gas species in the pore phase of the electrode the M-factor yields: e pore 1:15 b pore 0:37 where e pore is the porosity of the MIEC-phase and b pore the constrictivity and t pore the tortuosity of the pore-phase.
Note that morphological analysis are also available for the gas flow permeability as e.g. reported by Holzer et al. 30

and
This journal is © the Owner Societies 2021 Neumann et al. 33 As the results of the dusty-gas model (see Section 2.3.7) are less sensitive to the permeability as to the diffusivity, we only report the permeability extracted from the flow simulation on the voxel mesh.

Charge transport by drift and diffusion.
To model the transport of electrons (eon) and vacancies (vac), the corresponding drift diffusion equations 11 read: J eon = ÀD eon,eff rc eon À z eon m eon,eff Fc eon rf.
Here, J vac and J eon are area specific the molar fluxes, c vac and c eon the molar concentrations, f the electric potential, D vac,eff and D eon,eff the effective diffusion coefficients, z vac and z eon the number of elementary charges and m vac,eff and m eon,eff the effective mobilities of the vacancies and electrons, respectively.Furthermore, F denotes the Faraday constant.In addition, the Nernst-Einstein relation is used for the relation between the diffusion coefficient and the mobility: 37 where T is the temperature, R gas the ideal gas constant and i denotes vac or eon.The effective diffusion coefficient D i,eff is determined according to eqn (1): 2.3.4Continuity equations for the charge carriers.The continuity equations for electrons and vacancies reads: The reaction rates for the vacancies R vac and electrons R eon are given by: where i SR is the surface reaction rate at the MIEC/pore interface.Note that i SR is expressed as a source term because of the used 1D setup.The assumption for the HOR occurring at the MIEC/ pore interface has been shown experimentally for CGO-based anodes by several authors. 11,20,22,38ollowing Kishimoto et al. 14 i SR is expressed as a Butler-Volmer like expression with a surface reaction overpotential Z SR described below: where i 0 is the exchange current density.Note that the Butler-Volmer formulation is also applied to traditional Ni-YSZ anodes, where the charge transfer process is assumed to occur at the three-phase boundaries.Thereby, the negatively charged oxygen ions are transported through the YSZ to the three-phase boundaries where they are consumed and two electrons are transferred over the phase boundary to the Ni-phase.This charge transfer can be associated with a welldefined overpotential between the YSZ-and the Ni-phase. 39owever, MIEC-anodes do not exhibit a classical charge transfer process over a phase boundary like in Ni-YSZ anodes.Hence, oxygen ions are transported trough the MIEC to the MIEC/pore interface, where they are consumed.Electrons are produced at the MIEC/pore interface and have to be transported away from it.Oxygen ions (or vacancies respectively) and electrons are transported in the same MIEC-phase.However, the use of an overpotential formulation according to the Butler-Volmer equation is justified in the sense that an electrical double layer is formed at the MIEC/pore surface due to the charged surface reaction species.But the relation between the surface potential and the measurable overpotential of the surface reaction (or activation overpotential) Z SR can be complex, as derived by Fleig 40 for a MIEC cathode.Moreover, the surface reaction overpotential cannot be directly associated to an electric potential difference within the same MIEC-phase.We therefore choose the approach that the surface reaction is driven by the chemical potential drop, which is associated with the surface reaction as suggested by Kishimoto et al. 14 However, this surface reaction overpotential Z SR needs to be formulated appropriately in order to reflect the rate determining reaction step for the particular material system.An appropriate formulation for CGO is suggested in Section 2.4.4.
Note that the overall driving force for the HOR is the chemical potential difference between the reactants H 2 þ 1=2O 2 and the product H 2 O.However, as we use a prescribed current density in our model formulation, we will directly account for the different overpotentials needed to drive this prescribed current.The available voltage of the fuel cell is then the chemical potential difference between reactants and product minus the sum of all overpotentials.Nonetheless, we will restrict our discussion on the overpotentials and ASR, because the discussion of current-voltage curves is out of the scope.
Moreover, the Poisson equation needs to be fulfilled, i.e.

rÁ(Àe
where r v is the space charge density and f the electric potential.The space charge density is the sum of all charges: 2.3.5 Oxygen partial pressure.Many MIEC materials undergo a change in stoichiometry in the reducing anode atmosphere.The nonstoichiometry for these materials is thus a function of the oxygen partial pressure p O 2 .The nonstoichiometry can have an important impact on material properties like the conductivity of electrons and vacancies and the reaction kinetics, as will be discussed in further detail for the case of CGO in Section 2.4.In order to use these p O 2 dependent material laws, the local p O 2 is estimated from the gas-phase equilibrium of the local atmospheric conditions (see ESI, † Section A).
2.3.6 Model for the gas concentration impedance.Depending on the cell architecture, the gas concentration impedance can be a significant or even the dominant part of the impedance spectrum.In this section the approach and assumptions to model the gas concentration impedance of a button cell are elaborated.The general term gas concentration impedance can be used for the change of the Nernst potential due to concentration changes (nomenclature from Bessler 41 ).Concentration gradients can be a result of limited convective transport of the fuel along or by diffusion limited transport perpendicular to the electrode as illustrated in Fig. 3.
Convective transport along the cell is limited in technically relevant systems in order to achieve a high fuel utilisation.This results in a concentration gradient along the cell (see Fig. 3) due to the consumption of hydrogen and the production of water vapour.This so-called gas conversion impedance 41 is usually the dominant contribution to the gas concentration impedance for technically relevant large area cells and stacks.However, the gas conversion impedance can be suppressed in button cell configurations by using an excess fuel supply and is therefore neglected in our model.
Concentration gradients perpendicular to the electrode are because of diffusion limitations from different contributions and can be summarized by the term gas diffusion impedance. 41 possible diffusion limitation within the gas channel is associated with a stagnant gas layer, where the convective fuel flow does not contribute to mass transfer because of the no slip condition at the electrode surface (see Fig. 3).In contrast to the gas conversion impedance, there is no simple measure to avoid this effect.Nevertheless, Primdahl and Liu 4 presented a setup with a MIEC-electrode mounted between auxiliary electrodes to avoid significant impedance contributions from gas diffusion outside the porous electrode structure.However, as our model intends to describe a classical button cell setup with an axial fuel supply as depicted in Fig. 1, a separate computation domain is introduced in the 1D-model for the stagnant gas layer.There, the gas diffusion is modeled with a Stefan-Maxwell approach, described in Section 2.3.8.At the boundary to the non-resolved gas channel region beyond the stagnant gas layer, the gas species concentrations are fixed, which mimics excess fuel supply without any concentration gradients along the button cell (i.e. the gas conversion impedance is suppressed).
Within the porous electrode, the gas transport is also dominated by diffusion and hence contributes to the gas diffusion impedance, which is illustrated in Fig. 2. In this simulation domain, the gas transport is modeled after the dusty gas model (DGM) which accounts for bulk and Knudsen diffusion, see Section 2.3.7 and Section B of the ESI.† The determination of the appropriate diffusion coefficients for the models is described in Section C of the ESI.† Note that in button cells a metal mesh is often placed on top of the anode to ensure a good electric contact.This additional diffusion resistance is neglected in our model.Likewise, also the diffusion resistance in the not spatially resolved ideal CCL is neglected.
2.3.7 Dusty-gas model (DGM) for the gas transport in the MIEC-pores.Within the porous MIEC anode, the transport and the reaction of the gas species needs to be described.The continuity equations for the gas species reads: where c 1 , c 2 are the molar concentrations, N 1 , N 2 are the molar flux densities and R 1 , R 2 are the reaction rates for hydrogen (species 1) and water vapour (species 2), respectively.Thereby, the consumption rate of hydrogen is equal to the production rate of water vapour, which is equal to the production rate of oxygen ion vacancies: where R vac is defined in eqn (11).The corresponding boundary conditions for the continuity eqn ( 16) and ( 17) are documented in Section 2.3.10.
For the molar flux densities N 1 and N 2 , the effect of the porous microstructure needs to be considered.In nanoporous materials, gas diffusion is controlled by collision of gas molecules with the pore walls.In this case the mean free path of the gas molecules is longer than the characteristic length of the pores so that we enter the Knudsen diffusion regime.Modern high performance MIEC electrodes are usually very fine-grained in order to increase the specific surface area and associated reactivity.In most cases, the pore size in modern MIEC electrodes is in the mm to sub-mm scale, where both Knudsen and bulk diffusion become relevant.It is generally agreed 14,42,43 Fig. 3 Illustration of gas concentration impedance with contributions from gas conversion and gas diffusion.
that the dusty-gas model (DGM) is the most convenient approach for modeling combined bulk and Knudsen diffusion.Therefore, in our approach, the diffusion in the porous CGOlayer is modeled with the DGM.The DGM is an extension of the Stefan Maxwell model, where the pore phase (= ''dust'') is treated as an additional, stationary species. 42In this work, an explicit version of the DGM for a binary system derived by Liu et al. 44 is used.A detailed presentation and discussion of the used equations can be found in Section B of the ESI.† The DGM needs to be parametrized appropriately.The bulkdiffusion contribution is governed by the effective binary diffusion coefficient D 12,eff : where D 12,0 is the intrinsic binary diffusion coefficient between hydrogen and water determined according to Fuller et al. 45 and the microstructure factor M pore for the pore phase defined in eqn (4) accounts for the effect of the microstructure.The effect of the Knudsen diffusion is governed by the effective Knudsen diffusion coefficients for hydrogen D 1Kn,eff and water D 2Kn,eff , respectively: where D iKn,0 is the intrinsic Knudsen diffusion coefficient according to the kinetic theory of gases, 46 M pore,Kn is the microstructure factor for the Knudsen diffusion mechanism and i denotes species 1 (hydrogen) or 2 (water vapour).It is worth noting that the intrinsic (and thus also the effective) Knudsen diffusion coefficient is proportional to a characteristic pore diameter d pore of the porous structure, while the effective binary diffusion coefficient D 12,eff is size independent and does only depend on the morphology of the microstructure.The relative contribution of the bulk and Knudsen diffusion on the overall diffusion resistance can be estimated by the Knudsen number, which is the ratio of the bulk and Knudsen diffusion coefficients.However, the bulk microstructure factor M pore and the Knudsen microstructure factor M pore,Kn might be similar, but not necessarily equal, as different processes are hindering the transport for the two transport mechanisms.Therefore, we distinguish two kinds of Knudsen numbers in this work: the intrinsic Knudsen number Kn i,0 comparing the intrinsic properties and the effective Knudsen number Kn i,eff comparing the effective properties.The definitions of both versions of the Knudsen number are presented and discussed in further detail in Section D of the ESI.† Finally, the DGM includes a convective part, which is governed by the gas-flow permeability k flow of the microstructure and the dynamic viscosity of the gas mixture m visc .
All the mentioned parameters needed for the DGM are tabulated in Tables 1 and 2. A detailed presentation of the equations and discussion of the binary and Knudsen diffusion coefficients can be found in Section C of the ESI.† 2.3.8Stefan-Maxwell formulation for the gas transport in the stagnant gas layer.There is a stagnant gas layer above the electrode, where the gas transport is governed by diffusion and not by convection of the excess fuel supply.The thickness of this stagnant gas layer strongly depends on the flowrate and on the experimental setup of the button cell.For our standard simulations, a quite small stagnant gas layer thickness of L Stag = 0.2 mm is used.This value was estimated from a flow simulation of a particular stagnation point-flow like setup (see ESI, † Section F).
The continuity equation for the stagnant gas layer reads: where N 1 , N 2 are the molar flux densities for hydrogen and water vapour, respectively.Boundary conditions are documented in Section 2.3.10.The molar fluxes densities N 1 and N 2 in the stagnant gas layer are modeled with the Stefan-Maxwell 42 equations as documented in detail in Section E of the ESI.† Note that for the used simplified version of the Stefan-Maxwell model, the only governing parameter is the intrinsic binary diffusion coefficient as defined in eqn (23) in Section C of the ESI.† 2.3.9Gas diffusion overpotential.The gas concentration overpotential is the change of the Nernst potential due to concentration changes of the gas species. 41According to the model assumptions stated in Section 2.3.6,only diffusive contributions are considered.Therefore, the term gas diffusion overpotential is used, which describes the effect of the gas diffusion resistance in the stagnant gas layer and in the porous MIEC-electrode.Hence, the gas diffusion overpotential is the basis to calculate the gas diffusion impedance in eqn (38).
The gas diffusion overpotential 41 Z gas can be calculated as follows: where x 1,0 and x 2,0 are the mole fractions of the fuel supply at the stagnant gas layer on the channel side (i.e. at x = ÀL Stag ).
The gas diffusion overpotential results from averaging over the species concentration present in the electrochemically active MIEC-layer.

Boundary conditions
Boundary conditions for the drift-diffusion equation.At the MIEC/CCL interface at x = 0, the vacancies are blocked and there is a prescribed flux for electrons: where J charge is the prescribed current density.At the MIEC/ electrolyte interface at x = L, the electrons are blocked and there is a prescribed flux for vacancies: Boundary conditions for the Poisson equation.The Poisson eqn ( 14) is coupled with the drift diffusion eqn ( 5) and ( 6) over the space charge density eqn (15).The voltage drop over the electrode therewith results from the drift transport losses of the charge carriers.Note that the total fluxes of the charge carriers are prescribed at the boundaries in eqn ( 24)- (27).
Per definition, we set the MIEC potential at the interface to the CCL to ground:  The electric potential at the MIEC/electrolyte interface results from the full voltage drop over the electrode: where E = Àrf is the electric field.
Boundary conditions for gas transport.The goal of the simplified model for the gas concentration impedance described in Section 2.3.6 is to mimic an ideal fuel supply, but respecting the diffusion limitations in the stagnant gas layer.Therefore, the concentrations at the end of the stagnant gas layer at x = ÀL Stag are fixed to the given gas concentrations of the fuel supply of hydrogen c 1,0 and water vapour c 2,0 : A zero-flux boundary condition is implemented for all gas species at the CGO/electrolyte-interface, as the non-porous electrolyte is assumed to be impermeable for the gas species: 2.3.11Frequency domain perturbation for the simulation of the EIS-spectra (AC-mode) and implementation in Comsol Multiphysics.In order to simulate the EIS-spectra, a small sinusoidal perturbation with angular frequency o = 2pf is superimposed on a steady current density and the response of the corresponding overpotentials is calculated.This is exactly the same approach as used for experimental EIS measurements.In the time-domain, the perturbed current density reads: where J ˆcharge is a large steady current density, J ˜charge is a small perturbation current density and J charge (t) is the total current density.This perturbation provokes a time dependent response of the dependent variables fðtÞ ¼ b f þ e f expð jotÞ and c i (t) = c ˆi + c ˜i exp( jot), where the variables with the b are the stationary solutions, the variables with the e are the small perturbed components and i denotes eon, vac, 1 or 2, respectively.
However, in our simulation approach a linear perturbation in the frequency domain is used in order to compute the EISspectra.Therefore, the equations are linearized around the stationary operating point c ˆi, b f.By substituting fðtÞ ¼ b f þ e f expð jotÞ and c i (t) = c ˆi + c ˜i exp( jot) in eqn ( 9), ( 10), ( 16), ( 17), ( 21) and (22) and calculating the first order Taylor expansions, a time-independent linear system for the perturbed components c ˜i, f is obtained, which is solved for each frequency. 47For the linearized system, the perturbed components c ˜i, f are the linear deviations from the linear perturbation.The boundary conditions for the perturbed system are obtained analogously.The linear perturbation of the current density is applied at the boundary conditions of eqn (25) and (26).The calculations are performed using the ''Frequency Domain Perturbation'' of the electrochemistry module of Comsol Multiphysics. 48For the DGM-model, the general form implementation of Comsol Multiphyiscs was used and the frequency domain perturbation terms were added as source terms, as the they are not included in the standard implementation of the general form.
From the linear deviation results, the different impedance contributions can be computed.The impedance according to the charge transport is computed from the linear deviation of the electric potential at the interface to the electrolyte f(x = L) as the electric potential is set to zero at the CCL: The surface reaction resistance/chemical capacitance impedance is calculated from the average surface reaction overpotential: The average surface reaction overpotential ZSR,avg is calculated from the local surface reaction overpotential ZSR according to eqn (37) (derivation is documented in the ESI, † Section G): Note that the local surface reaction overpotential Z SR is a property specific to the used MIEC material.For CGO it is formulated in eqn (55) of Section 2.4.4.
The gas diffusion impedance is calculated from the average gas diffusion impedance overpotential from eqn (23): The total impedance of the anode is finally the sum of the three impedance contributions:

Application of the model for CGO10
In Section 2.3, a general anode SOFC model valid for MIEC materials in general has been proposed.However, the corresponding material laws for the nonstoichiometry, the numbers of charge carriers, the charge transport and reaction mechanisms are specific to the considered MIEC material and material composition.In order to achieve a consistent model description, in the following, the general MIEC model is specified and parametrized for gadolinium doped ceria with the composition Ce 0.9 Gd 0.1 O 1.95Àd , commonly abbreviated with CGO10.As has already been stated, the general model formulation is valid for many types of MIEC electrodes.If another MIEC electrode material is used, the formulation of this Section 2.4 needs to be adapted appropriately, which is further discussed in Section 2.4.7.
2.4.1 Crystal chemistry and determination of equilibrium densities.Gadolinium doped ceria owns the following general chemical formula: 49 Ce 1Ày Gd y O 2Ày/2Àd (40)   where y is the stoichiometry of gadolinium, i.e. the amount of dopant per amount of ceria, and d is the oxygennonstoichiometry, which depends on the oxygen partial pressure.
The crystal structure of CGO is schematically visualized in Fig. 4. Throughout this paper we use the concept of relative charges for the description of charge transport within the MIEC crystal.According to this concept it is not the actual charge of an ion which is counted.The charge of each species is always compared with the charge(s) of the species that is 'normally' occupying the considered lattice cite in a unit cell of an ideal crystal (i.e. in CeO 2 ).For example, on each oxygen site two positive charges are attributed to the oxygen vacancy relative to O 2À .In analogy, on the cation site, Ce 3+ and Gd 3+ have one negative charge relative to Ce 4+ .For more information about the relative charge concept see Kro ¨ger-Fink notation. 50,51n the following, we determine the equilibrium charge densities, i.e. the densities of charge carriers when no electric current is drawn.The crystal structure of CGO is schematically visualized in Fig. 4. As an intermediate step, the number density of cerium for pure ceria is expressed with the unit cell volume: where the factor 4 is due to the 4 cerium atoms in a unit cell, f MIEC is the solid volume-fraction of CGO in the porous MIEC electrode and the lattice parameter of the unit cell is L Lattice = 5.422 Å 52 for a gadolinium doping of y = 0.1.By gadolinium doping, one Ce 4+ ion is replaced by one Gd 3+ ion.
To compensate the negative charge of Gd 3+ relative to Ce 4+ , an oxygen vacancy is introduced for two doped Gd-atoms, see Fig. 4. The number density of the oxygen ion vacancies due to the Gd-doping is: These positively charged vacancies (relative to oxygen ions at this lattice site) are mobile and enable an ionic conductivity, while the negative charges of the Gd 3+ -ions relative to the Ce 4+ions are immobile.Overall, the crystal remains neutral.To express this electro-neutrality, we can formulate the number density of this negative immobile charge carriers as follows.
In addition to the charge compensation of Gd-doping, oxygen ion vacancies can be introduced due to reduction of the ceria in a reducing atmosphere (low oxygen partial pressure).The additional oxygen ion vacancy number density n vac,red due to this change in stoichiometry can be expressed as: where d is the oxygen nonstoichiometry.To compensate the positive charge of these additional vacancies, the reduction state of two cerium-ions per vacancy is changed from Ce 4+ to Ce 3+ .In equilibrium, the corresponding number density of electrons n eon,eq yields: Although these electrons are localized to Ce-ions, they can be transported due to small polaron hopping 53 from a Ce 3+ -ion to a Ce 4+ -ion as visualized in Fig. 4. For convenience, the molar densities c i of the charge carriers can be obtained from the corresponding number densities as follows: where N A is the Avogadro constant.For the equilibrium molar concentrations of the charge carriers we can therefore write: c eon,eq = 2c vac,red = 2c Ce d (48) Therewith, the space charge density as the source term for the Poisson equation can be specified by inserting c eon,imob = c eon,Gd,imob in eqn (15): 2.4.2Nonstoichiometry of CGO10.To determine the equilibrium concentrations of the charge carriers as described in Section 2.4.1, the oxygen-nonstoichiometry d of CGO needs to be known.Wang et al. 54 studied the oxygen-nonstoichiometry of CGO10 experimentally for different temperatures and oxygen partial pressures.The manually extracted data of Wang et al. 54 are plotted in Fig. 5 and the result of the used fit-functions are dependency for moderate oxygen partial pressures.However, for very low oxygen partial pressures as present in the anode of SOFC, there is a considerable deviation from this p O 2 À1/4dependency.For temperatures above 800 1C, the temperature dependency can be described with a translation of the d( p O 2 )curve in the log-log plot along the p O 2 -axis with reasonable accuracy.However, the curve for a temperature of 700 1C deviates significantly.As the intended operation temperature of the system we study is above 800 1C, we simply stick to the experimental data starting from 800 1C.For the numerical implementation, a simple fit-functions for the oxygennonstoichiometry is desirable.We therefore restrict ourself to a bi-linear expression for the p O 2 -dependency and a parabolic expression for the temperature-shift (see ESI, † Section H). 2.4.3Electronic and ionic conductivity of CGO.In this section, the mobilities and diffusivities of the charge carriers used in the drift-diffusion eqn ( 5) and ( 6) are determined for CGO10 according to available literature data for the ionic and electronic conductivities.
CGO is able to conduct double negatively charged oxygen ions for two reasons: due to oxygen vacancies introduced by the gadolinium doping and due to the oxygen-nonstoichiometry resulting from low oxygen partial pressures.The oxygennonstoichiometry also results in the presence of reduced Ce 3+ ions, which are negatively charged with respect to the original Ce 4+ ions.Via small polaron hopping, electrons can hop from Ce 3+ to Ce 4+ ions providing the MIEC's electronic conductivity (visualization in Fig. 4).
According to Steele 49 the ionic (vacancy) conductivity s vac and electronic conductivity s eon of CGO10 can be approximated by the following equations: where s 0 vac = 1.09Â 10 7 , s 0 eon = 3.456 Â 10 11 are the conductivity parameters and DH vac = 0.64 eV, DH eon = 2.475 eV the activation energies for vacancies and electrons, respectively.Note that eqn ( 51) is valid for temperatures T 4 400 1C.
The conductivity of the charge carriers is related to the mobility through: where i denotes either eon or vac.As the charge carriers are assumed to be transported by drift and diffusion according to eqn ( 5) and ( 6), the mobility and diffusivity of the charge carriers is required.Therefore, eqn ( 53) is solved for the mobility using a reference oxygen partial pressure of p This anode reaction incorporates several intermediate reaction steps as e.g.described by Feng et al. 55 Note that a detailed description of the intermediate reaction steps is out of scope of this study.Nevertheless, the aim is to formulate a simplified surface reaction overpotential accounting for the rate limiting step.Tabish et al. 56 summarized the different rate-determining steps for ceria electrodes reported in the literature and concluded that the majority of the studies depict a charge transfer step to be the rate-determining step, though there are differences in the understanding of how and where this rate limiting charge transfer process takes place.We therefore use the following formulation for the surface reaction overpotential Z SR , which implies that the charge transfer is the ratedetermining step and that all other intermediate steps are very close to the equilibrium.As a consequence, we assume that the species of H 2 , H 2 O and O Â O do not influence the surface reaction Fig. 6 Central orthoslices of the virtual microstructures.The corresponding microstructure parameters are listed in Table 2. overpotential: where c eon,eq and c vac,eq are the equilibrium (i.e. at open circuit) concentrations of the electrons and vacancies and c eon and c vac are the concentrations of the electrons and vacancies under load.According to this formulation, the electrons and vacancies are depleted to drive the reaction, which will be discussed in more detail in Section 3.2.If we insert Z SR from eqn (55) in eqn (13), the exponential and logarithmic terms cancel, yielding a rate expression with two branches: where i 0 is the exchange current density, which is further elaborated in the following Section 2.4.5.
2.4.5 Exchange current density.The exchange current density i 0 defines the rate of the HOR and hence strongly influences the surface reaction resistance.The following expression is suggested by Kishimoto et al.: where S is the specific surface area and k 0 the surface specific exchange reaction rate constant.According to various ref.-dependency as well.However, k 0 is hard to determine from EIS measurements in a reliable way, because of the influence of the microstructure and because of the overlapping processes of e.g. the gas concentration impedance.Therefore, we estimate k 0 from measurements of Gerstl et al. 57 for a thin-film model electrode based on CGO20.For the estimation of k 0 in our CGO anodes we assume that the surface reaction resistance is proportional to the CGO/pore interface area.This is described in detail in Section J of the ESI † and the estimated value for the exchange reaction rate k 0 is reported in Table 1.
Moreover, it is reported in the literature that the reaction kinetics are enhanced for CGO-based anodes with increasing humidity.Riegraf et al. 58  -dependency.However, the positive effect of the humidity seems to be significantly stronger.This behaviour is also observed for dense SDC thin-film electrodes in the work of Chueh et al., 59 where a power low of p H 2 O 0.3 was reported for the reaction kinetics with embedded platinum electrodes.The quantitative dependency of the surface reaction resistance on the water content thus needs to be determined experimentally for a specific material system, which is beyond the scope of the current study.Therefore, in this study we will hold the water friction constant for all simulations and just lump the effect of the water fraction in the rate constant k 0 .Moreover, positive effects of increasing water contents due to water vapour concentration gradients from gas diffusion processes are not considered.The exchange current density can be strongly enhanced with the addition of catalyst metals.Often, CGO-based anodes are realized as a composite together with Ni. 3,4,14,19 If the Ni-phase percolates, the electronic conductivity of composite is drastically enhanced, which cannot be described by the current model with one solid phase.Nonetheless, we suggest that Ni-CGO anodes with an Ni-phase far below the percolation threshold can be represented quite well with the current model by simply considering the effect of the Ni-phase by scaling the rate constant k 0 .
However, the role of Ni in Ni-CGO anodes is very complex and some selected literature references shall be briefly discussed here.Primdahl and Liu 4 reported a decrease in the impedance by about a factor of 4 when adding a small amount of Ni of e.g.2.5 wt% to a CGO40 anode.This drastic decrease of the reaction resistance cannot be explained by the variation of the microstructure and therefore was associated with the catalytic activity of Ni.However, in the work of Kleinlogel et al., 60 the positive effect of Ni was interpreted to origin from an improved microstructure, and the impact of Ni was explained as a sintering aid.A further interpretation of the role of Ni was presented by Primdahl and Mogensen, 3 who reported for a Ni-CGO10 anode with a volume ratio of 60/40 that the impedance decreased about a factor of three after adding 3 wt% of additional fine dispersed Ni (i.e.Ni addition via impregnation of the porous, sintered composite).This implies that fine dispersed Ni particles on the surface of CGO show a higher catalytic activity than the ''bulk'' Ni phase, contacting the CGO only at the phase boundaries.However, the authors also suggested that this improvement might be because of an increase of the three-phase boundaries or possibly by an improved sintering of the contacts between Ni and CGO10.
For our model we assume that the non-percolating Ni-phase in the Ni-CGO anode has only a catalytic effect, which is considered by scaling the rate constant k 0 .However, it cannot be excluded that the addition of Ni might also add a small additional intermediate-frequency process.Such a process was e.g. reported by Riegraf et al. 19 for Ni-CGO anodes in the range of 20-100 Hz.In the literature, this process is either interpreted as a charge transfer process 3,61 or the oxygen ion transport across the MIEC/electrolyte interface. 3The latter effect is probably independent of the Ni-phase.However, both possible effects are neglected in our model.
2.4.6Chemical capacitance of CGO-electrodes.It is well known that ceria-based electrodes involve a chemical capacitance due to a change in the oxygen-nonstoichiometry. 11,13,16,22,23 The chemical capacitance can be generally expressed as follows: 62 where n i is the number density of mobile charge carriers, e the electron charge, A is the area and L is the thickness of the sample.For CGO, we have two charge carriers: This journal is © the Owner Societies 2021 Oxygen ion vacancies, abbreviated as vac with an elementary charge of z vac = +2.The corresponding number density can be expressed as n vac,eq = c vac,eq N A , where c vac,eq is the molar concentration of oxygen ion vacancies defined in eqn (47).Ce 3+ -Ions, abbreviated as eon with an elementary charge of z eon,eq = À1.The corresponding number density can be expressed as n eon,eq = c eon,eq N A , where c eon,eq is the molar concentration of Ce 3+ -ions (electrons) defined in eqn (48).Therewith, the chemical capacitance for the current CGO material system is given by: where the identity R gas = k B N A was used.However, the effect of this chemical capacitance on the impedance spectra is already included in the governing equations and does not need to be implemented separately.Note that eqn ( 59) is also valid with other dopants as e.g.samarium or even for undoped ceria, if the corresponding charge carrier concentrations are used.
In addition to the chemical capacitance, there is also a surface capacitance present at the CGO/pore interface.Chen et al. 13 included this surface capacitance in their model for SDC thin-film electrodes.They estimated the value of the surface capacitance on the basis of the experimental work about the capacitance in SDC thin-films of Chueh and Haile. 17They associated the capacitance independent of the electrode thickness to be the surface capacitance at the SDC/gas interface, assuming that a possible interface capacitance at the SDC/ metal interface is negligible.However, as there is no specific data available for the surface capacitance of porous CGO, we do not consider this effect in the present study.
2.4.7 Application of the model for other MIEC materials.The general MIEC-model formulated in Section 2.3 is specified for CGO10 in this Section 2.4.If another MIEC material than CGO10 shall be used, the following parameters and equations need to be modified: The charge carrier concentrations as a function of the oxygen partial pressure (eqn ( 47)-( 49)).Ionic and electronic conductivities as a function of p O 2 (eqn ( 51) and ( 52)) to determine the mobility of the charge carriers as described in Section I of the ESI.† Formulation of the material specific surface reaction overpotential (eqn ( 55)).Dependency of the exchange current density on the p O 2 and water content (eqn ( 57)).The model can be easily applied for other ceria-based anodes like CGO20, CGO40 or Sm-doped ceria, as reaction and charge transport mechanisms are similar.For composite anodes where both phases contribute significantly to the charge transport (as e.g. for Ni-CGO anodes where both phases percolate), an additional transport-phase is needed in principle.However, effective conductivities of the composite structure might be used as a first approximation.
For other MIEC materials as e.g.titanate-based perovskite anodes like LST, 63 LSCT, 5 STFO 2 , etc., the material laws need to be adapted more carefully as additional physicochemical processes are expected.Burnat et al. 63 e.g.presented an additional process in the EIS-spectra of LST-CGO composite anodes, which was associated with an additional reaction pathway at the three phase boundaries.

Simulation results
The focus of this paper is the description of the mathematical model and the discussion of the complex physico-chemical processes in CGO-based anodes.We therefore do not compare our model results to experimental EIS-spectra of distinct cell setups, nor do we discuss in detail the dependency on the operating conditions.Both points are discussed conceptually in a previous publication (Marmet et al. 15 ).However, realistic model parameters are used to match the order of magnitude for resistive/capacitive effects of the processes.

Parameter values
Note: the parameters in Table 1 are the standard parameters used for the model.Most of these parameters are estimated based on literature data as referred in the reference-column.The anode layer thickness L and the operating conditions can vary in a wide range and are therefore chosen according to illustrative aspects.The relatively large anode layer thickness of L = 100 mm is chosen in order to illustrate the effects of gas and charge carrier transport.Note that the anode functional layer thickness is typically thinner (e.g.L o 30 mm 19 ) for the commonly used Ni-CGO cermets, where the reaction resistance is lower due to the catalytic activity of Ni and the ionic conductivity is lower due to the lower CGO volume fraction.The current density is chosen according to an internal standard value (a current of 300 mA for an active button cell surface of 1.44 cm 2 ).Parameters which are changed for parameter studies are mentioned separately in the subsequent text.The effective properties of the virtual microstructures MS1, MS2, and MS3 are listed in Table 2.All the listed properties are determined from the virtual microstructures according to microstructure analysis and simulations described in Section 2.3.2.The microstructure MS 2 is used for all the simulations where not otherwise specified.

Steady state solution for an initial parameter set
Steady state results for the parameter set of Table 1 are obtained by solving the eqn ( 9), ( 10), ( 16), ( 17), ( 21) and ( 22) without the time dependent terms.In Fig. 7(a), the concentration of the double relative positive charged oxygen ion vacancies (hereafter simply called as vacancies) and the concentration of the single relative negative charged Ce 3+ -ions (hereafter simply called electrons) are plotted as a profile across a 100 mm thick anode layer.Note that there are less than twice as many mobile electrons as vacancies, because of the immobile negative charge of the Gd 3+ -ions.These charge carrier concentrations resulting for a prescribed current load deviate from the equilibrium at open circuit conditions.The deviation from the equilibrium, i.e. the difference of open circuit to an effectively drawn current, defined as Dc i = c i À c i,eq (where c i,eq is the equilibrium concentration which is established at OCV) is plotted in Fig. 8(a).We can observe a negative deviation over the whole CGO-layer.According to the eqn (54) for the overall anode reaction, the products (i.e.electrons and vacancies) need to be depleted in order to drive the reaction.Moreover, it isdespite the very simplified reaction mechanism used in this study -in qualitative agreement with the in operando study of Feng et al., 55 where the response of surface oxygen vacancies and electrons were tracked using synchrotron-based ambient pressure X-ray photoelectron spectroscopy (APXPS).With these APXPS measurements, Feng et al. documented a depletion of oxygen vacancies and Ce 3+ -ions at the surface of the MIEC material under current load.Furthermore, Fig. 8(a) also illustrates that the deviation from equilibrium increases towards the CGO/CCL and CGO/electrolyte interfaces.This is in accordance with the distribution of the reaction along the CGO-layer as shown in Fig. 8(b), which is enhanced towards the CGO/CCL and the CGO/electrolyte interfaces.The distribution of the reaction is illustrated with the electron and vacancy source terms.
As mentioned in Section 2.3.4,the overall driving force for the reaction is the chemical potential difference between the reactants H 2 þ 1=2O 2 and the product H 2 O.This overall driving force then provokes the local driving forces for the different processes.In order to discuss the different driving forces for the charge carrier transport, the flux components due to drift and diffusion are plotted in Fig. 7(c) for the electrons and Fig. 7(d) for the vacancies.Note that only the total fluxes are physical and the separated drift and diffusion fluxes are virtual fluxes in order to illustrate the driving forces.The driving forces of the diffusive fluxes are the concentration gradients.A positive flux means that the charges are transported towards the electrolyte, which is characteristic for vacancies, whereas a negative flux (towards the current collector) is the overall direction for electrons.However, the total flux of each charge carrier is composed of drift and diffusion components.In this context it is important to note that the direction of diffusive flux components may be opposite to the overall transport direction for the same species, depending on the local concentration gradients.As shown in Fig. 7(a) the concentration profiles exhibit a maximum, so that the diffusive fluxes tend to go in the direction of the neighbouring interface (electrolyte or current collector, respectively).This is in contrast to the more constant potential gradient (Fig. 7(b)), which induces constant transport directions of the drift components (i.e.vacancies towards electrolyte, electrons towards current collector).This  This journal is © the Owner Societies 2021 leads to a puzzling picture for the different transport species and components as illustrated in Fig. 7(c) and (d).For example, in the region close to the electrolyte, the negative concentration gradients for both species induce diffusive fluxes towards the electrolyte.For the electrons, the diffusive flux component in the region close to the electrolyte is opposite to the direction of the drift component.Therefore, close to the electrolyte, the diffusive flux component hinders the transport of electrons to the CCL.As a result, the total electronic current in this region is close to zero, consistent with the blocking boundary condition for electrons at the electrolyte (eqn ( 27)).Vice versa towards the CCL, the concentration gradients favor the diffusive transport of electrons to the CCL and hinder the diffusive transport of vacancies to the electrolyte.Therefore, in the region close to the CCL, the total flux of vacancies is almost zero, consistent with the blocking boundary condition for vacancies at the CCL (eqn ( 24)).The total flux of electrons is twice as big as the flux of the vacancies, because the vacancies carry two unit charges per charge carrier.However, the mobility of the electrons exceeds the mobility of the vacancies by approximately a factor of 10, depending on the oxygen partial pressure.As electrons and vacancies experience the same electric field, the concentration gradients are adjusted according to the energy minimization of the system to achieve low charge transport losses.Therefore, the zone, where the vacancy transport is favored by the concentration gradients is larger than the zone, where the electron transport is favoured.This explains the asymmetric shape of the concentration profile in Fig. 7(a).Of course, the charge transport resistance provide not the only energetic contribution but adds to the surface reaction resistance and gas impedance, which are all captured in the anode model.
Our model thus provides a quantitative description of the distribution of charge carrier concentrations within the MIEC and allows to study the complex interplay with various transport components (drift, diffusion) and with the local distribution of surface reaction rates.In a similar way, the model also provides a quantitative description of species concentrations in the gas phase based on the coupled effects of transport and surface reaction.In Fig. 9, the molar fractions of hydrogen and water vapour are shown in the stagnant gas layer and in the pores of the CGO-layer.The concentrations are fixed at the end of the stagnant gas layer (left side), which reflects the ideal convective transport above the stagnant gas layer due to the excess fuel supply without concentration gradients, as introduced in Section 2.3.6.As hydrogen is consumed in the CGO-layer, the mole fraction decreases towards and within the CGO-layer and vice versa for water vapour.The concentration gradients are higher within the CGO-layer, because of the lower effective diffusivity in the porous media due to the ''bulk'' microstructure effects (M-factor M pore ) and the Knudsen diffusion resistance.Moreover, the concentration gradients are non-constant in the CGO-layer because transport and reaction (HOR) of gas species are both present.In contrast, the concentration gradients are constant within the stagnant gas layer because gas species are only transported (no reaction).

Frequency dependent results for an initial parameter set
The linear perturbation in the frequency domain enables the calculation of the impedance spectra.The impedance according to the surface reaction Z SR , the impedance due to the transport of the charge carriers in the CGO-MIEC Z chrgtpt , the gas concentration impedance Z gas and the sum of the three are plotted as Nyquist plot in Fig. 10(a) and the imaginary part of the impedance in Fig. 10(b).
In the separate Nyquist plot Fig. 12(a), the charge transport impedance shows a characteristic mouse shape, which is the Fig. 9 Mole fractions of hydrogen x 1 and water vapour x 2 in the stagnant gas layer and in the porous CGO-layer.Fig. 10 (a) Nyquist impedance plot and (b) imaginary parts of the impedance for the surface reaction impedance Z SR , impedance associated with the charge carrier transport in the CGO-layer Z chrgtpt , gas diffusion impedance Z gas , decoupled gas diffusion impedance Z gas,decoupled and the total anode impedance Z anode,tot .Note that Z gas,decoupled is determined by a separate simulation without coupling to Z SR , while Z gas is coupled with Z SR , resulting in non-distinguishable processes in an experimental EIS analysis.
result of a frequency dependent optimization of the transport pathways illustrated in Fig. 12(b).The different aspects of this process are visible in more detail in Fig. 11, where profiles across the CGO-layer of the linear deviations from the dynamic equilibrium associated with the harmonic perturbation are presented for different frequencies for CCL and electrolyte, respectively.The linear deviations of the reaction rate R ˜eon and of the electron concentration c ˜eon show that the electrochemical activity is concentrated to the electrolyte interface with increasing frequencies.The transport pathways for electrons and vacancies are adapted, resulting in much lower average transport pathlength for the mobile vacancies with increasing frequencies.A measure for the mean transport pathlengths for electrons and vacancies at a frequency of 500 Hz are the colored surface areas in Fig. 11(c).The mean pathlength for electrons is about a factor of 3 higher than for vacancies at a frequency of 500 Hz.In contrast, the mean pathlengths for electrons and vacancies are almost identical at low frequencies o0.1 Hz.
The mouse-shaped spectra of the charge transport process (Fig. 11(a)) ends with a nonzero real part for the plotted frequency range up to 1 MHz.In an experimental EIS measurement, this offset of the charge transport resistance is hard to distinguish from the electrolyte resistance.Another half circle process would add to the transport impedance due to the dielectric capacitance of the CGO-layer.However, this process has a peak frequency of about 100 GHz and is far beyond the frequency range normally used in EIS measurements.We thus plot the impedance spectra only up to 1 MHz and therewith this process is not visible in the plots of the simulation as well, even if its physics are included in the Poisson-eqn (14).
For the initial parameter set, the gas concentration resistance and the charge transport resistance are in the same range, while the surface reaction resistance is significantly higher than the other two effects.The gas concentration and the surface reaction impedance are in the same frequency  This journal is © the Owner Societies 2021 range of about 0.05 Hz.From eqn ( 11)-( 13), ( 16) and ( 17) it follows that the (perturbation of) sink and source terms for gas species are directly coupled to the surface reaction.Therefore both impedance-components Z SR and Z gas have the same peak frequency.
To illustrate the dependency of the gas diffusion impedance on the surface reaction resistance/chemical capacitance impedance, a decoupled simulation was performed for the gas diffusion impedance.Therefore, the gas diffusion impedance was calculated separately with the following homogeneous source terms in the porous layer: Therewith, the isolated gas diffusion impedance process can be calculated independently from the surface reaction.The resulting process Z gas,decoupled is plotted along with the coupled processes in Fig. 10(b).The isolated gas diffusion impedance is shifted to a much higher frequency range of about 2 kHz compared to the coupled gas diffusion impedance Z gas .

Parameter set resulting in distinguishable gas impedance arc
To illustrate the relation of the gas impedance and the surface reaction resistance/chemical capacitance processes, a simulation with a parameter set where those processes can be distinguished (also in a simulation with full coupling) is performed.
The parameter set is listed in Table 3.Only the parameters varied with respect to Table 1 are listed.The major changes are (a) the use of a huge stagnant gas layer thickness resulting in higher gas diffusion resistance and capacitance and (b) the reduction of the charge carrier equilibrium concentrations and a reduction of the MIEC-layer thickness, which results in a smaller chemical capacitance.This situation could be possibly observed experimentally for a large area cell with an anode based on a conductive backbone impregnated with CGO, which will be discussed in more detail in Section 4.3.The results for this parameter set are summarized with the Nyquist-plot in Fig. 13(a) and the imaginary part of the impedance in Fig. 13(b).Due to the choice of the parameters, the surface reaction resistance/chemical capacitance processes is shifted to higher frequencies and the gas diffusion impedance process is shifted to lower frequencies.The two processes are now distinguishable in the total impedance Z anode,tot .To judge the remaining coupling of the processes, the decoupled gas diffusion impedance Z gas,decoupled , calculated in the same manner as described in Section 3.3 with eqn ( 60) and ( 61), are plotted as well in both diagrams.The decoupled gas diffusion impedance has a characteristic mouse shape form, which is also visible to some extent for the coupled gas diffusion impedance process Z gas .

Parameter study with different CGO-layer thicknesses
We vary the CGO-layer thickness in a range of 10-800 mm to study its influence on surface reaction impedance, charge transport impedance, gas impedance, total impedance and the distribution of the reaction.Fig. 14 shows profiles across the CGO-layer with the source term of electrons, which illustrates the spatial distribution of the surface reaction rates for different anode layer thicknesses.The surface reaction rate varies only slightly along the CGO-layer for small CGO-layer thicknesses.The larger the CGO-layer thickness, the more the reaction is shifted towards the interfaces CGO/electrolyte and CGO/CCL.For large CGO-layer thicknesses (e.g.400 mm) the reaction rate in the center of the CGO-layer is almost zero.
In Fig. 15, the real parts of the different impedance contributions obtained by linear perturbation are summarized.The charge transport resistance ASR chrgtpt and the gas diffusion resistance ASR gas increase and the surface reaction resistance ASR SR decreases with growing CGO-layer thickness.There is a minimal total resistance for a thickness in the range of 100-200 mm, which is a specific finding for this anode and thus depending on the chosen parameter set and microstructure.For small CGO-layer thicknesses, the surface reaction resistance scales with 1/L.For larger thicknesses, there is a deviation from this law and the resistance does not decrease anymore with growing layer thickness, as illustrated in Fig. 15 by plotting the 1/L-behaviour.In Fig. 16(a) the Nyquist plot and in Fig. 16(b) the imaginary impedance part of the total anode impedance are shown for selected CGO-layer thickness's L = [10,100,800] mm.For a small thickness, almost a perfect semi-circle results, while for a large thickness, a characteristic ''mouse'' shape can be observed.

Parameter study with different microstructures
We now investigate the effect of microstructure on DC-and ACperformances.Changing the microstructure will affect the effective transport properties (ion-and electron-transport in the solid and gas transport in the pores) via the M-factors for MIEC-and pore-phase (see eqn ( 1)-( 4)).For gas transport in fine-grained anodes, Knudsen diffusion is also directly affected by the pore radius (eqn (26) in Section C of the ESI †).Furthermore, varying the microstructure also has an impact on the exchange current density and associated reaction rates via the specific surface area (see eqn (57)).Experimentally, these effects are difficult to quantify separately because the underlying processes are interlinked with each other.For example, in coarse grained materials with low surface area it can be expected that the reaction is extended deeper into the anode layer in order to make use of more surface area.As a consequence also the transport distances of vacancies, electrons and gas species will change.Moreover, the same microstructure changes also have an impact on the M-factor of the MIEC-and pore-phase (e.g.effective diffusivities), which may increase or compensate the corresponding ASR-contributions.This example illustrates that microstructure effects in SOFC electrodes are always the result of a complex interplay.A thorough understanding of these coupled effects can be achieved with a modeling approach that combines DC-simulations, which provide distributions of currents and species concentrations across   the anode layer, with AC-simulations, from which the different ASR-components can be extracted.
In order to capture the influence of the microstructure, the performance of three anodes with different microstructures are presented.Virtual microstructure data are used in order to keep the focus on the model and for a better reproducibility of the simulation results.The virtual microstructures have been realized using GeoDict 27 software using mono-sized spherical particles with a random seed of 45.The microstructure parameters are summarized in Table 2 and the corresponding orthoslices are shown in Fig. 6.Note that the virtual microstructures are isotropic and also the intrinsic conductivity of doped ceria is isotropic as stated by Koettgen et al. 65 MS1 has a high porosity with a solid volume fraction of only f MIEC = 0.4 and is coarse grained with a monosized particle diameter of d mono = 1000 nm.MS3 is a dense (f MIEC = 0.75) and fine-grained (d mono = 400 nm) microstructure.MS2 is an intermediate microstructure (f MIEC = 0.6, d mono = 800 nm), which is used as standard microstructure for all simulations unless otherwise specified.
Effective transport properties for the different microstructures are determined in two ways: (a) indirectly, using numerical simulation on the voxel grid and (b) directly, via morphological analysis as described in Section 2.3.2.The results are listed in Table 2.For the direct numerical analysis on the voxel grid, GeoDict 27 software was used.The M-factors for the effective diffusivity M sim pore and for the effective conductivity M sim MIEC were calculated with bulk (Laplace) diffusion using the DiffuDict-module.Note that the equations for the diffusion and drift (i.e.electric conduction) are mathematically equivalent.The M-factor for the effective Knudsen Diffusivity M sim pore,Kn and the characteristic pore diameter d pore for the Knudsen diffusion were calculated with a random walk algorithm of the DiffuDict module.The flow permeability k sim flow is calculated using the Stokes (LIR) solver of the FlowDict module with an inand outflow region of 100 voxels.
The M-factors for the effective diffusivity in the pore and the MIEC-phase were determined in both ways, indirectly (a) and directly (b).The ratio's M sim MIEC /M pred MIEC and M sim pore /M pred pore of the values extracted from simulations on the voxel mesh to the values predicted by morphological analysis are reported in Table 2.They show a deviation in the range of AE7%, except for MS 1 where a larger relative deviation is observed for the small absolute values of the MIEC-phase.For the parametrization of the FEM-model, the simulated properties were used because they are expected to provide the higher accuracy.However, the predicted M-factors are accurate enough to study the impact of different microstructure characteristics on the effective properties and thereby to identify the dominant microstructure effects.The microstructure properties t, b, e pore and f MIEC and the corresponding components t À4.39 , b 0.37 , f MIEC 1.15 and e pore 1.15 for the computation of the M-factor according to eqn (2)-( 4) are reported in Table 2.
The DC-results for three different microstructures are plotted in Fig. 17 as profiles across the porous CGO-layers showing distributions of (a) molar fraction of hydrogen, (b) potential drop and (c) reaction rates.The different contributions to the anode ASR as a function of the CGO-layer thickness obtained from AC simulations are plotted in Fig. 18.The highly porous and coarse grained MS1 shows the smallest gradient for the hydrogen molar fraction x 1 (Fig. 17(a)) and the highest potential drop along the CGO-layer (Fig. 17(b)).In accordance, MS1 shows the smallest gas diffusion resistance ASR gas (Fig. 18(a)) and the largest charge transport resistance ASR chrgtpt (Fig. 18(b)) of the three studied microstructures.The dense and fine-grained MS3 shows the steepest gradient for the hydrogen molar fraction x 1 (Fig. 17(a)) and the smallest potential drop along the CGO-layer (Fig. 17(b)).Accordingly, MS3 shows the largest ASR gas (Fig. 18(a)) and the smallest ASR chrgtpt (Fig. 18(b)).
The distribution of the reaction zone is shown in Fig. 17(c) for the different microstructures.For MS1, the electrochemical activity is enhanced towards the electrolyte and attenuated towards the CCL.The inverse behaviour can be observed for MS3.The surface reaction resistance ASR SR (Fig. 18(c)) is largest for MS1 and smallest for MS3, and therewith inverse  proportional to their specific surfaces.In Fig. 18(d), the total area specific resistance ASR tot is plotted.The highly porous MS1 shows the largest ASR tot for all calculated CGO-layer thicknesses.The dense, fine-grained MS3 shows the smallest ASR tot up to a CGOlayer thickness of about L = 75 mm, while for larger CGO-layer thickness's the intermediate MS2 shows the lowest ASR tot .
In addition, the simulation MS2 + Ni is plotted along, where the rate constant k 0 is increased by a factor of three, approximating the effect of additional fine dispersed but nonpercolating Ni-particles on the surface of the CGO.This case shows the lowest surface reaction resistance ASR SR (Fig. 18(c)) and also the lowest total resistance ASR tot (Fig. 18(d)) for all the reported CGO-layer thicknesses.The ASR gas and ASR chrgtpt for MS2 with Ni differs only very slightly in comparison with MS2 without Ni and are therefore not shown.

Charge transport resistance
The charge transport resistance is basically an ohmic resistance and therefore, one could expect to only observe a real impedance contribution.However, there is an imaginary component, because the transport pathways change with frequency, as illustrated in Fig. 12(b).With increasing frequency, the statistical mean transport pathlength for vacancies is minimized at the expense of a longer mean transport pathlength for the more mobile electrons.This results in a lower total charge transport resistance above the characteristic frequency of 0.45 Hz for the charge transport impedance.This can be explained with the nature of the chemical capacitance.The chemical capacitance is a volumetric effect in contrast to the dielectric capacitance, where the charges are stored at the interfaces.To access the chemical capacitance, the charge carriers need to travel for some distance into the bulk of the CGO-layer.Above the characteristic frequency of 0.04 Hz, the surface reaction is more and more bypassed by just buffering the charge carriers in the chemical capacitance of the CGO-layer.Note that the characteristic frequency of 0.45 Hz for the charge transport impedance is considerably higher than the characteristic frequency of 0.04 Hz for the surface reaction/chemical capacitance process as shown in Fig. 10(b).
The different aspects of this processes can be discussed in more detail by considering the profiles across the CGO-layer of the linear deviations from the dynamic equilibrium associated with the harmonic perturbation for different frequencies in Fig. 11.In Fig. 11(b) the profiles of the linear deviation of the electron concentration c ˜eon , show a higher depletion towards the electrolyte.Also the deviation from the dynamic equilibrium of the surface reaction rate is higher towards the electrolyte as shown in Fig. 11(a).The change of the distribution of pathlengths for electrons and vacancies with increasing frequency is shown in Fig. 11(c) by the normalized deviations of the electron and vacancy fluxes J ˜eon,norm and J ˜vac,norm .The colored areas are measures for the mean transport pathlengths at a frequency of 500 Hz, which is about a factor of 3 higher for the electrons as for the vacancies.At low frequencies ( f o 0.1 Hz), the pathlengths are approximately equal for electrons and vacancies.These results can be explained as follows.With increasing frequency, the chemical capacitance and associated storage of charges is shifted closer to the electrolyte interface.In this way, the mean transport distance for vacancies is reduced and at the same time, the surface reaction at longer distances from the electrolyte interface is suppressed.
By minimizing the transport pathlength of the less mobile vacancies associated with by-passing of the surface reaction, the transport resistance is decreased, resulting in a mouseshaped spectrum for the charge transport impedance with a remaining real part for the plotted frequency range of up to 1 MHz (Fig. 12(a)).In a measured impedance spectrum, this remaining charge transport resistance (visible in the high frequency limit) might not be distinguishable from the other ohmic series resistances such as the resistance of the electrolyte.However, the frequency-dependent optimization of the charge transport pathways can only take place when the surface reaction is mainly bypassed and a smaller active part of the CGO-layer is needed to provide sufficient reaction surface and chemical capacitance.Therefore, the peak frequency of 0.45 Hz for the charge transport process is considerably higher than the peak frequency of 0.04 Hz for the surface reaction resistance/ chemical capacitance process as shown in Fig. 10(b).
Note that the CGO-layer also owns a very small dielectric capacitance, where the charges are stored at the interfaces to the CCL and electrolyte respectively.To access this dielectric capacitance no transport into the bulk is necessary and therefore this additional charge transport process would end at zero in the high frequency limit.However, the characteristic frequency of this process is in the range of 100 GHz and therewith far beyond the frequency range measured in classical impedance spectroscopy and will not be observed in experiments.
The charge transport resistance depends also strongly on the microstructure.For the studied microstructures MS1-MS3, the microstructure factor varies in-between M MIEC = 0.04660-0.4826.Accordingly, the voltage drop differs remarkably for the different microstructures, as shown in Fig. 17

Distribution of the reaction zone
The distribution of the reaction zone and associated concentration profiles are driven and controlled by the following coupled processes and effects: Charge transport by drift and diffusion of electrons and vacancies in the CGO MIEC.Gas transport within the porous CGO-layer.Surface reaction resistance for the HOR.The combined effects of gas and charge transport limitations (but also taking into account the coupled HOR resistance) can be illustrated by comparing the impact of microstructure variation on local distributions of species concentrations, electrical potential and surface reaction rates (Fig. 17).For the highly porous MS1, the gas transport in Fig. 17(a) is less limiting than the charge transport in Fig. 17(b).Therefore, there is a higher electrochemical activity towards the electrolyte in Fig. 17(c) in order to reduce the transport path for the less mobile vacancies.These results are in very good agreement with the study of Nakamura et al., 20 where a higher electrochemical activity towards the electrolyte in a CGO anode was reported from experiments in a semiquantitative way.For the very dense and micro-porous MS3, where the gas transport is essentially governed by Knudsen diffusion and the conductivity is comparably high, the situation is inverted.Accordingly, an electrochemical more active zone towards the CCL Fig. 17(c) is observed, where the transport distance for the gas is low.For MS2, the gas transport and charge transport resistances are in the same range, resulting in a similar electrochemical activity at both interfaces, CGO/CCL and CGO/electrolyte.However, the reaction is not homogeneously distributed but is enhanced towards the interfaces.This is because of the complex mechanism of the drift-diffusion transport of the charge carriers.The charge transport by diffusion is enhanced by the concentration gradient as visible in Fig. 7(a).Therewith, the

Gas concentration impedance and surface reaction impedance
The gas concentration impedance process can be classified by the following contributions: Gas conversion impedance, especially relevant for large area cells.Gas diffusion impedance in the stagnant gas layer above the electrode.Gas diffusion impedance within the porous electrode, where bulk and Knudsen diffusion (described with dusty-gas model) can be relevant.The gas conversion impedance is neglected in our study, as we aim to describe a button cell system with an excess hydrogen flow rate.However, a very large stagnant gas layer might have a similar effect as a gas conversion impedance, as will be discussed later in this section.
The magnitude of the gas diffusion impedance in the stagnant gas layer depends on the thickness of the stagnant gas layer and therewith on the button cell setup and the flowrate.For our standard simulations, a quite small stagnant gas layer thickness of L Stag = 0.2 mm was used, which was estimated for a particular setup (see ESI, † Section F).The effect of a larger stagnant gas layer will be discussed later in this section.
The magnitude of the gas diffusion impedance in the porous electrode depends on the microstructure and on the thickness of the electrode.While the Knudsen diffusion coefficient is directly proportional to the pore size, the bulk diffusion depends on the dimension-less M-factor.Therefore, in contrast to the Knudsen diffusion, the bulk diffusion does not change when a certain microstructure is scaled (e.g. by changing the coarseness but at the same time maintaining volume fractions and particle and pore shapes).Both, pore sizes and M-factors vary strongly for the three microstructures and thus, these microstructure variations also have an impact on the effective Knudsen and bulk diffusivities.The Knudsen number reflects the ratio of bulk and Knudsen diffusion.We distinguish between the intrinsic Knudsen number Kn 0 = D 0 /D Kn,0 , defined as the ratio of the unitary bulk diffusion coefficient to the intrinsic Knudsen diffusion coefficient and the effective Knudsen number Kn eff = D eff /D Kn,eff , defined as the ratio of the effective unitary bulk diffusion coefficient to the effective Knudsen diffusion coefficient (see ESI, † Section D).The intrinsic Knudsen diffusion coefficient Kn 0 varies significantly for the three studied microstructures (Kn 0 o 1 for MS, Kn 0 4 1 for MS2 and Kn 0 c 1 for MS3).The Knudsen numbers are specific for different species and are listed in Table 2.Note that a Knudsen number of Kn 0 = 1 would reflect an equal contribution of bulk and Knudsen diffusion to the gas diffusion resistance.As the microstructure hinders the transport by bulk and Knudsen diffusion in a different way, the effective Knudsen numbers Kn eff are slightly larger than the intrinsic Knudsen numbers.
Our study suggests that in porous CGO anodes, the gas concentration impedance is only one of several processes that affect the low frequency spectrum.The RC-behavior of the chemical capacitance and surface reaction resistance might even be the larger contribution to the total impedance.These two low-frequency contributions (i.e.gas concentration impedance and surface reaction resistance/chemical capacitance impedance) are coupled with each other in two ways.There is a weak coupling between the two processes, as the concentration change due to the gas diffusion involves also a change in oxygen partial pressure, which impacts the exchange current density according to eqn (57).However, this effect is only significant if the gas diffusion limitation has a large impact on the p O 2 .This is the case for a very dry fuel supply, a setup with a large stagnant gas layer thickness or for thick, dense and fine porous electrodes.The second coupling becomes clear by comparing the characteristic frequencies of the processes.As already mentioned, the gas diffusion impedance and the surface reaction impedance show exactly the same characteristic frequency, which can be observed in Fig. 10(b).This can be understood in the way, these processes are coupled.An RC-like arc in an impedance spectrum always results by bypassing the resistance by just oscillating the charge in the capacitor.For the surface reaction impedance that means that the surface reaction resistance is bypassed above the characteristic frequency of the process.Significantly above the characteristic frequency there is no change in the surface reaction due to the applied perturbation at all, because the harmonic current is entirely stored in the chemical capacitance.If there is no effect on the reaction anymore, there is also no effect on the gas concentration at the electrode surface.Furthermore, if the concentration does not change anymore, also the gas diffusion impedance is bypassed.Due to this strong coupling associated with a frequency dependent bypassing, the two processes cannot be distinguished directly in an experimental EISspectra.For illustration, the decoupled gas diffusion process Z gas,decoupled is plotted in Fig. 10(b).The characteristic frequency f 0 of a process is generally given by f 0 ¼ 1 2pRC with associated resistance R and capacitance C. Because the process is now not cut off by the vanishing surface reaction, the arc is built by the capacitance of the gas transport process itself, while the associated resistance of the gas diffusion remains the same.As the gas transport capacitance is much smaller than the chemical capacitance, the process Z gas,decoupled is located at a much higher characteristic frequency around 2 kHz compared to the coupled gas impedance Z gas .In Section K of the ESI, † experimental examples from literature illustrating this observation are discussed and the order of magnitude of the characteristic frequency is confirmed with an analytical estimation.
With the parameter set of Section 3.4, a situation is constructed, where the gas diffusion impedance has a lower characteristic frequency than the chemical capacitance and associated surface reaction process.The huge stagnant gas layer results in a high 'RC' product for the gas diffusion process (larger diffusion resistance due to larger diffusion length and larger capacitance due to larger gas volume).In this situation, the two processes are clearly distinguishable from each other in This journal is © the Owner Societies 2021 the total anode impedance Z anode,tot in Fig. 13(a) and in Fig. 13(b), respectively.A weak coupling of the two processes still persists, which becomes clear by comparing Z gas with the decoupled Z gas,decoupled impedance in the imaginary part plot in Fig. 13(b).Between 10-100 Hz, Z gas,decoupled impedance still shows a small contribution, resulting in characteristic mouse shape in the Nyquist plot in Fig. 13(a), which is typical for diffusion limited processes.However, the coupled Z gas process is cut above 10 Hz because of the vanishing surface reaction in Fig. 13(b), which results in a distorted mouse shape in the Nyquist plot in Fig. 13(a).
The goal of the simulation with this parameter set was to illustrate 'theoretically' the dependency of the gas impedance with the surface reaction resistance/chemical capacitance process.However, the same simulation with this dedicated parameter set can also be used to briefly discuss the correspondence with two potentially existing physical cell setups, as follows: (a) The huge stagnant gas layer thickness (which leads to a large Z gas with a relatively low characteristic frequency) could e.g.represent a flow-situation, where a button cell is clamped between two impermeable disks with an axial fuel supply.Therewith, the mass transport across the button cell needs to be largely achieved by diffusion.If a large amount of fuel is provided, the gas conversion impedance can be avoided, but a large gas diffusion impedance remains.Measurements with a setup of this kind have been reported by Aravind et al. 66 However, the effect of this setup is similar to the gas conversion impedance, especially important for large area cells, where the order of magnitude of the impedance for this parameter set is easily reached.
(b) The lower charge carrier concentration and the associated lower chemical capacitance (with a relatively high characteristic frequency) could be realistic in a composite material, where e.g. a conductive matrix is impregnated by CGO.Because the chemical capacitance is a volumetric effect, the lower CGO volume fraction would then result in a lower capacitance.Of course, a lower charge carrier concentration could possibly also be observed for a different MIEC-materials.In summary, these hypothetical cases used for the illustration of the gas impedance process could possibly be observed in a large area cell or in an anode with different material composition.
In the following we will discuss how low frequency (LF) arcs in impedance spectra of MIEC anodes are interpreted in the literature.In general, there are two main LF processes, which are either associated with the gas concentration impedance or with the surface reaction resistance/chemical capacitance impedance.As will be discussed subsequently, depending on the authors but also depending on the cell architecture, the low frequency arc is usually interpreted as being either gas impedance or surface reaction impedance, but rarely as a mixture of both.For example, for studies with patterned thin-film electrodes as e.g.performed by Chueh et al., 22 the low frequency process is attributed to the surface reaction resistance/chemical capacitance impedance process.In thin-film electrodes (massive, no pores) the gas concentration impedance is negligible, simply because the surface area is relatively small and thus the surface reaction resistance is dominant.Additionally, as no net current is drawn in the symmetric cell setup, also the possible effect of p O 2 change along with concentration changes due to gas diffusion is probably negligible for the whole common parameter range.Furthermore, in thin-film electrodes transport distances in the MIEC are usually very short so that also the charge transport resistance becomes negligible (unless especially designed to enforce a significant charge transfer contribution as e.g.studied by Chen et al. 13 ).Therefore, thinfilm electrodes typically exhibit EIS-spectra with a perfect semicircle and low characteristic frequency.
Riegraf et al. 19 suggested for their study of an Ni-CGO10 anode (4 Â 4 cm 2 active area) that there is a thermally activated electrode process in addition to a gas conversion process for the low frequency arc around 0.2 Hz.We suggest that this electrode process is the process of surface reaction resistance/chemical capacitance described in this work and hence in these anodes we are faced with a mixture of both LF processes.
Price et al. 5 assigned the low frequency process in the impedance spectra of a button full cell (1 cm 2 active area) with a CGO/Ni-impregnated LSCT A backbone anode as a gas concentration effect.Their argument is that the process is not thermally activated, which is in fact expected for a gas impedance process.However, in the work of Riegraf et al. 19 with Ni-CGO anodes, the thermal activation was only visible in the low frequency arc for temperatures below 750 1C.So the thermal activation might just not be visible in the studied temperature range of 800-900 1C.Therefore, there might be a surface reaction resistance/chemical capacitance process included in the low frequency arc.
Burnat et al. 63 performed a study with composite MIEC anodes of variable LST-CGO ratios (in button cells).The resulting low frequency effect was attributed to a chemical capacitance process, which in our nomenclature would be a surface reaction resistance/chemical capacitance process.They argued that this process depends on the LST-CGO ratio and also shows a thermal activation, which is very unlikely for a pure gas concentration impedance process.However, there will be for sure a gas concentration impedance effect included in the low frequency arc, but this effect might be small in comparison with the surface reaction resistance/chemical capacitance process.For the anode compositions with CGO contents above 50%, which show a good performance, the gas diffusion impedance may become significant.In such composite anodes with low surface resistance, the gas impedance can become even more significant when the nitrogen content in the fuel gas composition is increased, because this leads to a stronger gas diffusion impedance.However, the size of the gas diffusion impedance depends strongly on the actual test setup and cannot be judged here for the experiments by Burnat et al. 63 Lomberg et al. 67 performed EIS measurements for CGO10 anodes infiltrated with 14 wt% Ni.They interpreted a low frequency arc measured at T = 750 1C as a combination of gas diffusion impedance and chemical capacitance.The interpretations of Lomberg et al. are compatible with our model-based interpretations, according to which the low frequency arc would be attributed to surface reaction resistance/ chemical capacitance process in combination with a gas diffusion process, where the latter is forced to the same frequency (i.e.Z gas ), as explained before.Aravind et al. 66 measured the impedance spectra of a Ni-CGO anode using a setup with a large stagnant gas layer and identified a gas diffusion impedance process around a frequency of 0.5 Hz and a second, not further characterized process around 8 Hz.In this configuration, a low frequency process still persists, which therewith has to be attributed to the surface reaction resistance/chemical capacitance process.
Nakamura et al. 38 performed EIS measurements for CGO10 anodes and associated the low frequency arc to a surface reaction resistance/chemical capacitance processes because of it's clear thermal activation.Although the effective transport properties of the microstructure were not characterized in detail within this publication, the reported SEM-pictures show a relatively high porosity with pores sizes about 1 mm.It thus seems reasonable that the gas impedance might be small compared to the surface reaction resistance/chemical capacitance process.
Using a special setup, Primdahl and Liu 4 studied a Ce 0.6 Gd 0.4 O 1.8 anode, mounted between auxiliary electrodes to avoid significant impedance contributions from gas diffusion outside the porous electrode structure.The remaining gas diffusion in the porous electrode was found to be negligible compared to the other processes.With this procedure, they assured the reliable interpretation of the EIS-spectra in an experimental way.
Riegraf et al. 58 performed an EIS-study on symmetrical Ni-CGO10 electrodes for different temperatures and the gas compositions.They attributed a low frequency arc around 0.1-1 Hz as a combination of a surface process attributed to a charge transfer and a gas impedance process.The gas conversion and gas diffusion effects of the setup were quantified using another measurement set with Ni-YSZ anodes, where the gas impedance could well distinguished from other processes.The gas diffusion impedance within the electrode was assumed to be negligible, which might be appropriate for the relatively porous electrode with a thickness of 25 mm.For this particular case, they could isolate the surface reaction resistance/chemical capacitance process by using a simple model for the gas impedance.
Linder et al. 18 performed EIS measurements for a Hexis short stack with Ni-CGO anodes for different fuel flow rates.The low frequency process o0.5 Hz is associated with a gas concentration process, as this process is strongly affected by a variation of the fuel amount.In the frequency range 40.5 Hz, two smaller not fully formed arcs are associated with the anode and cathode electrodes.This setup is comparable with our simulation with very large stagnant gas layer shown in Fig. 13 in Section 3.4.The very large stagnant gas layer has a similar effect like the gas conversion impedance due to fuel consumption along the cell in the setup of Linder et al. 18 In contrast to our simulation study, which is parametrized in order to achieve well distinguishable gas and surface reaction resistance/ chemical capacitance processes, the two processes largely overlap in the EIS-spectra of Linder et al. 18 However, we suggest that the first small and incomplete arc in the frequency range 40.5 Hz can be associated with the surface reaction resistance/chemical capacitance process.
Overall, this short discussion of literature highlights that the interpretation of the low frequency arc for MIEC-based anodes (mostly CGO or CGO-composites) strongly depends on the cell architecture (i.e.thin-film electrodes, button cell, large area cell/stack).In thin-film electrodes, surface reaction resistance/ chemical capacitance is the dominant LF phenomena and gas impedance tends to be negligible.In button cells, the LF-arc results in most cases from a combination of surface reaction resistance/chemical capacitance and gas impedance.Thereby the relative contribution of each process depends on several aspects, such as cell architecture, anode composition, microstructure and working conditions (fuel, temperature).Finally, in large area cells, the gas conversion impedance typically dominates over the surface reaction/chemical capacitance impedance.

EIS-spectra of CGO-based anodes for different applications
As demonstrated before, the CGO-anode arc, which appears as only one single process in the EIS-spectra, deconvolutes in three processes: (a) surface reaction resistance/chemical capacitance process, (b) gas diffusion impedance process and (c) charge carrier transport impedance process.Because of the strong coupling and overlapping of these processes, they can hardly be captured by analysing experimental EIS-spectra only.The contribution of the three processes to the total anode process depends on the actual cell design.In the following we provide some application examples and discuss the importance of the different contributions.Please note that the specified parameters are not absolute but depend on the actual material system and microstructure.If for example an additional catalyst like Ni or Pt is used, the surface reaction resistance will be considerably lower, and the transport of gas and charges will become important already at smaller CGOlayer thicknesses.The same is true for the different microstructures.For very dense and fine-grained microstructures with a relatively high surface area, the surface reaction resistance and the charge transport resistance decrease, while the gas diffusion resistance increases and will be relevant already for lower CGO-layer thicknesses and so on.Therefore, the different depicted cases will be observed for different absolute parameter ranges for different microstructures and material systems.
(a) For very thin CGO-electrodes on the button cell level, the surface reaction resistance/chemical capacitance process dominates the anode process and the EIS-spectra will form an almost perfect semi-circle like in Fig. 16 for CGO-layerthickness of L = 10 mm.The small impedance contribution This journal is © the Owner Societies 2021 from the gas and charge transport are negligible.An extreme case of this situation are thin-film electrodes.
(b) For CGO-electrodes with an intermediate thickness on the button cell level, the surface reaction resistance/chemical capacitance process still dominates the anode process, but the impedance contribution from the gas and charge carrier transport are no more negligible.However, the gas impedance is forced to the same frequency range because the gas impedance process is cut off when the surface reaction vanishes, which is the case above the characteristic frequency of the surface reaction resistance/chemical capacitance process.The two processes are not distinguishable by studying the total anode impedance experimentally.The small contribution of the transport resistance associated with relatively short transport pathways are usually not notable in the experimental spectra for high temperature applications in the range of T = 850 1C.However, for lower operating temperatures, the charge transport will be more significant according to the temperature dependency of the conductivity of eqn ( 51) and (52).If the anode EIS-spectra does not overlap with the cathode-process, a distortion of the arc-like process in a mouse-shaped end towards the higher frequencies might be visible like in Fig. 16 for CGO-layer-thickness of L = 100 mm.This configuration is typical for electrolyte supported button cells.
(c) For very thick CGO-electrodes, the transport resistances can become dominant because of the long transport ways for the charges and the gas species.For dense and very finegrained electrodes, the gas impedance can already become limiting for intermediate layer thicknesses due to the transport limitation by Knudsen diffusion.It is worth noting that a part of the charge transport resistance is included in the anode process arc, while another part is only visible as an ohmic resistance as discussed in Section 4.1.This fact should be kept in mind if a MIEC-spectrum is analysed.However, the described behavior is only valid for the case of CGO-based anodes, where the CGO-phase is the only percolating phase in the composite.The EIS-spectra for this case will look like in Fig. 16 for CGO-layer thickness of L = 800 mm with a pronounced mouse shape and a large ohmic offset.This configuration might be important for anode supported button cells.
(d) For large area cells, the gas impedance process, which then is dominated by the gas conversion impedance, is usually the dominant part of the anode process.If the characteristic frequency of the gas impedance is below the characteristic frequency of the surface reaction resistance/chemical capacitance process, the two processes can even become distinguishable in experimental EIS-spectra, such as the one in Fig. 13.Note that for the case in Fig. 13(a), a dedicated parameter set for button cells is used, which is characterized by a huge stagnant gas layer.This case can be experimentally realized, but is not of importance for technical applications.However, the behaviour of this case is very similar to the gas conversion impedance configuration, relevant for SOFC large area cells.

Towards a systematic optimization of MIEC electrodes with digital materials design (DMD)
The optimization of MIEC electrodes is a very complex task with many parameters involved.As an example, in Fig. 18 different contributions to the anode ASR for different microstructures as a function of the layer thicknesses are shown.
The total anode ASR is shown in Fig. 18(d) and is largest for MS1 for all CGO-layer thicknesses.For a CGO10 anode with MS1 microstructure, the surface reaction resistance is too high for the porous and coarse grained microstructure with its low specific surface area.With increasing thickness of the CGOlayer also the surface area increases and therewith the surface reaction resistance is reduced.At the same time, this causes a large charge transport resistance (Fig. 18(b)) because of the low M-factor for the MIEC (i.e.low effective ionic and electric conductivities) in combination with long transport pathways for charge carriers.
In contrast, the very dense and coarse grained MS3 has a large surface area and thus a low surface reaction resistance.For small CGO-layer thicknesses, MS3 therefore shows the best performance of the three microstructures.However, for the dense microstructure with Kn 4 1, the gas diffusion resistance in the pores is very high even for intermediate CGO-layer thicknesses as can be seen in Fig. 18(a).Note that the gas diffusion resistance does not grow linearly with the CGO-layer thickness.This is because of the shift of the reaction zone towards the CCL as visible in Fig. 17(c), which limits the additional gas transport distance.
The intermediate microstructure MS2 shows a higher total ASR for small CGO-layer thicknesses (compared to MS3), because of its lower specific surface area resulting in higher surface reaction resistances.However, as the charge and gas transport resistance are both not too high the transport resistances increase moderately with increasing CGO-layer thickness as shown in Fig. 18(a) and (b).Therefore, the MS2 shows the lowest total ASR for L 4 75 mm.Of course, the trade-off changes as well if the material system is changed.If for example a fine dispersed catalyst like Ni is added on the pore surface of the CGO-layer, the surface reaction resistance decreases, which also decreases the optimal layer thickness for the same microstructure (Fig. 18(d)).
The effective transport properties in the porous electrode can be described in a reliable way by running numerical simulations on the voxel mesh obtained from 3D imaging or from virtual microstructures.However, the effective properties (and the M-factors, respectively) do not provide a detailed understanding of the microstructure influence on the various transports.Additional information about possible optimization approaches can be deduced from morphological analysis described in Section 2.3.2.In order to illustrate the detailed level of morphological information, we compare in the following the various microstructure limitations for transport in the solid phase of MS1 with those in the pore phase of MS3.Apart from the rather intuitive dependency of the transport properties on the phase volume fraction, the constrictivity b and the tortuosity t document additional transport restrictions due to bottlenecks and tortuous transport pathways, respectively.For example, the M-factor of the MIEC-phase M sim MIEC = 0.04660 for MS 1 with f MIEC = 0.4 shows a considerably lower value than the M-factor of the pore-phase M sim pore = 0.09292 for MS 3 with e pore = 0.25, despite the considerably larger phase volume fraction.This can be explained by the fact that the constrictivity b MIEC = 0.2010 of MS 1 is much lower than the constrictivity b pore = 0.4756 for MS 3.This is in agreement with the visual impression of limited connectivity of the only slightly overlapping particles in Fig. 6.In addition, the tortuosity t MIEC = 1.2114 of MS 1 is larger than the tortuosity t pore = 1.1136 of MS 3.This fact is less obvious from the visual impression from Fig. 6.Because of the different weights of the b 0.37 and t À4.37 in eqn (2), the tortuosity effect is even slightly higher for the transport limitation than the constrictivity effect as reported in Table 2.
Note that the virtual microstructures employed do not correspond to typical MIEC-electrodes.Such simplified microstructure were chosen for the purpose of illustration and easy reproduction.Realistic microstructure data can be obtained e.g. by FIB-tomography.Moreover, experimentally obtained microstructures can be varied by using microstructure simulation.For this purpose stochastic, digital twins of the experimental microstructures from tomography are created.They serve as a basis for a stochastic model (as e.g.performed by Neumann et al. 68 for Ni-YSZ anodes), which is then capable to provide a large number of different virtual microstructures with realistic properties (e.g.parametric variations of porosity or particle sizes).The output from these microstructure simulations can then be combined with our numerical model for MIEC anodes, which opens new possibilities for virtual materials testing (VMT) and digital materials design (DMD).It must be emphasized that a systematic optimization of MIEC-anodes is only possible, if the material system with its associated bulk material properties, the microstructure and the design-parameters as e.g. the anode layer thickness are combined into a multiphysics model including all relevant physical processes.Such a systematic digital materials design (DMD) process for MIEC anodes could have the following procedure: the bulk properties of a specific material system are collected with literature data or measurements.Moreover, a specific cell-design and production process is chosen.A couple of button cells are then fabricated and characterized with EIS measurements, where the missing material parameters can be fitted.In a post mortem analysis, the microstructure can be reconstructed using 3D-tomography (e.g.FIB-tomography) and its effective properties can be determined using morphological analysis or direct simulation on the voxel microstructure using e.g.GeoDict. 27All the material and microstructure parameters can then be inserted in our MIEC model, which includes all relevant physical effects.The parametrized model can be validated and calibrated by comparison with the EIS-spectra.In the next step the numerical model can be combined with a microstructure model, which leads to a framework for digital materials design (DMD).Thereby, digital twins are deduced from the real microstructures as realistic start or edge-points for achievable microstructures.Virtual microstructures can then be used to vary microstructure properties and to analyze their impact on the cell performance (virtual materials testing VMT).Also the cell-design parameters (like e.g. the anode layer thickness, etc.) or the material system (e.g.different MIEC materials or addition of catalysts) can be varied.The MIEC model can then be used to predict the corresponding cell or anode performances.If one or more promising design-points are found, the realizability of these design points can be checked and validated by further experimental work.

Conclusions
A comprehensive model revealing the complex physicochemical processes for CGO-based anodes on the button cell level is presented in this paper.The 1D-FE model can be operated in DC and AC mode and therewith the behaviour of the normal DC operation can be directly linked to the corresponding EIS-spectra.The fully implemented Nernst-Planck Poisson equations enables an appropriate description of the transport of electrons and vacancies in the CGO-MIEC.For fine-grained microstructures, the gas diffusion is limited by both, bulk and Knudsen diffusion.To capture these effects, a dusty-gas model is implemented, which is believed to be the most appropriate description for combined bulk and Knudsen diffusion.For the parametrization, literature data for CGO10 was used to estimate realistic values for the equilibrium charge carrier concentration, which has a strong influence on the chemical capacitance, the ionic and electronic conductivity and for the surface reaction resistance.It must be emphasized that it is beyond the scope of this paper to present a match with experimental measurements.Therefore, in the present paper, the microstructure effects are demonstrated on simple, well reproducible virtual structures.
The CGO-anode arc, which appears as only one single process in the EIS-spectra, is shown to deconvolute in three processes (a) surface reaction resistance/chemical capacitance process, (b) gas diffusion impedance process and (c) charge carrier transport impedance process.Because of the strong coupling and overlap of these processes, they can hardly be captured by analysing experimental EIS-spectra only.Therefore, in the literature the anode arc, which typically has a low characteristic frequency (around 0.2-2 Hz), is often either depicted as gas impedance or as chemical capacitance process, without conclusive evidence.The contribution of the three processes to the total anode process depends on the cell design such as mechanical support (and associated thickness of anode layer) or design of gas distribution (and e.g.associated thickness of stagnant gas layer).In addition, also the microstructure and associated effective materials properties play an important role.This results in a large number of anode configuration, whereby the same material concept (i.e.CGO-based MIEC anode) can reveal very different performances characterized with different impedance spectra.The controversial interpretation of rate limiting processes in literature indicates that experimental methods alone may not be suitable for a reliable interpretation of such EIS-spectra from different anode configurations.Fortunately, our numerical model can be used to simulate the EIS-spectra for all these different MIEC/CGObased anode configurations.In addition, the model also enables us to identify the limiting effects and/or quantifying the resistive/capacitive contributions from different processes.The model can thus be used to illustrate the impact of different CGO-anode configurations (i.e.microstructure and associated effective properties, anode layer thickness, gas distribution system and stagnant gas layer thickness).Fig. 19 provides an overview of the simulated EIS-spectra and identifies the limiting processes for many different configurations discussed throughout this paper at an operating temperature of T = 850 1C.In the following we briefly emphasize four specific configurations, which represent the most important cases in For very low CGO-layer thicknesses in button cells or for thin-film electrodes, the surface reaction resistance is dominant.For intermediate electrode thicknesses, which are e.g.used for electrolyte supported button cells, the surface reaction resistance is still dominant, but the resistances due to gas and charge transport are no more negligible, even if they are hardly distinguishable in measured EIS-spectra.
For larger CGO-layer thicknesses, which might be used for anode supported cells, the resistances for gas or charge transport -dependent on the actual material system and microstructure -might be dominant.
For large area cells, the gas impedance process, which then is dominated by the gas conversion impedance, is usually the dominant part of the anode EIS-spectra and might be even distinguishable if the gas process shows a considerably lower characteristic frequency than the surface reaction resistance/chemical capacitance process.Although the impedance according to the transport of electrons and ions is basically an ohmic effect, it shows a mouse-shaped arc with a non-zero ohmic contribution at the high frequency limit as a result of changing transport pathways for higher frequencies (in combination with by-passing of the surface reaction).A part of the transport resistance is included in the anode process arc, while another part is only visible as an ohmic resistance adding to the electrolyte resistance.This behaviour should be kept in mind for an appropriate interpretation of experimental EIS-spectra for MIEC (CGO) electrodes.
With our MIEC model in DC mode, the extension of the reaction zone for different CGO-layer thicknesses can be illustrated.It is observed that thick CGO layers are also characterized by a large inactive zone.A good qualitative agreement with experimental studies from literature could be shown and a more precise mechanism for the electrochemically active region could be suggested.For very dense fine-grained microstructures, the electrochemical most active zone might be found towards the CCL while for porous, course-grained microstructures, it might be found towards the electrolyte.The extension of the reaction zone can have a direct impact on the anode optimization in order to reduce the ASR.It is obvious that a too small CGO-layer thickness will posses a very high surface reaction resistance and a too large CGO-layer thickness will own a very high transport resistance.The optimal layer thickness in between will depend on how the extension of the reaction zone will behave with increasing CGO-layer thickness.The optimal thickness thus depends on various effects and coupled processes, which can be predicted reliably with our MIEC-CGO model.
Our simulation study enables the appropriate interpretation of the EIS-spectra and the understanding of the complex physico-chemical processes, which is a crucial prerequisite for a systematic materials optimization for CGO-based anodes.The model thereby serves as a platform for the integration of information from bulk material properties, fabrication parameters, microstructure effects, material laws and operating conditions and relates them to the cell performance on button cell level.Hence, with the integration of detailed physicochemical properties into a single model framework, findings from basic and applied research, from theoretical models and experimental characterization can be directly used for a systematic optimization of CGO-based anodes.
We also discussed the methods of digital materials design (DMD) for the efficient optimization of anode materials and microstructures.In such a DMD framework, the MIEC (CGO) anode model is combined with other analytical and simulation methods (tomography, image analysis, stochastic microstructure modeling and characterization e.g. with impedance spectroscopy).The DMD framework allows virtual testing of a very large parameters space (e.g.many different virtual microstructures, different anode architectures/thicknesses, etc.) and predicting the corresponding anode performances.Such extensive parameter studies would not be achievable with experimental methods only.As a perspective for the future, the DMD approach that includes our MIEC-model will enable a systematic, fast and cost-effective research and development process for MIEC-anodes with improved properties.However, the DMD approach is not limited to fuel cells, but can be adapted for any application where microstructure effects are important, like e.g.batteries, filter materials, sensors, porous membranes, etc.
In a nutshell, the most important conclusions are the following: EIS-impedance only on experimental evidence is difficult to interpret.The presented model helps to understand important effects in MIEC-and especially CGO-based anodes and to distinguish different limiting processes as e.g. the gas impedance and the surface reaction resistance/ chemical capacitance process.The model shows that the dynamics of these two processes are coupled and hence occur in the same frequency range.
Based on the model it can be illustrated that configurations with the same anode material but with different cell architectures (thin film, button cell, large area cell/stack) result in very different impedance signatures.The model implies that the extension of the reaction zone also depends on the transport limitations.If the charge transport is limiting (e.g. for highly porous and course microstructures), the most reactive zone is towards the electrode/electrolyte interface.If the gas transport is limiting (e.g. for very dense and fine-grained electrodes), the most reactive zone is towards the electrode/CCL interface.Using effective properties from microstructures (e.g.impact of tortuosity) and intrinsic properties from materials (e.g.MIEC intrinsic conductivities, impact of catalyst on reaction kinetics) as input, the model provides a basis for a systematic anode optimization.Thereby it turns out that the microstructure needs to be optimized specifically for a given architecture and material system.E.g. for a thick anode the microstructure optimization needs to face the transport limitations, whereas in thin anodes sufficient reaction surface has to be provided.
Solid-phase constrictivity contribution to the predicted M-factor b MIEC 0Solid volume fraction contribution to the predicted M-factor f Ratio predicted to simulated solid-phase M-factor M sim MIEC /This journal is © the Owner Societies 2021

Fig. 4
Fig. 4 Schematic visualization of the crystal structure of CGO, with cerium ions in oxidation states Ce 4+ and Ce 3+ , gadolinium ions Gd 3+ , oxygen ions O 2À and oxygen ion vacancies V O .
O 2 ,ref = 10 À19 bar and a reference temperature of T ref = 850 1C (see ESI, † Section I for the detailed calculations and assumptions).2.4.4 Formulation of the surface reaction overpotential for CGO.In Section 2.3.4,the surface reaction rate was introduced in eqn (13) in a general Buttler-Volmer form.The surface reaction overpotential Z SR , which drives the surface reaction rate, shall now be formulated to capture the rate limiting step in a CGO-anode.The overall anode reaction reads:

Fig. 5
Fig. 5 Oxygen-nonstoichiometry of CGO10 as a function of temperature and oxygen partial pressure.Markers are manually extracted data from Wang et al., 54 the solid lines are the results of the fit function used in this work.

2 À1/ 4 -
11, 12, 14 and 22, the exchange current shows a p O dependency.This indicates that the exchange current density is proportional to the number of available charge carriers and hence with the oxygennonstoichiometry d, which shows a p O 2 À1/4 reported a power law of approximately p H 2 O 0.08 for the reaction rate of an Ni-CGO10 symmetric anode button cell.Note that the oxygen partial pressure also increases for increasing humidity, which has a negative effect on the reaction kinetics according to the p O 2 À1/4

Fig. 8
Fig. 8 (a) Deviation of the electron and vacancy concentrations with respect to the equilibrium concentration (no current), (b) distribution of the reaction along the CGO-layer, illustrated with the electron and vacancy source terms.

Fig. 7
Fig. 7 Profiles across a 100 mm thick anode layer (x-axis) under current load, illustrating the effects of transport by drift and diffusion on (a) molar concentrations of electrons and vacancies, (b) electric potential and electric field, (c) molar electron flux, (d) molar vacancy flux.CCL = current collector layer, ELY = electrolyte.The dashed lines indicate the location of the vertex of the charge carrier concentrations.

Fig. 11
Fig. 11 Profiles of local deviations from the dynamic equilibrium due to a harmonic perturbation across a 100 mm thick anode layer (x-axis) for different perturbation frequencies to illustrate the charge transport impedance: (a) electron reaction rate R ˜eon , (b) electron concentration c ˜eon and (c) normalized electron and vacancy fluxes J ˜eon,norm and J ˜vac,norm where the colored areas are measures for the mean transport pathlengths of electrons and vacancies (directions indicated with arrows) for the case of f = 500 Hz.Note: the linear deviations of the reaction rates and the charge carrier concentrations are very similar for electrons and vacancies and essentially differ by a scaling factor.Therefore, only the profiles for the electrons are shown in (a) and (b).

Fig. 12
Fig.12Visualization of the charge transport impedance process: (a) Nyquist plot of the charge transport impedance Z chrgtpt , (b) schematic illustration of the transport pathways of electrons and vacancies for low and high frequencies.At low frequencies, the reaction rate is almost homogeneously distributed over the CGO-layer, resulting in similar transport distances for electrons and vacancies.At high frequencies, only a small reaction zone in the vicinity of the electrolyte persists, resulting in longer transport distances of the more mobile electrons and therewith an increase of the total charge transport resistance.Note: the red and violet dots indicate typical locations of electron and ion sources (i.e.statistical mean values related to frequency dependent distributions, which are illustrated in Fig.11).

Fig. 13
Fig.13Nyquist (a) and imaginary part of the EIS-spectra (b) for a parameter set resulting in distinguishable impedance signatures from gas transport and surface reaction resistance/chemical capacitance processes.Surface reaction impedance Z SR , impedance associated with the charge carrier transport in the CGO-layer Z chrgtpt , gas concentration impedance Z gas , gas concentration impedance Z gas,decoupled decoupled from the surface reaction (separate simulation), and the total anode impedance Z anode,tot .

Fig. 15
Fig.15Different contributions to the anode area specific resistance (ASR) as a function of the CGO-layer thickness: surface reaction resistance ASR SR , resistance associated with the charge carrier transport in the CGO-layer ASR chrgtpt , gas diffusion resistance ASR gas and the total anode area specific resistance ASR anode,tot .

Fig. 16
Fig. 16 (a) Nyquist plot and (b) imaginary impedance plot of the total anode impedance for selected CGO-layer thicknesses.

Fig. 14
Fig. 14 Surface reaction rate (source term of electrons) across the CGOlayer for different CGO-layer thicknesses.

Fig. 17
Fig. 17Results from DC simulations plotted as profiles across the porous CGO-layers for three different microstructures showing distributions of: (a) molar fraction of hydrogen, (b) electric potential drop and (c) reaction rates.The three anodes have constant thickness but different microstructures (MS1 = high porosity and coarse, MS2 = intermediate, MS3 = dense, fine-grained).
Fig. 17Results from DC simulations plotted as profiles across the porous CGO-layers for three different microstructures showing distributions of: (a) molar fraction of hydrogen, (b) electric potential drop and (c) reaction rates.The three anodes have constant thickness but different microstructures (MS1 = high porosity and coarse, MS2 = intermediate, MS3 = dense, fine-grained).

Fig. 18
Fig. 18 Results from AC simulations for different microstructures.Different contributions to the anode ASR for different microstructures are shown as a function of the CGO-layer thicknesses: (a) gas diffusion resistance ASR gas , (b) resistance associated with the charge carrier transport in the CGO-layer ASR chrgtpt , (c) surface reaction resistance ASR SR and (d) the total anode area specific resistance ASR anode,tot . (b).

Fig. 19
Fig. 19 Summary of typical EIS-spectra of CGO-based MIEC anodes for cells on different scales (a) thin-film electrode, (b) button cell, (c) button cell with thick stagnant gas layer which shows a similar behaviour as large area cells, different microstructures MS and different anode layer thicknesses L.Effective transport properties vary with microstructure, which is quantitatively described by the M-factor (M pore = D eff /D 0 for the gas species transport, M MIEC = s eff /s 0 for the charge transport).The specific surface are S is a measure for the available two-phase boundaries for the surface reaction.For the thin-film electrode the parameter CS = d/L is a normalized measure for the contact spacing, where d is the contact spacing and L is the dense MIEC-layer thickness.The red coloured spectrum (button cell with small stagnant gas layer, MS3 and L = 10 mm) is plotted in each figure as a reference.The operating temperature is T = 850 1C for all cases.

Table 1
Standard parameter values used for the simulation.For the effective properties, the values for MS 2 listed in Table 2 are used as standard parameters vac 3.9228 Â 10 À9 m 2 s À1 Intrinsic diffusivity of vacancies Estimated from Steele 49 as documented in the ESI Section I D eon 3.3996 Â 10 À8 m 2 s À1 Intrinsic diffusivity of electrons Estimated from Steele 49 as documented in the ESI Section I k 0 5.5 Â 10 À4 A m À2 Surface specific exchange reaction rate constant Estimated from Gerstl et al. 57 as documented in the ESI a D a These parameters have been chosen according to illustrative aspects as discussed in Section 3.1.

Table 2
Microstructure parameters and effective properties from image analysis and transport simulations for the virtual microstructures MS1, MS2 and MS3 (Fig.6).sim = effective transport property determined by transport simulation, pred = effective transport property prediction determined by microstructure parameters from image analysis

Table 3
Parameter set resulting in distinguishable gas impedance arc