Computational engineering of the oxygen electrode-electrolyte interface in solid oxide fuel cells

The Ce0.8Gd0.2O2−δ (CGO) interlayer is commonly applied in solid oxide fuel cells (SOFCs) to prevent chemical reactions between the (La1−xSrx)(Co1−yFey)O3−δ (LSCF) oxygen electrode and the Y2O3-stabilized ZrO2 (YSZ) electrolyte. However, formation of the YSZ–CGO solid solution with low ionic conductivity and the SrZrO3 (SZO) insulating phase still happens during cell production and long-term operation, causing poor performance and degradation. Unlike many experimental investigations exploring these phenomena, consistent and quantitative computational modeling of the microstructure evolution at the oxygen electrode–electrolyte interface is scarce. We combine thermodynamic, 1D kinetic, and 3D phase-field modeling to computationally reproduce the element redistribution, microstructure evolution, and corresponding ohmic loss of this interface. The influences of different ceramic processing techniques for the CGO interlayer, i.e., screen printing and physical laser deposition (PLD), and of different processing and long-term operating parameters are explored, representing a successful case of quantitative computational engineering of the oxygen electrode–electrolyte interface in SOFCs.


INTRODUCTION
Owing to the advantages of high-energy efficiency, broad fuel options, low pollutant emissions, and scalable stacks, solid oxide fuel cells (SOFCs) are attracting increasing attention in the clean and efficient generation of electrical power from renewable sources. However, worldwide commercialization of SOFCs is hindered due to problems induced by the high operating temperature, i.e., high material cost and materials degradation during long-term operation. Within the past two decades, significant progress has been made in the selection of suitable materials of SOFCs to reduce operating temperature from~1000 to 500-800°C, of which the latter is known as the intermediate temperature SOFCs (IT-SOFCs) 1,2 .
Generally speaking, the structure of an SOFC cell consists of a dense electrolyte sandwiched by two porous electrodes (cathode and anode). The electrolyte transports oxygen ions from the cathode, where the oxygen is being reduced into oxygen ions, to the anode. High ionic conductivity is an important prerequisite when selecting the electrolyte materials. 8 mol.% yttria-stabilized zirconia (8YSZ) is the most wildly used electrolyte material so far. Within the intermediate temperature range, gadolinia-doped ceria (CGO) with high ionic conductivity is also regarded as a good electrolyte material for IT-SOFCs, though it may exhibit a minor degree of electronic conductivity at temperatures above 550°C, reducing the SOFC efficiency 3 . The cathode is devoted to reducing oxygen molecules into O 2ions. For IT-SOFCs, the (La 1−x Sr x )(Co 1−y Fe y )O 3-δ (LSCF) perovskite with mixed ionic and electronic conductivity (MIEC), good catalytic properties, and good stability has been widely adopted as the cathode material 4,5 . Unfortunately, LSCF reacts with YSZ, forming insulating phases of La 2 Zr 2 O 7 and SrZrO 3 6 . A reaction barrier layer made from CGO is often required. The anode material should possess high electronic conductivity and good electrocatalytic activity for fuel oxidation and electrons release. The most commonly used anode material for IT-SOFCs is Ni-YSZ cermet since its invention in 1970, and this situation seems to remain until considerable improvement is reached in the future development of oxide-based anode materials 2 . In addition to the outstanding performance of the individual materials mentioned above, good mechanical stability and chemical inertness between each component of IT-SOFCs are stringent requirements to develop high-performing and durable full cells 7 . Considering these, the anode-supported Ni-YSZ (anode)/YSZ (electrolyte)/CGO (barrier)/LSCF (cathode) cells are the state of art IT-SOFCs having promising initial performance in the desired temperature range 5,8 .
The LSCF cathode suffers from performance degradation, which to a large extent is caused by phase transformation or microstructural changes. Apart from the degradation processes in LSCF itself, such as Cr poisoning, sulfur poisoning, coarsening of the microstructure, and loss of conductivity etc. [8][9][10][11][12] , interactions of the YSZ electrolyte with the CGO barrier as well as the LSCF cathode are two important contributors to the degradation [13][14][15][16] .
In an IT-SOFC with the LSCF cathode, the CGO barrier layer is introduced to limit the interdiffusion and reaction between the LSCF cathode and the YSZ electrolyte, which otherwise will lead to formation of insulating phases, e.g., La 2 Zr 2 O 7 (LZO) and SrZrO 3 (SZO). CGO is chosen as the material for the barrier layer, not only due to its high ionic conductivity at intermediate temperature range as mentioned above but also due to its good chemical and thermomechanical compatibility with YSZ and LSCF 13,[17][18][19][20] . The literature reported CGO barrier layer is effective, but not yet sufficient to completely stop the interdiffusion and formation of zirconate phases. The above degradation phenomena still happen in IT-SOFCs with LSCF cathode during either cell fabrication or long-term operation as described further below.
In order to fabricate the YSZ (dense)/CGO (preferably dense)/ LSCF (porous) structure, wet deposition techniques, e.g., screen printing, spray deposition, and dip coating, are commonly employed. Usually, a two-step sintering is performed, during which the interdiffusion and reactions already take place, as illustrated in Fig. 1. In this work, the fuel electrode side is not taken into consideration. On the oxygen electrode side, we divide the relevant cell production and operation into three steps. In Step 1, a CGO layer is deposited on the pre-sintered YSZ electrolyte followed by high-temperatures sintering (1200-1350°C) to achieve an effectively dense diffusion barrier layer of CGO. Such high temperatures can activate the interdiffusion processes between CGO and YSZ, which eventually leads to the formation of a CGO x YSZ (1−x) solid solution with low ionic conductivity, the well-known Kirkendall voids, or even dopant migration 13,14,17,[20][21][22] . In Step 2, a LSCF layer is deposited on the CGO layer, followed by sintering at a slightly lower temperature (950-1150°C), where the interdiffusion across the YSZ-CGO interface will still continue, but in a much less degree due to much slower diffusion kinetics. Despite of the existence of the CGO barrier layer, formation of the SZO phase at the YSZ-CGO interface is confirmed by a large number of experimental investigations so far 13,15,[22][23][24][25][26][27][28][29][30] . Sr accumulated at the YSZ-CGO interface is believed to diffuse from the LSCF cathode through the CGO barrier layer mainly during the cathode-sintering process. During the long-term operation process, i.e., Step 3 in Fig. 1, the kinetics of Sr diffusion (at 650-750°C) is much slower than in Step 2. But the growth of SZO at the YSZ-CGO interface continues and contributes to the cell longterm degradation when the operation time reaches a magnitude of >10 3 h 16,24,28 . A dense CGO barrier layer realized without hightemperature sintering can effectively limit these interdiffusioninduced degradation processes. Pulsed laser deposition (PLD), or physical vapor deposition (PVD), is applicable to produce thin and dense CGO epitaxy on YSZ substrate 16,22,23 . Other methods, such as co-firing of the YSZ and CGO layers at reduced temperature utilizing a sintering aid, provides the possibility of compositional and microstructural engineering 31 . Nevertheless, screen printing is so far the most widely adopted method for producing the CGO barrier layer (and the LSCF cathode as well), considering its high deposition rate and low cost, though it is difficult to achieve a fully dense CGO barrier layer as by PLD or PVD.
The area-specific resistance (ASR) is a common variable to represent the performance of SOFC cells. The ASR can be split into polarization and ohmic resistance. The ohmic resistance (ASR Ω ) can be ascribed to mainly the electrolyte resistance 32 . The interdiffusion between YSZ and CGO and the formation of the CGO x YSZ (1−x) solid solution causes a significant increase in the ohmic resistance. Recent studies indicate that this part of performance degradation depends strongly on the fabrication method of the CGO barrier layer (e.g., screen printing or PLD, etc.) as well as the fabrication conditions (e.g., sintering temperature and time) 16,28,29 . The PLD method avoids the barrier layer sintering step (Step 1) where most of the YSZ-CGO interdiffusion takes place. In addition, it can produce a thin dense CGO layer with a large grain size to slow down the Sr diffusion, resulting in a very little increment in ASR Ω after cathode sintering (Step 2) and after the long-term operation (Step 3) 16 . For a screen-printed CGO barrier layer, achieving high densification often requires a high sintering temperature (>1200°C) in Step 1. This results in significant YSZ-CGO interdiffusion and thus increment of ASR Ω . The resulting YSZ-CGO interdiffusion zone can further influence the Sr diffusion and formation of SZO in Steps 2 and 3 29 . The higher the CGO sintering temperature, the denser the CGO layer, the more the CGO-YSZ interdiffusion, however the less the Sr diffusion and SZO formation. This makes it difficult to identify appropriate fabrication parameters in order to reach an optimum condition to minimize both the YSZ-CGO interdiffusion and the SZO formation. Moreover, it is even harder to derive a general rule by making a crosswise comparison among the experimental results existing in the literature, since the fabrication conditions of the tested cells from different labs are hardly the same. Considering the above, computational numerical modeling of the microstructure evolution in SOFCs shall be performed not only to reproduce the kinetic and microstructural changes observed experimentally but also to correlate with the cell performance, hence be able to provide a reliable prediction of cell performance as functions of fabrication and long-term operation conditions. If successful, this will open a new route of computational engineering of the oxygen electrode-electrolyte interface, targeting both superior initial performance and good long-term durability.
Recently, the current authors have tried modeling the YSZ-CGO interdiffusion and thin layer SZO formation at the YSZ-CGO interface 33,34 . Several simplifications were made in the model, e.g., consideration of no CGO grain growth during high-temperature sintering and treating the SZO phase as a continuous layer instead of dispersed distribution. These simplifications do not reflect what actually happens in IT-SOFCs. To the best of our knowledge, no other computational study has been reported which models the CGO-YSZ interdiffusion, Sr diffusion, and SZO formation during cell production and long-term operation in a systematic way. This is thus the aim of the current work: (i) to develop a set of thermodynamic and diffusion kinetic databases providing an accurate description of the phase relations and elemental diffusivities in the YSZ-CGO-LSCF system based on the calculation of the phase diagram (CALPHAD) approach; (ii) to kinetically model the experimentally observed 16,24,28 YSZ-CGO interdiffusion, Sr diffusion across the CGO layer and formation of SZO by performing both one dimension (1D) solid-state reaction and three-dimension (3D) phase-field simulations; (iii) to identify the correlations between the cell fabrication and operation parameters, the microstructures, and the ohmic resistance according to the current simulation results, and compare with the previous experimental findings 29 .

RESULTS AND DISCUSSION
Simulation of SZO formation in screen-printed CGO Figure 2 plots the simulated composition profiles in comparison with previous experimental data measured at the CGO-YSZ interface ranging from 5 to 11 μm 24 , for the case of a SOFC with screen-printed CGO barrier layer during the different steps in cell production and long-term operation. In Fig. 2a, the distribution of Ce and Zr after the barrier layer sintering step at 1523 K for 2 h is presented, showing a ∼1-µm thick interdiffusion zone at the interface. We have previously reported a systematic investigation on the interdiffusion process between the YSZ electrolyte and the CGO barrier layer, including both experimental investigations and 1D kinetic modeling 33 . The detail will therefore not be repeated here. Figure 2b presents the simulated composition profile after the cathode-sintering step at 1373 K for 2 h compared with the experimental data reported by Kiebach et al. 24 . An accumulation of Sr on the Ce-rich side of the ZrO 2 -CeO 2 interface can be clearly observed from the experiment data and is also predicted by the simulation. Why is Sr not accumulated on the ZrO 2 -rich side can be well explained by investigating the Sr diffusion route on the isothermal section of the ZrO 2 -CeO 2 -SrO system. As can be seen from the phase diagram shown in Supplementary Fig. 1, the ZrO 2 -SrO binary subsystem represents a direct contact between YSZ and LSCF, where SZO forms due to the reaction between them. When a layer of CeO 2 is placed in between, diffusion of Sr in the Gibbs triangle shall follow the path starting from the SrO corner and going through the SrO-CeO 2 subsystem to reach the ZrO 2 -CeO 2 side, as shown in Supplementary Fig. 1. Since there is no intermediate reaction in the SrO-CeO 2 system, Sr will first enter the two-phase region of Cubic + SZO. In order to reach the twophase equilibrium, Sr diffuses from the right boundary of the CeO 2 layer to the Ce-rich part of the interdiffusion zone (IDZ), and further reacts with the ZrO 2 dissolved in CeO 2 (during the barrier sintering step), and forms SZO in the Ce-rich part of IDZ. As the amount of Sr on the Ce-rich side of IDZ increases, the diffusion front of Sr will finally reach the equilibrium between the miscibility gap of Cubic ZrO 2 + CeO 2 and SZO, indicating growth of SZO into the YSZ region. As can be seen in Fig. 2c, when the SOFC cell is tested at 973 K for 2000 h, the diffusion front reaches the original interface between YSZ and CGO. So for the screen-printed CGO barrier layer, the property of CGO and more specifically the IDZ layer which can be controlled by varying the processing conditions, is of great importance with regard to stop or slow down the Sr diffusion and the SZO formation.
Formation of SZO is the result of Sr diffusion across the CGO barrier layer to the Ce-rich side of IDZ and reaction with ZrO 2 dissolved in CeO 2 , and the process continues during long-term operation at intermediate temperature ( Step 3). Figure 3 shows the simulated amount of formed SZO along the diffusion distance for different operation time at 973 K in the case of screen-printed CGO barrier layer, compared with the systematical investigation by Rinaldi et al. 28 . The numerical modeling was performed considering sintering at 1300°C for 3 h for the 6-µm thick CGO barrier layer screen-printed on an 8-µm YSZ substrate, followed by screen printing and sintering of the LSCF layer at 1100°C for 3 h. The long-term operation was performed at 700°C. The simulated changes in the volume percentage of SZO during the long-term test is recorded at 0, 1900, 4700, and 10,700 h, achieving a reasonable agreement with the experimental data based upon the segmentation of the low kV ESB and InLens grayscale electrons images and complemented with 3D EDX elemental mapping for the 10,700 h sample. It can be seen that the peak value (i.e., the maximum volume percentage of formed SZO) increases from 3% at 0 h (after cell production) to 25% at 10,700 h (after 10,700 h operation at 973 K). The peak position within the IDZ region moves towards the Zr-rich side, which is consistent with the current finding of the Sr diffusion path. A wider distribution of SZO is predicted based on the current simulation, implying a deeper expansion of SZO into the CGO layer. In addition to the simplifications adopted in current modeling, the discard of elastic strain energy could be one reason for the suppression of SZO expansion in the real case. Our simulation result further demonstrates the applicability of current modeling to simulate the formation of SZO in IT-SOFC in the case of the screen-printed CGO barrier layer. Figure 4 plots the simulated SZO formation at the YSZ-CGO interface in the case of PLD CGO, after cathode sintering or during long-term operation. In Fig. 4a, b, formation of SZO is predicted at the CeO 2 -ZrO 2 interface, which is of a reasonably small amount in Fig. 2 Composition profiles across the oxygen electrode-electrolyte interface (screen-printed CGO). Simulated composition profiles (dotted lines) in comparison with previous experimental data (scattered symbols) 24 , for the case of a SOFC with screen-printed CGO barrier layer during the different steps in (a, b) cell production and (c) long-term operation. The distance at 0 μm corresponds to the left boundary of the ZrO 2 layer, and the original ZrO 2 -CeO 2 interface is located at 8 μm. agreement with the experimental data measured at the CGO-YSZ interface ranging from 5 to 9 μm. Comparing with the case of the screen-printed CGO barrier layer, i.e., Figure 2, the PLD CGO layer shows a much better effect on hindering the Sr diffusion and the SZO formation. This can also be seen in Fig. 4c that the accumulation of Sr within the simulated IDZ is very slow (<0.3% after 2 h sintering at 1373 K and <0.4% after 1500 h operation at 973 K).

Simulation of the formation of SZO in PLD CGO
The dense CGO barrier layer deposited by the low-temperature PLD technique eliminates CGO-YSZ intermixing due to no hightemperature sintering, and also prevents extensive SZO formation. Only the grain boundaries in the epitaxial CGO barrier layer facilitate Sr diffusion to the YSZ electrolyte to form SZO 16 . It is also reported that the SZO grain generated at the quadruple junction of the YSZ and epitaxial CGO grains represents an elongated spreading shape along the CGO-YSZ interface 16 . Figure 5 presents the 3D phase-field simulation of the microstructure evolution during SZO formation at 1373 K for 2 h. The simulation results indicate the random formation of SZO at the quadruple junction begins at around 1200 s, during which Sr should diffuse across the CGO barrier layer. The particles merge to form a contiguous region after growth for 1 h, and then grow into an elongated spreading shape along the CGO-YSZ interface due to faster grain boundary diffusion of Sr 2+ . The distribution of Sr is scanned along two lines of L1 and L2, where the Sr content differences in the bulk and at grain boundaries is~1.2% and is much lower than that in the regions of SZO and SrO. The simulation result reproduces the experimental findings reported in 16 reasonably well. In Fig. 5, the amount of Sr accumulates within the IDZ is calculated and compared with the result of 1D DICTRA simulation, showing good agreement. These results further prove the current modeling of the SZO formation is also valid in the case of the PLD CGO barrier layer.

Evolution of the ASR ohm
Our current model can also be used to estimate the electrolyte resistance after cell production and during cell long-term operation by using Eq. (9). The previous experimental study 29 demonstrates that the CGO sintering temperature strongly influences the microstructure evolution in the YSZ/CGO/LSZF region. It was shown that the SZO formation is largely reduced by a screen-printed CGO diffusion barrier layer sintered at sufficiently high temperature (1400°C), forming a dense IDZ at the YSZ-CGO interface, resulting in low ohmic resistance after cathode sintering 29 . The current work tries to simulate these processes as well as to further study the effect of the barrier layer sintering condition. The numerical modeling is chosen to be identical to Fig. 2, and the simulation conditions are selected to keep consistent with the experiment 29 , where a 5.5-μm thick screenprinted CGO was sintered at 1100, 1250, and 1400°C for 3 h on a 200-μm thick YSZ substrate. The subsequently screen-printed LSCF layer (on top of the CGO layer) was sintered at 1100°C for 3 h, and the cell long-term operation was performed at 800°C. In this work, DC1, DC2, and DC3 refer to the samples where the CGO layer is sintered at 1100, 1250, and 1400°C, respectively. Figure 6 gives the calculated ASR ohm in DC1, DC2, and DC3 at the three stages of barrier layer sintering (Step 1), cathode sintering (Step 2), and long-term operation (Step 3). As for Step 1, compared with the theoretical ASR ohm (433 mΩ cm 2 ) estimated when no IDZ or SZO forms, DC1 and DC2 show only a slight increase in ASR ohm , i.e., 435 mΩ cm 2 for DC1 with a barrier layer sintering temperature of 1373 K and 438 mΩ cm 2 for DC2 with a barrier layer sintering temperature of 1523 K, indicating weak interdiffusion between YSZ and CGO at lower sintering temperature. The ASR ohm of DC3 sintered at 1673 K reaches 461 mΩ cm 2 after the 3-h barrier layer sintering stage. Figure 7a gives the simulated IDZ thickness for DC1, DC2, and DC3 after the barrier layer sintering at 1100, 1250, and 1400°C, respectively. It is clear that the IDZ thickness increases with the sintering time, and the growth rate increases with the sintering temperature. Higher temperature for barrier layer sintering leads to faster interdiffusion between CGO and YSZ, and hence results in thicker IDZ and higher ASR ohm after Step 1.
After Step 2, DC1 then exhibits a dramatic increment in ASR ohm : from 430 to 483 mΩ cm 2 , which is close to the experimentally measured value of 481 mΩ cm 2 (see ref. 29 ). This increment in ASR ohm is mainly due to the formation of SZO at the CGO-YSZ interface. The amount of Sr accumulated in the IDZ is calculated and presented also in Fig. 7. At this stage, it can be seen that DC1 shows the fastest increment of Sr content, while DC2 accumulates slightly more Sr than DC3, implying that DC3 has the best ability to resist SZO formation. This is consistent with the predicted ASR ohm of DC2 (442 mΩ cm 2 , increased by 7 mΩ cm 2 ) and of DC3 (461 mΩ cm 2 , nearly unchanged) in Step 2. The current modeling results show a reasonable trend in the case of developing a universal computational model, where all the initial conditions and processing parameters are taken into account. The deviation from the experimentally measured results of DC2 (485 mΩ cm 2 ) and DC3 (431 mΩ cm 2 ) 29 can be partly attributed to the simplifications adopted in our model, where no microstructure information (porosity, size, and distribution of SZO, etc.) is taken into account. In addition, the calculation of ASR from Eq. (9), which considers the parallel distribution of cubic solid solution and SZO to calculate the effective resistivity, could be another source of the ASR deviation. Last, but not least, the lack of accurate experimental detail in the published literature (due to  16 . c Estimated amount of Sr in the simulated IDZ and the corresponding ASR ohm at 700°C. Distance at 0 μm corresponds to the left boundary of the ZrO 2 layer, and the original ZrO 2 -CeO 2 interface is located at 8 μm. confidentiality and IPR) is another important factor contributing to the deviation from the experimental results.
During long-term operation ( Step 3) at 1073 K, DC1 exhibits again the fastest SZO formation. The resulting ASR ohm increases from 483 to 530 mΩ cm 2 after 10,000 h. Apparently, DC1 with CGO barrier layer sintered at 1100°C shows the worst ability to resist the formation of the SZO phase in both Steps 2 and 3. Although DC2 shows slightly faster Sr accumulation than DC3 at the IDZ, it turns out to show the lowest ASR ohm after 10,000 h (464 mΩ cm 2 ), since the overall ASR ohm increment in Steps 1 and 2 is the lowest among the three samples. DC3 shows the least Sr diffusion, SZO formation, and thusly the lowest ASR ohm increment from 461 to only 469 mΩ•cm 2 after 10,000 h, indicating the best performance of preventing ohmic loss caused by SZO formation. The major degradation procedure of DC3 should be the very thick IDZ generated in Step 1.
For all the three samples of DC1-3, the simulation conditions in Steps 2 and 3 are identical. The different behavior of Sr diffusion and accumulation in IDZ and ASR ohm increase is caused by the different barrier layer sintering temperatures applied in Step 1. Considering that the difference in the diffusion rate of Sr in DC1-3 is caused by the different microstructure at the IDZ after Step 1, the grain size can be selected as a microstructure indicator. Our calculated grain size of CGO after Step 1 is shown in Fig. 7b, representing an increasing trend from 0.15 μm at 1100°C to 0.63 μm at 1400°C, in good agreement with the experiment values 29 . In Fig. 7b, the effective interdiffusion coefficient of Sr in CeO 2 after the barrier sintering (Step 1) is estimated using Eq. (4) at both 1100°C (the cathode-sintering temperature) and 800°C (the long-term operation temperature). It can be clearly seen that the Sr diffusion slows down by a factor of about 10 5 when temperature decreases from 1100 to 800°C, and the degree of interdiffusion decreases when the barrier sintering temperature (T Step1 ) increase. That is to say that the higher sintering temperature in Step 1 can result in larger grain size in IDZ, better prevention of Sr diffusion, and thus less SZO formation.
However, based on the current simulation result, the thicker IDZ shown in Fig. 7a generated at higher sintering temperature can lead to a significant increase of ASR ohm , e.g., DC3. So both of the two variables, i.e., the thickness of IDZ and diffusion of Sr in CGO, shown in Fig. 7a, b have a complex influence on the formation of SZO and the resulting ASR ohm increase, which should be carefully controlled to reach optimum cell performance. In the case of PLD CGO barrier layer, the absence of Step 1 avoids the formation of  thick IDZ, and the large CGO grain size of~2 μm results in about ten times slower Sr diffusion than that in screen-printed CGO, as shown in Fig. 7b. Together with the thin PLD CGO layer, the current simulation gives an ASR ohm below 100 mΩ cm 2 , as shown in Fig. 4c. In spite of the advantages of the PLD CGO barrier layer, the screen-printed CGO layer remains the most economic and promising method for the purpose of commercialization. In such case, the current work gives a technical route of computational aided optimization of cell fabrication conditions for the oxygen electrode-electrolyte interface and also a model to predict cell degradation in solid oxide fuel cells.
In summary, models are developed to reproduce the YSZ-CGO interdiffusion and the SZO formation at the YSZ-CGO interface in this work. The YSZ-CGO-LSCF system is simplified to the ZrO 2 -CeO 2 -SrO system, accounting for interdiffusion of the major elements (Zr and Ce) and the transportation of Sr from LSCF and reaction with YSZ forming SZO. The CALPHAD-type thermodynamic and diffusion kinetic information of the ZrO 2 -CeO 2 -SrO system is carefully evaluated and used as unique input to all the numerical modeling. Our 1D and 3D models reproduce the experimentally observed element distribution and microstructure evolution during the relevant steps in cell production and operation, for both screen-printed and PLD CGO barrier layers, indicating the credibility of our model. The calculated ASR ohm suggests a trackable performance evolution within the oxygen electrode-electrolyte interface during the entire fabrication and long-term operation of SOFCs with the help of computational materials modeling. The two key factors of the cell degradation, i.e., IDZ caused by CGO-YSZ interdiffusion and Sr diffusion induced SZO formation, are analyzed to be strongly dependent on the fabrication condition of CGO barrier layer and a trade-off should be made when choosing the optimum processing route. It is expected the current route of computer engineering of the oxygen electrode-electrolyte interface in solid oxide fuel cells can serve as quantitative guide to the selection of fabrication conditions and prediction of cell degradation.

Construction of thermodynamic and diffusion kinetic database
One set of CALPHAD-type thermodynamic and diffusion kinetic databases is constructed and applied as the accurate numerical input for all the simulations performed in this work. The basic idea of CALPHAD approach is to describe the thermodynamic properties of each phase in a system with a mathematical model containing adjustable parameters, which can be evaluated by the optimization method to fit the model to all available experimental information. It is then possible to calculate phase diagrams as well as thermodynamic properties of all the phases 35 . The compound energy formalism (CEF) 36 is widely used in CALPHAD modeling, which is introduced here to describe the Gibbs energy for the phases in the current YSZ-CGO-LSCF system. The general formula of a Gibbs free energy function 0 GðTÞ ¼ G T ð Þ À H SER (lattice stability) for stoichiometric compounds can be expressed as, where H SER is the weighted molar enthalpy of the stable element reference (SER), the pure element in its stable state at 298.15 K, and T the absolute temperature. The subscript i represents the phase described and A i~Fi are parameters to be optimized. There are nine elements included in the current YSZ-CGO-LSCF system, i.e., Zr-Y-Ce-Gd-La-Sr-Co-Fe-O, and it is an extensive project to establish the thermodynamic database for the entire nine-element system, which is also unnecessary for the current simulation.
In the current work, the YSZ-CGO-LSCF system is simplified to a model system of ZrO 2 -CeO 2 -SrO. Specifically, the interdiffusion between YSZ and CGO can be simplified to the generation of interdiffusion zone (IDZ) between cubic ZrO 2 and CeO 2 , where the two boundaries of IDZ are defined when detecting 0.1 at. % ZrO 2 on the CeO 2 -rich side or 0.1 at. % CeO 2 dissolved in the ZrO 2 part. The cubic ZrO 2 and CeO 2 phases can be thermodynamically treated as a single Cubic phase following our previous simulation result 33,34 . We then use SrO to represent the LSCF perovskite that provides Sr 2+ (or Sr) to react with ZrO 2 in IDZ. Such a simplification can reasonably and sufficiently account for the main purpose of studying the interdiffusion between YSZ and CGO, the diffusion of Sr 2+ (or Sr) through the CGO layer, and the formation of SZO at the YSZ-CGO interface. The end-members of the current ZrO 2 -CeO 2 -SrO system will then be used to represent the corresponding component in the YSZ-CGO-LSCF system in the kinetic modeling. The lattice stabilities of the Cubic, SrO, and SZO phases are taken from previous thermodynamic modeling [37][38][39] , and are carefully evaluated to construct the current thermodynamic database as described in Supplementary Methods. No higher-order parameter is used except for the Cubic phase, in order to describe their mutual solid solubility. The molar Gibbs energy of the cubic phase can be described as, where R is the gas constant, 0 L Cubic is the interaction parameter, c Zr 4þ ð Þ O 2À ð Þ 2 and c Ce 4þ ð Þ O 2À ð Þ 2 are the mole fraction of cubic ZrO 2 and CeO 2 , respectively. The thermodynamic description of the ZrO 2 -CeO 2 -SrO system is listed in Table 1, where the sublattice model with charged species is adopted. The diffusion of Ce and Zr between CGO and YSZ can be treated as in our previous work 33,34 . On the other hand, the Sr migration through CGO layer is rather complex, and different mechanisms have been proposed in the literature, including gas-phase diffusion through the pore phase, surface diffusion, or bulk/grain boundary diffusion 40,41 . For the gas-phase diffusion, it is generally believed that SrO derived from the strontium segregation in LSCF reacts with water vapor to form the volatile Sr(OH) 2 gas species, and afterward transports in the vapor form and deposits again at the reactive sites as SrO. The fast surface diffusion indicates the transport of Sr 2+ (or Sr) along the crack wall or the inner surface in the CGO barrier, while the bulk/grain boundary diffusion driven by the chemical potential of Sr 2+ (or Sr) appears to be slower, especially for the bulk diffusion which is significantly slower than the other diffusion mechanisms. These different diffusion mechanisms play different roles, depending on the experimental conditions, such as porosity of the CGO barrier layer, temperature, P(O 2 ), and P(H 2 O), etc. However, there is still a lack of experimental proof and quantification in the literature on how the different diffusion mechanisms contribute to the overall Sr diffusion. In this work, instead of quantitatively distinguishing the contributions from different diffusion mechanisms or the form of Sr to transport (as gas species, neutral atoms, or cations) which is more or less impossible based on available experimental information, we developed a quantitative description of the overall Sr 2+ migration and SZO formation by using an "apparent" diffusion coefficient. In terms of the diffusion kinetic data, the Arrhenius expression is sufficient to describe the temperature dependence of diffusivity of each charged species in different phases, i.e., where D 0 is the prefactor and Q is the activation energy. The elemental diffusivity data of the current ZrO 2 -CeO 2 -SrO system evaluated using previous experimental data [42][43][44][45] are listed in Table 2, presented as the CALPHAD-type atomic mobility functions as input data for the current 1D simulation. The methodologies to evaluate atomic mobility and the quality of resulting parameters are detailed in Supplementary Methods, Supplementary Fig. 2a-d  In the current database, the "effective" diffusivity formulated as Hart's equation 46 is used, which consists of two contributing items corresponding to the bulk and grain boundary (gb) diffusion, respectively. As shown in Eq. (4), the grain boundary item is used to account for the faster Sr migration via grain boundary, surface, and vapor diffusion without further distinguishing between them, due to lack of precise experimental information.
where δ and d are the grain boundary thickness and the grain size, respectively. As for δ, a nominal value of 5-10 atomic diameter is suggested as in ref. 47 , and δ is thus assumed to be 1 nm in the current simulation. As for d, only the grain growth in the CGO layer is included in the present database due to the model of 1D simulation, which will be described in detail in the modeling part. The sintering temperature for the YSZ electrolyte layer is in general a few hundred degrees higher than that of CGO. We hence assume no grain growth in the YSZ layer during sintering of barrier layer and cathode and during long-term operation (Steps 1-3 in Fig. 1). Previous studies [48][49][50] performed systematic experimental investigations on the grain growth behavior of the CGO system as functions of temperature and time, and proposed the following analytical equation, where d is the average grain size at time t, d 0 is the initial grain size, K 0 is a pre-exponential constant, R is the gas constant, T is the absolute temperature, and Q g denotes the activation energy for grain growth. According to the previous work 50 , the parameters in our simulation are set to be d 0 = 1.5 × 10 −7 m, K 0 = 0.04 m 2 s −1 and Q g = 477 kJ mol −1 , respectively.

1D kinetic modeling
As the current YSZ-CGO-LSCF system is simplified to the ZrO 2 -CeO 2 -SrO model system, the phase relations can be computed based on the current thermodynamic description listed in Table 1 51 . Figure 8a shows the calculated isothermal section of the ZrO 2 -CeO 2 -SrO system at 1523 K. The SZO phase is located in the middle of the ZrO 2 -SrO binary line, while a miscibility gap between cubic ZrO 2 and CeO 2 can be found on the ZrO 2 -CeO 2 side. This miscibility gap of the phase Cubic extends to the ZrO 2 -SrO side to reach a two-phase equilibrium with the SZO phase.
We first construct the 1D numerical model to simulate the elemental interdiffusion and formation of SZO in the ZrO 2 -CeO 2 -SrO system, considering different methods for fabricating the barrier layers, i.e., the screen-printed CGO and PLD CGO. Figure 8b gives the 1D numerical model in the case of the screen-printed CGO barrier layer. The three-step (Steps 1-3) simulation is in correspondence to the three processes described in Fig. 1. In order to guarantee the comparability of the current simulation results with previous experimental data, the simulation conditions are set according to the experimental work 24 . In Step 1, the initial thickness of ZrO 2 and CeO 2 layer is set to be 8 and 6 μm, respectively, representing the initial YSZ electrolyte and the CGO barrier layer. The sintering temperature and time are set to be 1523 K and 2 h. The grain growth in CeO 2 region should be considered due to the dramatic increase in grain size during barrier layer sintering at high temperature 24,[48][49][50] , which has a strong influence on the elemental diffusion behavior 29 . The reason for no grain growth considered in the ZrO 2 layer is that the YSZ electrolyte layer has been pre-sintered at~1673 K 28 , resulting in quite a large grain size. Co-sintering together the CGO barrier layer at 1523 K will thus not lead to significant grain growth in the YSZ layer. Grain size of 2 μm is adopted for the ZrO 2 region in our model. At the right boundary of the CeO 2 region (Fig. 8b), the appearance of LSCF cathode is represented by a virtual SrO layer with no thickness and one boundary condition on the activity of SrO phase (a(SrO)), which is set to be 0 in Step 1. In Step 2, i.e., the cathode-sintering process, the a(SrO) is then set to be 1 and the simulation condition is set to be at 1373 K for 2 h. The grain growth in the CeO 2 region should still be considered in this step. During long-term operation, i.e., Step 3, the grain growth in CeO 2 region is ignored due to the low operating temperature at 973 K. The simulation time is selected to be 2000 h to be the same as reported in ref. 24 . Figure 8c gives the 1D numerical model in the case of the PLD CGO barrier layer. The simulation conditions are the same as in the reported experimental work 16 . A CeO 2 barrier layer (representing CGO) with a thickness of 0.6 μm is placed on the right side of the ZrO 2 region (8 μm, representing YSZ). Since the PLD CGO barrier layer exhibits an epitaxy microstructure, the grain size of CeO 2 is set to be identical to that of ZrO 2 i.e., 2 μm. The deposition of the PLD CGO layer requires no further high- temperature heat treatment (sintering) afterward. In this case, only two processing steps are considered in our simulation, i.e., Step 2-Cathode sintering and Step 3-Long-term operation. The simulation conditions are set to be 1373 K for 2 h and 973 K for 1500 h, for Steps 2 and 3, respectively. Similar to the case of screen-printed CGO, the CGO grain growth is considered during the cathode-sintering step, and the boundary condition of a(SrO) = 1 is applied to represent the presence of the LSCF cathode. The current 1D numerical model is also used to reproduce the change of volume percentage of SZO in IDZ during long-term operation for up to 10,700 h determined by focused ion beam-scanning electron microscopy (FIB-SEM) serial sectioning (Fig. 3) 28 , as well as to estimate the effect of barrier layer sintering temperature on the ohmic loss to compare with the previously reported data (Figs. 6 and 7) 29 . In these cases, the modeling parameters (i.e., thickness of YSZ and CGO, temperature, time, etc.) are adjusted according to different experimental conditions adopted in the previous work 28,29 . In this work, the above numerical kinetic modeling is carried out using the DICTRA software package implemented in Thermo-Calc 51 . The homogenization model 52 in DICTRA is activated to reveal a good coupling between thermodynamic and diffusion kinetic information in Tables 1 and  2. The non-uniform grid is sufficiently refined close to the ZrO 2 -CeO 2 interface region and the right boundary of the CeO 2 region to guarantee numerical stability during simulation, as illustrated in Fig. 8b, c. It should be further pointed out that due to the simplification of using an apparent diffusion coefficient to represent the Sr migration, i.e., Eq. (4), in this paper, we refer Sr diffusion all in the form of Sr 2+ irrespective of the different diffusion mechanisms. The different "effective" diffusion behavior of Sr 2+ for the 1D cases of screen-printed CGO and PLD CGO is reproduced by considering the evolution of grain size, with no further consideration on the change of porosity. Meanwhile, the current model only includes the effect of temperature and time in the simulation of Sr 2+ migration and SZO formation during the 1D long-term operation stage. The effect of electric boundary conditions on the local P(O 2 ) at the oxygen electrode-electrolyte interface is not considered in the current frame of calculation, since it would lead to the complex phenomenon of simultaneous anion and cation migration induced by both electric and chemical potential gradient, with thermodynamic and kinetic phase transformation due to element redistribution 53 . This is also based on considerations that in reality, the temperature has an even bigger effect than P(O 2 ) simply due to fast diffusion kinetics at elevated temperatures. These simplifications adopted in current modeling may cause deviation from experimental data, and continuous efforts are being made to include more variables of real conditions in our future computational engineering studies.

3D phase-field modeling
Apart from 1D simulation, the current work takes advantage of the multiphase field model 54,55 to simulate microstructure evolution and SZO formation in the case of the PLD CGO barrier layer. Comparing with the 1D DICTRA simulation, which is suitable for the simulation of SZO formation on a very large time scale, the phase-field simulation mainly focuses on reproducing the details of microstructure evolution within a short period of time. Since the phase-field simulation is computationally expensive, in this work, it is employed only to describe the SZO formation and growth during Step 2 in the case of the PLD CGO layer. In the multiphase field theory, each phase grain α is distinct from the other phases by its individual phase field ϕ α . Three phases are considered in the current ZrO 2 -CeO 2 -SrO model system, i.e., the Cubic, SrO, and SZO phases. The sum constrain exists as P α¼1 3 ϕ α ¼ 1. The governing equation for each phase-field can be  Due to lack of experimental data, the grain boundary diffusivity is assumed to be 10 5 times higher than the bulk diffusivity.
b Diffusivity of Sr 2+ in YSZ is assumed to be the same as in CGO due to lack of experimental data. c Diffusivity of Ce 4+ in CGO is assumed to be the same as in YSZ due to lack of experimental data.
K. Cheng et al.
expressed as, where μ αβ is the interfacial mobility, σ αβ is the interfacial energy, Δg αβ is the local deviation from thermodynamic equilibrium and η is the interface thickness. Evolution of the conserved composition field c α is coupled to phase-field and is governed by, Figure 8d shows the computational domain in our SZO formation model using the phase-field method. The cathode-sintering step (Step 2) in the case of the PLD CGO barrier is selected here to study the typical microstructure evolution during formation of SZO phase at the quadruple junction. Considering the grain size of the YSZ phase and the thickness of the CGO layer in the experimental work 16 , a cuboid region of 2 × 1 × 0.5 μm is modeled including three layers of ZrO 2 (1 μm), CeO 2 (0.6 μm), and SrO (0.4 μm) from the left to the right side. Grain boundaries are defined in the middle of the ZrO 2 and CeO 2 region forming a quadruple junction to keep consistent with the experimentally observed CGO epitaxy 16 . The current 3D modeling is carried out using MICRESS (http://www.micress.de.), and the grain boundary diffusion model is activated to detect the faster Sr diffusion along CGO grain boundary and the formation of SZO at the quadruple junction. Under such circumstances, the total diffusivity D α in Eq. (7) should contain an extra flux contributed by grain boundary diffusion, i.e., where D bulk α and D gb α are bulk and grain boundary diffusivities with the prefactor and activation energy shown in Table 2. It should be noted that D gb α contributes to D tot α only at the phase interface since it is also proportional to ϕ α ϕ β and will cancel out in the bulk region. In order to allow for the formation of the SZO phase, the seed density model 56 is applied by placing randomly distributed inactive SZO particles in the simulation area with an initial grain radius of 0. The SZO particles at the quadruple junction can be stabilized when a sufficient amount of Sr is accumulated to reach the two-phase equilibrium with the Cubic phase as shown in Fig. 8a.
In this work, Δg αβ and D α can be directly obtained from the thermodynamic and diffusion kinetic database developed for the ZrO 2 -CeO 2 -SrO system, as shown in Tables 1 and 2. Integration of these currently developed thermodynamic and diffusion databases into the 3D phase-field model can be revealed owning to the good coupling between Themo-calc and MICRESS 51,56 . For simplification, all the phases/grains were assumed to be isotropic. The interface mobility μ αβ for the CGO and YSZ grain boundaries are acquired from literature to be 1.7 × 10 −8 and 6.5 × 10 −10 cm 4 J −1 s −1 , respectively 57,58 . And the interface energy σ αβ for both of the CGO and YSZ grain boundaries are set to 3.0 × 10 −5 J cm −2 as reported in literature 57,58 . The interface thickness η is set to be 3 × 10 −8 m. For all the other interfaces with no experimental data available, μ αβ and σ αβ are adjusted to be 1.0 × 10 −10 cm 4 J −1 s −1 and 2.0 × 10 −5 J cm −2 to simply reach diffusion control of the phase transformation.

Prediction of ohmic loss (ASR ohm )
The ohmic loss at the oxygen electrode-electrolyte interface of an SOFC cell can be ascribed to mainly the electrolyte resistance. In our simulation, an additional contribution to the ASR ohm comes from both the YSZ-CGO IDZ and the SZO ion insulator. We consider an YSZ-CGO electrolyte region in 1D with a thickness of l, which is further split into n small grids. The theoretical ASR ohm can be computed by simply summing up the resistivity of YSZ and CGO layer multiplied by each layer width. When the formation of IDZ and SZO is considered, the phase (volume) fraction and electrical resistance of cubic solid solution and SZO should be first computed separately within each grid. The electrical conductivity of a YSZ-CGO solid solution depending on composition and temperature can be found in previous experimental work 59 , and the data at 700 and 800°C fitted using polynomial equations are adopted in this work. The temperaturedependent conductivity of SZO has also been reported 60 . Then, we consider a parallel distributed cubic solid solution and SZO to calculate the effective resistivity R(x) of each grid. The ohmic loss of the entire bi-layer electrolyte region is estimated by integrating R(x) along the thickness as follow, where σ Cubic and σ SZO are the electrical conductivity of cubic solid solution and SZO, respectively. N Cubic and N SZO are the phase fraction of cubic solid solution and SZO within each grid, which can be calculated based on the currently simulated composition profiles. In this work, the Ohmic loss is calculated based solely on the results of phase fraction obtained from Dictra 1D simulations of Steps 1-3 for both cases of the screen-printed and PLD CGO barrier layers.