Percolation modeling investigation of TPB formation in a solid oxide fuel cell electrode–electrolyte interface

A Monte Carlo percolation model has been developed and utilized to characterize the factors controlling triple phase boundary (TPB) formation in an SOFC electrode. The model accounts for (1) electronic conductor, ionic conductor, and gas phase percolation, (2) competition between percolation of gas and electronically conducting phases, and (3) determination of continuous, though not necessarily fully percolating, paths from TPBs to the bulk phases. The model results show that physical processes near the TPB, such as sorbate transport, signiﬁcantly affect TPB formation in a composite electrode. Active TPB formation is found to be most signiﬁcantly dependent upon continuous and competing percolation of multiple phases. Simultaneously requiring continuous paths and accounting for non-continuous boundary conditions results in lower active TPB formation levels (up to 8% of possible sites) than presented in the literature (75% of possible sites). In addition, the varying ratio of active to potential TPB sites predicted by the current model (up to 80%) differs signiﬁcantly from the constant reported in the literature (80%), which lacks analyses of three-phase percolation, gas phase paths, and gas/current collector boundary conditions. This dependence of active TPB formation on percolation of all three phases is important to understand as a basis for determining SOFC performance and optimization.


Introduction
Fuel cells for portable and stationary power generation have garnered significant interest in recent years largely due to their potential for particularly high energy conversion efficiency and low emissions characteristics. The high efficiency hinges upon relatively controllable and often predictable redox reactions, which occur at the triple phase boundaries (TPBs) within the cell. These TPBs are characterized by the adjacency and electrochemical interaction of electron conducting, ion conducting, gaseous, and (if necessary) separate catalyst phases in each of the electrode-electrolyte interfaces. Thus, in order to achieve appreciable functionality, power density, and high conversion efficiency, one must be able to construct cells with a microstructure that imparts a high volumetric density of these TPBs, ideally engineering the placement of elements of each phase with high accuracy. In addition, percolation of the ionic and electronic phases to the bulk electron and ion conductors at either end of the interface is required to make a TPB electrochemically active.
However, during the operation of a fuel cell, these TPBs at the electrode-electrolyte interface do not all equivalently contribute to the conduction of these species. Long conduction paths and poor interconnectivity between the electrode and the current collector can lead to TPBs which, though active, are not significant contributors to the overall current developed. In addition, some sites may be limited by morphology or electrocatalytic inadequacies that contribute to activation losses. Therefore, the conductivities of the various phases in the electrode structure are important in determining electrode performance, especially for phases with low conductivities (typical of ionic conductors). These key conductivity issues, however, fundamentally depend upon the formation of the TPB and the percolation of the respective phases since these conducting physics require a physical location for the electrochemical reactions to occur. Thus, this work focuses on characterization of the probabilistic physics of the formation of active TPB sites within the electrode. The extension to electrochemical performance and exploration of the interplay between these two key features of electrode performance can then be built upon the knowledge of this work in a future investigation.
Understanding how to achieve high TPB density with good interconnectivity to bulk phases is not a simple matter, and is dependent upon many parameters of the design, materials, and manufacturing methods of the fuel cell. In the case of a solid oxide fuel cell (SOFC), the ionic and electronic phases are of solid-state design, and operational temperatures are typically high enough to eliminate the need for a separate electrocatalyst phase. The manufacturing methods most applicable to this fuel cell include ceramic materials synthesis, particle manufacturing and sizing, mixing of particles of various phases, ball milling, calcining, and sintering, intermixed with processes such as screen printing, tape casting or calendaring, magnetron sputtering, flame and plasma spraying, and aerosol-assisted vapor deposition [1][2][3][4][5][6][7][8]. Although the methods vary widely, the production of electrode-electrolyte interfaces for all of these methods share a degree of randomness, whether it is due to mixing of constituent particles, processing materials to form powders of non-uniform particle size, or simply depositing material without spatial control.
Investigators have attempted to identify those processes of manufacturing that have a significant impact on cell performance. Studies of strontium-doped lanthanum manganite (LSM)/yttria-stabilized zirconia (YSZ) cathode-electrolyte interfaces have identified multiple pertinent aspects. Of primary importance are the particle size distribution [9][10][11] (of which a larger variation in size is preferable [10,11]), thickness of the electrode [9,12], and the temperature at which calcining and sintering steps are carried out [9,12]. Processing temperatures have also been identified as pertinent to the strontium-and iron-doped lanthanum cobaltite (LSCF)/gadolinia-doped ceria (CGO) interface performance [13]. Interestingly, the effects of each of these parameters can be altered through changes in the other parameters; for instance, altering the calcining temperature can result in differing optimal particle size distributions, as revealed by a comparison of [9,10]. Therefore, complex interproperty and interface manufacturing parameter correlations can be expected when investigating fuel cell performance. In addition, the inclusion of mixed electronic-ionic conductors (MEIC) in the structure of an electrode has been shown to potentially improve cell performance, indicating a dependence on the intrinsic material properties themselves [13]. Furthermore, a recent visualization study of a Ni-YSZ anode demonstrated a clear dependence of TPB length on the degree of intermixing of the various materials in the anode structure [14]. It was easily observed that such parameters as tortuosity, adhesion, and degree of interconnection, all of which are difficult to control and characterize in an experimental system, are important features in the formation of active TPBs.
Thus, the premise that the microstructure of the cell, which is itself intimately related to the manufacturing process, critically affects the performance of the cell has been developed and shown in the literature. However, since all current manufacturing processes are stochastic in nature, an effective theoretical under-standing of the relationships and the extents of their impacts compliments any empirical investigations. Thus, researchers have looked to modeling electrode-electrolyte interfaces in order to understand these relationships. Through the use of Monte Carlo methods and a variety of model formulations, certain relationships from the empirical investigations have been expanded upon and new relationships discovered. Particle size and electrode thickness have been confirmed through theory to be of particular importance [15][16][17][18][19], as has the relative amount of electronic and ionic conducting particles in the electrode [15][16][17]19,20]. Through these investigations, however, it has been noted that the optimal volume fractions of phases also depends on the conductivities, reinforcing the inter-property dependence of performance observed in the empirical investigations.
In general, there has also been a demonstrated correlation between the optimal volume fraction of either conductor and its percolation limit, thereby accentuating the applicability of percolation theory to studying the random nature of the fuel cell microstructure [15][16][17]. Novel observations that are more easily observable through models include the location of minimum overpotential in the interface, which tends to be present towards the center of the interface thickness with a preference for the electrolyte side [15], a dependence of optimal electrode composition on the ratio of electronic and ionic conductivities [16], improvement in performance with grading of various interface parameters such as particle diameter and conductivity [17,19], and the result that a given TPB length can be approximated as three times the radius of a particle [20]. These results demonstrate the advantage of modeling the system, as more complex interactions and more detailed phenomena are observable through modeling than with experimentation, providing much-needed insight into the pertinent physics and potential performance improvements.
In spite of the successes of previous models to discover new insights regarding the performance of fuel cells as dictated by interfacial microstructure and properties, there are still certain considerations that are not yet fully characterized or understood. One of the major gaps in understanding thus far is the effect of macropores and the structure of the gaseous delivery and removal pathways on performance. In the construction of an SOFC, pore formers melt out of the structure during heat treatment steps, leaving large pores, but the structure and dispersion of these is not mentioned in the available models. In addition, most of the previous Monte Carlo formulations require the full percolation of electronic and ionic phases across the electrode, even though this is not a necessary condition for the formation of active TPB sites. There is also a distinct lack of information regarding the effect of current collector geometry on cell performance and TPB formation in the available modeling work. Finally, many previous models ignore the gas phase or do not address the specific requirement of percolating the gas, ionic and electronic phases from each TPB to the bulk gas, ion, and electron phases in order to make it functional. The current model and investigation address these limitations and provide insight into these aspects of the problem for the specific case of an SOFC composite electrode microstructure. The treatment of the per-colation of the gas phase and the activity of this phase (through processes such as sorbate transport) are assumed at the model's inception to be of critical importance to the formation of TPBs; their effects are then verified through various means of data analysis. This particular feature presents a novel understanding of the microstructure requirements for TPB formation in composite electrodes, and provides an improvement to the knowledge previously developed. Although the effect of this gas phase is most directly linked to the formation of TPBs, as mentioned above, the conductivity of any given electrode is fundamentally dependent on the TPB formation characteristics. The significance of percolation physics related to the gas phase as reported in this work can be subsequently used to determine effects on overall electrode performance.

Model description
The current model uses a Monte Carlo simulation method for investigation of the fuel cell interface, developed in MATLAB. The Monte Carlo method depends upon analysis of a statistically significant number of randomly generated interfaces to estimate the performance as a function of the interfacial design parameters. The two major portions of the model are comprised of (1) the formation of a randomized composite electrode structure and (2) the analysis of the randomized structure, which itself consists of two types of analysis: (1) morphological analysis and (2) electrochemical performance analysis, though this analysis is not presented in the current investigation. The software was developed to create and analyze two-dimensional representations of the fuel cell microstructure with certain features defined by user input. These user-defined values included the relative volume fractions of the included phases of materials including gas phase and all conducting phases, the current collector and gas delivery boundary structure, the particle size and shape distribution for each of the materials included, and the conductivity of any conducting phases. For most cases, the particles were assumed to be monodispersed in size and shape and user inputs could therefore always be met exactly by randomized structures. However, for the cases of non-uniform particle size and shape, every combination of particle volume fractions and size and shape distributions was not guaranteed to enable a multitude of random realizations in the interface structure. Thus, the model allows for discrepancies between the user inputs and the realizations, accounting for and reporting errors between inputs and simulated parameters. In empirical investigations and commercial production of fuel cells, control of the relative amounts of each material is more easily achieved than control of particle sizes, which is often realized as a distribution of sizes about the target size. Thus, the model was coded to exactly preserve the user's input volume fractions, but not the size and shape distributions. A sample of a randomized monodispersed interface structure that is characteristic of those utilized in the majority of this study is shown in Fig. 1, with boundary conditions attached. It should be noted from Fig. 1 that the inclusion of gas phase elements in the electrode-electrolyte interface microstructure allows for the existence of large-sized pores, which can significantly alter the morphology and performance of the interface.

Morphological microstructure analysis
Although a TPB is defined locally as the point where all three phases necessary to support electrochemical reaction are present, on the larger scale of the entire electrode structure, there are additional requirements for the formation of an active TPB. In order to guarantee the continuous availability of electrons and ions at the TPB, there must be a path present from their source (or sink) to the TPB. Thus, electronically and ionically conducting phases must percolate through the randomized electrode structure to the point of the TPB. This percolating path must also have an origin with at least one particle of the appropriate type on the boundaries of the electrode, which becomes critical when considering varying current collector geometries. In addition, it is not assumed in this model that the gas phase is readily available throughout the entirety of the electrode. Instead, it too must percolate between the current collector boundary and the TPB in order to ensure the presence of the gaseous species, account for macropores, and avoid concentration polarizations. Thus, the formation of an active TPB in the current model is governed by the successful percolation of these three phases of materials to a common point in the randomized electrode structure. In this investigation, percolation has been assumed to occur not only through particle edges as in the classical case but also through particle corners. This functionality has been included as a means to approximate the necking and interconnectivity between like particles that occurs during high temperature manufacturing steps of the cell. Thus, the model was developed to count two classes of TPB. The first is termed a "potential" TPB, which simply satisfies the criterion of adjacency of necessary phases. Given the relatively simple nature of this calculation, it was possible to identify these TPBs immediately after the structure of the electrode was randomly generated by the model. The software was coded to recognize TPBs within the bulk of the electrode matrix as well as on the edges with the boundary conditions and the lateral edges in specialized cases. In this investigation, TPBs were only counted at corners of particles, which is the two-dimensional analog of the three-dimensional TPB length. For the cases when TPBs could be formed by only two particles, e.g., when an MEIC particle contacts a gas particle, TPBs were still counted only at the corners.
After identification of the potential TPBs, the model then searched the electrode to identify percolating paths for each of the three phases. Whenever MEIC particles were included in the structure of the interface, they were treated as having the possibility to be a part of both the electronic and ionic conducting paths. In order to facilitate coding of the software, paths were denoted in matrices with a similar structure to the original randomized electrode matrix, thereby allowing the same algorithm to be executed for locating potential and active TPBs, with the input matrices determining the type of analysis. Therefore, active TPBs were located within the electrode structure by comparing adjacent cells in the path matrices, as if the path matrices were overlaid to reveal only the percolating portion of the structure. As with the potential TPBs, active TPBs were located within the bulk and along the boundaries of the electrode-electrolyte interface.
Overall, two versions of the morphological analysis were developed. The model as described above represented the first, simple version of the code. An additional version was created to account for the effect of sorbate transport. Sorbate transport is a well-documented phenomenon [21] whereby various phases of the gases involved in the electrochemical reaction can adsorb or absorb onto or into the solid phases and move by various processes to participate in electrochemical reactions at a nearby TPB. Examples of the sorbate transport processes are shown in Fig. 2. The overall effect of these mechanisms is to extend the TPB in the three-dimensional case from an active TPB line length to an effective TPB area. In order to provide a simple means of accounting for this phenomenon, the modified version of the software allowed a TPB to be formed (in both the potential and the real case) when a gas particle was one cell removed from the interface between the electronic and ionic phases. Thus, this second version of the model checked every particle corner loca-tion in the electrode structure twice in order to determine the total number of TPBs.
The final morphological feature analyzed was the depth of the TPB within the composite electrode. For this investigation, the matrices that denoted the locations of TPBs were analyzed so that the distance between the current collector (or bulk ionic phase) and the TPB could be tracked. The analysis stored the layer number of each TPB, with the convention that layer 1 corresponds to the interface between the current collector and the electrode, with layer number increasing with every horizontal grid line in the direction of the electrolyte. Thus, for a 10 × 10 particle electrode, the layer numbers ranged from 1 to 11.

Model execution
Given the randomized nature of the structure that the model simulates, a Monte Carlo solution method was chosen for this investigation. In order to obtain statistical significance, 10,000 random interface instances were analyzed for every set of userdefined inputs. Thus, for each of these 10,000 iterations, a randomized electrode structure was formed, TPBs and paths were located, conductivity was calculated, and TPB depth was tracked. The selection of 10,000 iterations is discussed below. For the investigations shown here, 10 × 10 matrices were utilized for every realization of the matrix, the selection of which is also discussed below. Unless otherwise noted, the model was run assuming a 30% porosity in the electrode, with electronic and ionic conductor volume fractions composing the remaining 70%, either individually or summed together. Finally, various current collector boundary conditions were investigated as indicated in Fig. 3. The convention for the naming of the boundary structures indicates the relative number of gas and electronic particles lying over the randomized interface (the central 10 particles) followed by descriptive symbols of the order of the particle types, and in some cases, the location (middle or side) of key particles. Note that although the model was developed with  the capability to analyze particle size and shape, these aspects are not included in this paper for brevity. A discussion of the impacts of particle size will be presented in a future investigation.
Additional analysis of averages and standard deviations were performed in Excel, while statistical analysis involving ANOVA and significance tests was performed in Stat-Ease's Design Expert Software [24], using general factorial analysis types. In addition, for all statistical tests, any transformation employed was chosen to best satisfy the Box Cox test [24] for all factors simultaneously, utilizing the software's suggestions as guidelines.

Domain size
As mentioned previously, the simulations carried out for this investigation were based on 100-cell square matrices of particles. The determination of this domain size was not arbitrary; rather, it was a result of a preliminary investigation into the percolation aspects of the problem. Prior to acquiring the majority of the data for either TPB counts or conductivity of the composite electrodes, a sample set of TPB counts was collected for both a 100-cell square matrix and a 400-cell square matrix. Note that for proper execution, the algorithm required a number of particles which is a square multiple of 10; thus, 400 cells was the smallest matrix larger than 100 cells which was appropriate for the investigation. Averages of the real TPB count for a 30% porous electrode, with the idealized 1:1 boundary condition, shown in Fig. 4, were compared for the two domain sizes over a range of electronic conductor volume fractions from 0 to 70%. The 1:1 condition was chosen as it is the least limiting of those investigated, and it was desired that the results for this portion of the investigation be dependent only upon the volume fraction and domain size parametric variations. By ANOVA, it was found that the electronic volume fraction for the 100-cell domain was a significant contributor to the active TPB count, whereas domain size was not, with a p-value of 0.9201. This is shown in Table 1, where it can be seen that the statistics for the 400-cell domain do not differ substantially from the 100-cell domain, as this was the only factor included in the model for this analysis. Note that, per the standard ANOVA analysis, the sum of squares, degrees of freedom (d.f.), and mean square, were utilized to calculate an F-value for the parameters included in the model. This Fvalue was then used to find the corresponding p-value, which indicated the probability that the variance was significantly due to the factor itself rather than random error. A p-value of less than 0.1 was taken as significant for this study.
When considered together, as in Table 2, the two factors contributed to a model with significant dependence on the include variables, as is shown for the very low p-value of the model where factors A and B (Models A and B) are included. However, domain size remained insignificant as an individual factor, indicated by the high p-value (0.7488) of this factor (factor A). By contrast, the composition of the electrode (factor B) was shown to be a highly significant factor, given its extremely low p-value. This demonstrated that the electronic volume fraction was a major contributor in TPB formation and was dominating over the domain size as a pertinent factor. It was therefore concluded that the 100-cell domain size was sufficient for the current investigation and did not artificially limit the effectiveness of percolation in the composite electrode of interest.
In order to illustrate this conclusion, the active TPB counts for the 100-, 400-, and select 1600-cell matrices were normalized to their individual peak values and compared with respect to electrode composition. As can be seen in Fig. 5, there was a consistent dependence of the TPB counts on electrode composition, but the normalized values for each domain size were extremely similar, indicating trends that were not significantly dependent upon the domain size itself. It should be noted that the potential TPBs could also have been investigated, but these rely less heavily on percolation features, and it is therefore implicit that if active TPB counts are independent of domain size, then so should potential TPB counts. In addition, the domain size found to be sufficient in this investigation is in agreement with previous findings [20], although there is currently some contention over the validity of this finding noted in previous studies [15].

Number of randomized iterations
Many previous investigations have studied similar features of the randomized fuel cell electrode-electrolyte interface microstructure utilizing a small number of randomized seed variables on the order of 10 or less [16,19,20,22]. Given the size    of the simulation domains and the number of variables in the previous studies and the current study, this number of random seeds is quite low. In order to establish the number of iterations required to assess the impacts of material volume fractions on TPB counts, the average and standard deviation of TPB counts over a full range of volume fractions were investigated for varying numbers of simulations. A similar ANOVA analysis on the results described above was conducted to determine sensitivity to the number of iterations. It was found through this analysis that the number of iterations, though not a highly significant fac-tor, was only marginally so, with a p-value just over 0.1. Thus, the results were compared directly. As shown in Table 3, the percentage error in the average TPB counts was still significant for certain random mixtures utilizing 1000-point data compared to the 10,000-point data. For example, at 60% electronic volume fraction, there was still ∼22.5% error between the 1000 and 10,000 iteration averages, which would conventionally be regarded as significantly different. However, there was an overall trend for the error to decrease between increasing orders of magnitude of the number of iterations, indicating improvement with increasing numbers of iterations. In addition, the same trend is present in the standard deviation results, which accounts more for the iteration-to-iteration variability. Note that the 1000-point deviation data nearly approaches that of the 10,000-point data, with significant error between the 1000 and 10,000 iteration cases only present at 60 and 70% electronic volume fraction. However, it is clear from the results of Table 3 that neither 10 nor 100 random iterations adequately capture the same degree of variation as the 10,000 random iterations. Therefore, it was decided that although 1000 executions of code could possibly reduce execution time and closely match the results that would be obtained with 10,000 executions, it would still lead to slightly different average TPB counts when compared to results for 10,000 random iteration cases. Thus, all data points presented in the current paper represent the average over 10,000 randomized electrode structures, with error bars representing one standard deviation.
It should be noted that this verification process for the physical domain size and the number of randomized seeds provides confidence in the statistical significance of the results presented in this work. As mentioned, many previous works have used much smaller numbers of random iterations, often without discussion as to the selection process for this value. By verifying the model execution parameters as they have been above, the execution and quality of the data of the current study have been confirmed. It should be noted that other studies have been carried out on models with higher dimensionality (3D), but the choice of a 2D model with increased numbers of iterations and an appropriately sized simulation domain have imparted the ability to gather highly accurate statistics, with a reasonable execution time for the simulation software.

Relationship between real and potential TPBs
The same preliminary investigation that revealed the independence of active TPB formation on domain size at the 100-cell level and the need for a large number of randomized seeds provided useful insights into the formation of active and potential TPBs. Fig. 6 shows the trends in active and potential TPB counts with respect to electronic volume fraction (note that ionic volume fraction is equal to 70% minus the electronic volume fraction at any given point). As mentioned above, these results are produced by 100-cell, 10,000 random iteration simulations. It is immediately apparent that potential TPB counts peak at a 50/50 mixture of ionic and electronic conductor as expected, not requiring any special consideration of percolation or path formation. Thus, the most interfaces between these particles, at a constant gas volume fraction of 30%, are formed when the electrode is composed of just as many electronic conductor particles as ionic conductor particles. However, this is clearly not the case for active TPB counts. For the particular, uniform 1:1 boundary condition, the optimal number of active TPBs is actually formed when there are no electronically conducting particles in the electrode, thereby requiring that all TPBs are formed directly on the current collector. It should be noted that even at their peaks, the number of active and potential TPB counts are quite low. With the 10 × 10 cell geometry used in this study, the total number of corners where a TPB could form was 121. Only slightly more than onethird of these locations become potential TPBs at the optimal point, and the highest count of active TPBs is less than one-tenth of the number of available locations. These results are at odds with a recent study by Schneider et al. [19], which found that, in three dimensions, the number of TPBs in a randomized electrode was a constant around 1.3 times the total number of solid particles in the simulation domain. Schneider et al. found that this behavior is independent of the volume fraction of ionic or electronic conductor, which seems to be counterintuitive to the idea of percolating phases to form a TPB. The Schneider study also showed that the ratio of active to potential TPBs was a constant around 80%, which is a level that is only in the current study results for the case when the electrode is entirely constructed of an ionic conductor (see Fig. 6). The fact that results match when each case considers only two percolating phases suggests that the differences are not due to model dimensionality (2D versus 3D).
The differences between the current study and previous literature results can be explained by a set of considerations concerning the formulation of the model. Possibly of primary importance, the Schneider study, as well as many others, was carried out in three dimensions, whereas the current study was only carried out in two. This difference can significantly alter the effectiveness of percolation. In single-phase percolation, the percolation limit (concentration at which a randomly percolating network is guaranteed to span a geometry) for a two-dimensional square lattice with site percolation (characterized by placement of connected particles or sites, rather than bonds, as in this and other studies) is ∼59%. For three-dimensional studies on simple cubic geometries, it is ∼31%. Thus, in three dimensions, percolation of a single phase is more readily achievable than in two dimensions. Similar results are observed in multi-phase percolation [23]. Note that percolation limits are defined for an infinite domain size, so their realization within finite domains is constrained by the finite nature of the domain. In addition, multi-phase percolation is not well-understood, and there are exceedingly few cases in which percolation limits have been found, and the study of three phases percolating at once is particularly limited. Although the inability to find these percolation limits does not imply that they do not exist, it can be surmised that, given the limited effect of single-phase percolation limits on finite domains, the effect of multi-phase percolation limits will be even smaller in these finite domains. Thus, although there will be a difference in the percolation characteristics between two-and three-dimensional domains, it is expected to be small.
However, even more fundamental to the percolation physics is the consideration of the gas phase percolation, which is only included in the current work. Increasing the number of percolating phases necessarily causes percolation to be more difficult, thereby shifting the percolation limit for each phase. Thus, the percolation limits for the two models are expected to differ primarily because of the increased number of phases that are required to percolate. The effect of the differing percolation requirements can be expected to have at least as large of an effect on the percolation limits as the dimensionality. The disagreement with the results of Schneieder et al. is an indication of the importance of percolating the gas phase, as considered in this study. The effect of gas phase percolation on TPB formation is more completely characterized and verified below.
In addition, logical consideration of the problem when including the gas phase percolation can reveal the source of the disagreement in the models and hints at interactions between phases. The ratio of real to potential TPBs can be considered to be an index of the effectiveness of overall percolation of the phases (as opposed to individual percolation characteristics). In two-phase systems, as the effectiveness of percolation of one phase decreases (due to a lower volume fraction of this phase), an increase in the percolation effectiveness of the other phase can be expected, with an equal magnitude. Thus, constant overall percolation effectiveness can be expected. However, in the three-phase percolation considered in the current work, the third phase, which is held at a constant volume fraction, can alter the effectiveness of percolation of each of the other phases. If the third phase is positioned within the structure such that it equally effects the two phases, then the same result is expected for the two phases; a constant percolation effectiveness should result, but the effectiveness would be lower than in the two-phase case. However, if the third phase affects one of the two remaining percolating phases more than the other, but still negatively affects both phases (it is impossible for it to have a positive effect on percolation of the other phases), then the overall percolation effectiveness would not be constant, because increased inclusion of the less-hindered phase will continually increase the overall percolation effectiveness to a larger degree than increased inclusion of the more-hindered phase, which must counteract the effect of the third phase.
It can also be seen that, for this boundary condition, there is no local maximum in active TPB counts aside from the point with no electronic conductor, but there is a local minimum. The minimum is slight, but is most evident in the active/potential TPB curve, which has a positive inflection around 60% electronic volume fraction, where active TPBs begin to increase while potential TPBs fall. This result is not counterintuitive; rather, it simply indicates that the percolation of phases at high electronic phase particle concentrations is slightly more effective than at middle-range volume fractions. Interestingly, the potential TPB curve is not exactly symmetrical, as there are slightly fewer TPBs formed with an all-electronic electrode than an allionic structure. The same effect is visible in the active TPBs, though to a much greater extent. As seen in the active/potential curve, this effect immediately indicates a difference in percolation effectiveness between these two extremes in electrode composition, an important result investigated below.

Role of sorbate transport
Because of the small numbers of TPBs, especially active TPBs, that were predicted by the model as shown in Fig. 6, the model refinement to account for sorbate transport was developed. With sorbate transport, the effective probable area for TPB formation should be extended for any given random distribution of particles within the electrode. It can be seen from comparison of Fig. 6 with Fig. 7 that this obviously occurred for the potential TPB count, where the peak TPB count with sorbate transport is almost 1.5 times that for the case without sorbate transport. However, the active TPB count is barely affected by the inclusion of sorbate transport physics. Interestingly, though, the active TPB counts for the sorbate transport case (Fig. 7) now show a local maximum between 15 and 20% electronic volume fraction. Although this is still lower than what is intuitively expected from a perspective of cell performance, its source relates well to percolation considerations, discussed below.
In order to determine the range of electronic volume fractions over which the sorbate transport feature caused a significant change, ANOVA analysis was carried out on the data sets with and without sorbate transport for different domains of the electronic phase volume fractions. It was found that the resultant p-value indicated sorbate transport was not significant only for simulations that contained less than 35% electronic phase volume fraction. Thus, for all volume fractions at or above 35% electronic conductor phase, the inclusion of sorbate transport physics significantly increases the average number of active TPBs. The ANOVA results at the 35% volume fraction cutoff are presented in Table 4, with the p-value of 0.0251 indicating the statistical significance of the result.

Dependence of TPB formation on depth in electrode
The dependence of active/potential TPB counts on electronic volume fraction, together with the composition of the electrode at peak active TPB counts, suggests a correlation between location in the electrode and likelihood of forming an active TPB. The composition of the electrode, and consideration of how well phases might percolate at this composition (as indicated by the active/potential curve), indicates that percolation was of particular importance and may not have been equivalent between phases. This would then lead to a location dependence of active TPB formation as less effective percolation of a phase would require TPBs to be formed closer to that phase's bulk. To determine whether or not there is a correlation between the location of a particular particle within the electrode structure and the likelihood of forming a TPB at that particle depth, the location in the electrode structure of all active TPBs from an equivalent set of simulations was tracked. The location was stored simply as the depth at which the TPB existed within the electrode, defined on a cell edge basis as mentioned previously. It was found that there was a correlation between TPB count and position in the electrode with higher TPB counts occurring near the air/current collector boundary. This correlation is expected due to the presence of two phases at this boundary, compared to the single ionic phase at the opposite boundary. This correlation is also dependent upon the overall electrode composition. In addition, because the inclusion of sorbate transport physics affects TPB formation, some correlation between this feature and the most likely position of TPBs in the electrode was desirable. Since this effect is intimately connected to the availability of gas phase, it is beneficial to first study the effect of sorbate transport for a case that is independent of variations in solid phase mixture fractions. To do this, the likely location of TPBs in the electrode was first investigated with an MEIC as the only solid phase in the electrode for both cases with and without sorbate transport. Afterwards, the composite electrode was studied only with sorbate transport included. The TPB count curves for the MEIC cases are presented in Fig. 8.
It is immediately apparent that the sorbate transport feature has a significant effect, as seen by comparing the potential TPB curves of Fig. 8a and b. In addition, note the very significant differences between the TPB counts associated with an MEIC electrode (Fig. 8) and those associated with the composite electrode (Figs. 6 and 7). Although the peak active TPB count is much higher in the MEIC case, note that it is still only about half of the total available locations in the entire electrode; thus, even a single-phase, porous electrode is not completely active throughout its entire structure, because of the need to percolate the gas phase as well. In the MEIC case, not only is the peak at a much higher TPB count with sorbate transport, but the composition of the peak structure is also significantly higher. It can also be seen that the asymmetry in the potential TPB curves is much more pronounced when sorbate transport is included. Interestingly, for the MEIC electrode, the asymmetry is the reverse of the composite electrode case, indicating a slight preference for TPB formation with a lower porosity. This is due to the fact that without any solid phase in the electrode, no TPBs can be formed, but with an all-solid electrode, some TPBs can still be formed at a gas boundary through the guaranteed adjacency of the solid phase to the gas-containing boundary. It can also be seen that, as expected, when an MEIC is considered, percolation is more effective, but has a very different character from the composite electrode case, with a much more pronounced real TPB peak at an almost equal mix between gas and solid phase, and no TPBs forming at all until around 25-30% solid phase volume fraction, which is significantly lower than the percolation limit of ∼59.3% in site percolation on a square lattice [23].
Many of these trends are explained by the TPB position results, shown for the MEIC electrode case with and without sorbate transport in Fig. 9 and for the composite electrode case with sorbate transport in Fig. 10. As can be seen from Fig. 9, the inclusion of sorbate transport tends to increase the TPB count on all layers for a given solid phase volume fraction. However, it does not do so uniformly. It should first be noticed that at the peak TPB condition, 55% MEIC, the TPBs are roughly uniformly distributed throughout the electrode, but as the electrode is filled with increasing numbers of solid particles, TPBs tend to aggregate towards the gas/current collector boundary. This results from the fact that the limited amount of gas phase in the electrode structure causes gas phase percolation to be less effective, limiting the extent of gas phase paths through the electrode. Moreover, the gas phase particles on this boundary, which could only affect layer 1 without sorbate transport, cause the TPBs in layer 2 to be almost identical when there is a very high solid phase volume fraction. This effect, however, did not extend beyond the second layer, as expected given the means of modeling sorbate transport physics in this study. However, it should be noticed that sorbate transport, even in the MEIC case, only significantly affected TPB count for the more extreme MEIC volume fractions, and more for the layers closest to the gas/current collector boundary.
The relationship is slightly more complex when considering the composite electrode, where the electronic phase has to percolate from the same side as the gas phase. As seen in Fig. 10, the trend for the position of TPBs changes from all TPBs being present on the gas/current collector boundary when there is no electronic conductor in the electrode to all TPBs on the electrolyte boundary when there is no ionic conductor. This is as expected since these are the only locations where all three phases touch in the respective cases, since one phase is missing in the bulk. It is critical to note, however, that there are between one and two orders of magnitude difference in the TPB counts for these two cases. This clearly indicates a difference in the effectiveness of percolation. Moreover, since one phase is absolutely non-percolating in each case, it can be seen that the percolation of the electronic and gas phases down to the electrolyte is less effective than the percolation of ionic phase up to the current collector. Since the boundary condition at the current collector is uniform and provides the same opportunity for path initiation to the gas and electronic phases, and the ability to percolate for any individual phase is equivalent, single-phase percolation considerations cannot explain this phenomenon. The only effect that can explain this result is competition between percolation of the gas and electronic phases to a potential TPB to create an active TPB. Note that, given the random placement of particles, potential TPBs were expected to show no dependence on electrode depth, and a sample of results revealed this to be true of the data obtained.
This trend of asymmetry between electronic and ionic phases is mirrored in the complementary electronic volume fractions of 15 and 55%. The 15% electronic volume fraction case contains a higher concentration of TPBs near the gas/current collector boundary than the corresponding 55% electronic volume fraction case contains near the electrolyte. In addition, the TPBs in the 15% case are not present in as many layers as they are in the 55% case, indicating a difficulty in simultaneously percolating the gas and electronic phases. Moreover, approximately uniform TPB formation occurs at 40% electronic conductor volume fraction, slightly more than an even mixture between ionic and electronic phases, indicating that proper percolation requires slightly more electronic conductor. However, it should be noted that proper percolation does not necessarily guarantee optimal TPB formation, as it is actually the 15% electronic conductor case that has the highest number of TPBs formed. Thus, it is clear that TPB formation is not simply a matter of classical percolation of phases. This can also be seen in the MEIC results; even though a uniform TPB count (indicating successful percolation of all phases through the electrode) is exhibited as low as 25% MEIC concentration, the peak TPB count does not occur until 55% MEIC, which also presents a uniform TPB distribution. Thus, although full percolation is necessary for active TPB formation, it is not sufficient to provide optimal active TPB results, and is particularly sub-optimal for the composite electrode case.

The effect of boundary conditions on TPB formation
While both boundary conditions affect the computations, the boundary condition characteristics at the gas/current collector can particularly limit active TPB formation since it is the origin of any percolating set of both electronic and gas particles. In the investigations above, this limiting behavior was eliminated by using the idealized boundary condition that has gas and electronic conductor present throughout. This structure is not physically possible, though, and analysis of more realistic structures is required. In order to achieve this, the model was run with the gas/current collector boundary conditions shown in Fig. 3. Fig. 11 presents a comparison of the peak active TPB counts for each boundary condition case. Note that peak active TPB count was not always achieved at the same electronic conductor volume fraction, as shown in Fig. 12. Fig. 11 shows, as expected, that the active TPB counts significantly depend upon the gas/current collector boundary conditions. Note that for the 2:2 GEE case, data was collected for 20, 30, and 40% porous structures, so the high TPB count of the 2:2 GEE case with 20% porosity is not indicative of better performance than the ideal boundary condition, but rather that the decreased porosity allowed for higher active TPB formation, which is itself an interesting result.
In order to determine how the boundary condition affects TPB formation, a number of results and metrics were investigated. Fig. 13 compares the results in active TPB counts between a highly unbalanced boundary condition and one where a single particle has been changed towards a more balanced boundary. It can be seen that when both gas and electronic particles are in deficit at the current collector (1:9 GE(S) and 1:9 EG(S), respectively), inclusion of one more particle of that type (1:4 GE and 1:4 EG, respectively) increases TPB counts. Interestingly, when sorbate transport is included, the case with an electronic particle deficit gains more by inclusion of the deficit particle than when gas is in deficit; by contrast, without sorbate transport, the Fig. 12. Electronic volume fraction corresponding to the peak active TPB case for each of the gas/current collector boundary conditions cases. Fig. 13. Sensitivity of active TPB counts to an incremental increase in balancing the gas/current collector boundary condition phase composition. gas particle deficit shows a larger increase with inclusion of the deficit particle. However, this may be affected by the difference in peak TPB counts between cases with and without sorbate transport in general. In addition, the gas phase particles were held to a constant 30% of the overall volume fraction, which was not an optimized value, as it was for the electronic phase. Thus, inclusion of electronic phase in the boundary condition could be expected to not have quite as large of an effect as inclusion of the gas phase in the boundary. These results taken together indicate that a balance needs to be designed into the current collector gas delivery boundary. However, the gain that can be obtained by achieving a more balanced structure can be dependent upon the physics of the TPB reactions.
The composition of the gas/current collector boundary (i.e., relative amounts of electronic and gas phases in the boundary condition) is not the only factor of concern. As noted previously, active TPBs require percolation of both gas and electronic phases and can tend to aggregate at the boundaries. At the gas/current collector boundary, this requires that a gas particle be present adjacent to an electronically conducting particle. Thus, the structure of the boundary condition must also have an effect. In order to qualitatively visualize this, the results of Fig. 11 are rearranged in Fig. 14 in order of increased intermixing of the gas and electronic conductor phases. There is clearly a monotonic relationship between peak active TPB counts and the intermixing of the phases in the gas/current collector boundary condition. The degree of intermixing is determined by assessing the number of Fig. 14. Dependence of active TPB formation on particle intermixing in gas/current collector boundary condition. gas-electronic interfaces in the current collector and the spacing between these interfaces. High intermixing is characterized by high numbers of interfaces and smaller average spacing between interfaces. ANOVA analysis performed on these factors reveals that both of them are significant, as shown in Table 5. Furthermore, the results of Figs. 15 and 16 also show that both factors cause a nearly monotonic increase in the peak number of active TPBs with increased intermixing, both with and without sorbate transport. Thus, from the perspective of increasing active TPB count, the optimal structure of the gas/current collector design would introduce high intermixing of phases, which can be characterized and controlled by the independent factors considered in Table 5 and Figs. 15 and 16. However, this must be balanced with cost of manufacturing and gas flow restrictions, which both may be reduced with smaller numbers of gas channels (i.e., low intermixing in the gas/current collector boundary condition).

Summary and conclusions
A Monte Carlo model has been developed and utilized to identify and characterize the factors that control TPB formation in an SOFC electrode. The nature of SOFC electrode design can be analyzed by classical percolation theory, but the process of creating potential and active triple phase boundaries is not fully characterized by percolation physics alone. The requirement for (1) multiple phases to percolate to an intermediate point rather than across the domain to both edges, along with (2) the competition between the gas and electronically conducting phase percolation from the same boundary, and (3) the need to produce continuous paths from TPBs to the bulk gas, electronic and ionic phases to make each active each contributed independently to the unique physics apparent in the SOFC electrode. These requirements for percolating all three phases through the electrode, and the development of competition between phases are aspects that have not previously been characterized in the literature, yet have been shown in this investigation to have a significant impact. This study presents statistically significant results that thoroughly investigate active TPB formation that form the basis for future investigation of other pertinent physics such as those that govern overall contributions to current from the TPBs that are identified by the current model. This work therefore presents an important first step towards the understanding of full electrode performance characterization as governed by all requirements for three phase percolation to electrode-electrolyte interfaces.
The model was also used to show that physical processes in the fuel cell, such as sorbate transport, significantly alter the effectiveness of active TPB formation in a composite electrode. Sorbate transport always improves the formation of TPBs, regardless of electrode composition or gas/current collector boundary conditions. However, the improvement is not uniform across all cases, depending largely upon the amount of gas phase already present at the gas/current collector boundary. It also does not necessarily improve the percolation characteristics overall because active TPB formation in a fuel cell composite electrode is dominated by the characteristics of continuous percolation of multiple (3) phases to each TPB.
In addition, the model shows that gas/current collector boundary conditions significantly affect percolation. Structures that reasonably approximate practical interconnect (bipolar plate) designs reveal a need for proper design of gas channels and electronically conducting solid material in order to achieve high active TPB counts. This result is most closely tied to the need to percolate two phases from this boundary, which compete with each other in order to form continuously percolating pathways. Boundary conditions with significant phase disparity are shown to reduce TPB counts by up to an order of magnitude.
TPB formation in an SOFC composite electrode is most effective when the electronic volume fraction is quite low, usually 15% or less, and when the gas/current collector boundary contains highly intermixed gas and electronic conductor phases (i.e., a preference for narrower channels and contacts). Gas phase percolation can significantly alter the TPB formation, a factor that has not been thoroughly investigated in previous studies. Statistical analyses have proven the significance of these factors in determining active TPB counts. This optimum composition, however, is based solely on the consideration of the likelihood of the TPB formation. The choice of composition based on full performance characteristics requires the investigation of all pertinent electrode physics, which will be the subject of future study utilizing the current model.
In summary, this study presents statistically significant results from a Monte Carlo model that has revealed a multitude of methods by which features of electrode construction affect the formation of active TPBs. These features include: (1) accounting for the possibility that gas phase is not uniformly present at all locations in the electrode and can form structures just as significant as the structures of the solid phases, (2) development of percolating paths for the gas phases as well as the solid phases, (3) the enhanced effectiveness of the gas phase with respect to the solid phases, attributed to the action of sorbate transport, and (4) the effect of the current collector/gas boundary in minimizing the competition between gas and electronic phase percolation through the electrode.
While this TPB formation behavior has been characterized herein, previous investigations [15,19] suggest that optimal fuel cell performance characteristics (such as electrode conductivity) may not occur for conditions of maximum TPB formation. To this end, the current model has been developed (as mentioned previously) with capabilities to calculate the conductance associated with each TPB that will be used in future analyses.