The central premise of porous electrodes is to make more surface area available for reactions. However, the convoluted pore network of such reactors exacerbates the transport of reacting species. Tortuosity is a measure of such transport distortion and is conventionally expressed in terms of porosity (the fraction of electrode volume occupied by liquid-filled pores). Such an approach is overly simplistic and falls short of accounting for spatial variabilities characteristic of electrode samples. These networks are defined by multiple features such as size distribution, connectivity, and pore morphology, none of which are explicitly considered in a porosity based interpretation, thus limiting predictability. We propose a recourse using a two-point correlation function that deconstructs the pore network into its essential attributes. Such a quantitative representation is mapped to the transport response of these networks. Given the explicit treatment of pore network geometry, this approach provides a consistent treatment of three-dimensionalities such as inhomogeneity and anisotropy. Three-dimensional (3D) tomograms of Li-ion battery electrodes are studied to characterize the efficacy of the proposed approach. The proposed approach is applicable to abstracting effective properties related to different transport modes in porous fluid networks.

## INTRODUCTION

Electrochemical devices facilitate the interconversion of chemical and electrical energy. Such a transformation is identified as electrochemical reactions,^{1,2} which by their very nature occur at the electrode-electrolyte interface, often referred to as the active surface. To ensure enhanced interfacial contact between the solid and the electrolyte phases, the electrodes are made porous. These porous electrodes behave as electrochemical reactors. However, the porous nature of the electrodes not only affects the reaction kinetics (by providing more area for reactions in the same volume) but also impedes the ionic transport to reaction sites (referred to as transport distortion here), in turn limiting its functionality. The conundrum of increasing interfacial contact without convoluting pore network is common to porous electrodes^{3,4} and spans across a range of electrochemical systems.^{5–14}

Such porous electrodes exhibit coupled spatial variations in geometrical arrangement and physicochemical evolution (physical fields, e.g., temperature; chemical fields, e.g., concentration), with associated temporal complexities. The field inhomogeneities attenuate at smaller lengths. At very small scales, however, geometrical distribution is disorganized. Adhering to these two distinct variabilities, an intermediate length scale, *ℓ*, can be identified which is small enough for physicochemical fields to be locally uniform but large enough for statistically sufficient geometrical details to be sampled locally in volume ∼*ℓ*^{3}. In other words, *l*_{pore} ≪ *ℓ* ≪ *l*_{electrode}. For lengths larger than *ℓ*, fields vary due to the finite-time nature of transport processes. Such a construct, more commonly known as porous media theory or effective medium approximation,^{15–21} allows one to the abstract geometrical arrangement in terms of a few porous attributes, e.g., porosity, tortuosity, permeability, etc. *ℓ*^{3} is referred to as a Representative Elementary Volume (REV).^{19,22–24} It is possible for the electrode geometry to have long wavelength (larger than *ℓ*) spatial variations and lead to gradients in effective properties. Such an electrode is referred to as a heterogeneous electrode. (Philosophically, the REV concept is equivalent to the continuum definition in classical thermodynamics^{25–27} and associated field theories.^{28,29} For example, density—a continuum property—is meaningful over a large enough volume to accumulate enough molecular matter but small enough so that it can be treated as a point property. It allowed spatial variations over larger lengths.)

Following this chain of thought, it can be shown that the ionic current pathways experience distortion due to the specific geometrical arrangement of the porous electrode. The increased transport resistance relates to porosity and tortuosity as relevant effective properties at the REV scale, *ℓ*. Since porosity is easier to estimate, the conventional practice is to identify analytical expressions for other properties in terms of porosity—volume fraction of electrolyte filled pores^{30–36} (and, in general, as a function of phase fractions^{9,37–39} when more than two phases are present in the electrode structure; such an interpretation is characteristic of transport interactions in porous media and is not just limited to porous electrodes, e.g., Refs. 40–42). A classic example is the Bruggeman relation^{43} that correlates tortuosity, *τ*, to porosity, *ε*, for a porous arrangement of monodispersed nonoverlapping spherical solids,

The analytical simplicity of such correlations has two shortcomings:

Porosity is an adirectional property; hence, the extension for anisotropic electrodes is ambiguous.

A large number of structures

^{9,37,38}—typically of the order of 100—need to be analyzed to develop such a correlation. Of different tortuosity estimation approaches, direct numerical simulation of the solute diffusion (refer to Appendix A) is sufficiently valid across different pore morphologies, which multiplicates the computational overhead to explore tortuosity trends for each microstructure class.

Fundamentally, a pore network is defined by pore size distribution (in general, a four-dimensional histogram in three spatial coordinates—*x*, *y*, *z*, and pore dimension—*l*_{pore}), connectivity, and pore morphology. Hence, the transport distortion relates to these three attributes of any electrode pore network. The conventional approach of porosity correlations reduces the size distribution to a single value, i.e., porosity (a zeroth order estimate per statistical terminology) and implicitly accounts for connectivity via studying multiple realizations of similar structures (during correlation development stage). Morphology variations lead to different correlations, for example, an electrode made up of spherical particles^{37} has a distinct pore network as compared to a fibrous structure^{10} and correspondingly disparate expressions.

We propose an alternative treatment for abstracting electrolyte transport distortion in porous electrodes which expresses the pore network using a two-point correlation function.^{45–49} The two-point correlation function quantifies both the size distribution and connectivity (deconstructing pore network into its essential constituents). Given the detailed information about the pore network (as compared to porosity), one should be able to construct an accurate relationship between the transport distortion (as expressed by tortuosity) and the quantitative pore network (as described by the two-point correlation function). Since the two-point correlation function incorporates the directional information, it facilitates an intuitive route to abstracting anisotropy. We exemplify the strength of such an approach by considering tomograms of Lithium-ion battery cathodes^{44} (an illustrative reconstruction is shown in Fig. 1). These electrodes are large enough to contain spatial variations of particulate matter (pore morphology remains invariant within the same electrode). We show that an appreciably smaller number of REVs need to be studied explicitly to learn the aforementioned relationship. The uniqueness of such an approach lies in the fact that the transport distortion is explicitly codified in terms of pore network, which eliminates the need for large data to ascertain realistic predictability.

## QUANTITATIVE DESCRIPTION OF PORE NETWORK

The three dimensional electrode structures are obtained either through imaging^{50–52} or stochastic generation.^{53–55} The essential geometrical information is contained as a phase (i.e., material) identifier on a voxelated grid (voxel ≡ smallest resolved volume). For example, the electrode microstructure shown in Fig. 1 contains two phases: active material and pores (from the open source tomograms of Ebner *et al.*,^{44} the binarized dataset corresponding to 90 wt. % active material and 600 bars calendaring pressure is used to discuss the mathematical preliminaries). The indicator function can be consistently translated to a probability map assigning the likelihood of finding a particular phase at a given voxel location. Such a probabilistic representation is referred to as a one-point correlation function, $\varphi x\u2192|h$, and defined as

Here, $x\u2192$ is the vector identifying spatial points (it has a discrete representation on the voxelated grid), and h identifies a particular material phase. If the structure has (*N*_{x}, *N*_{y}, *N*_{z}) voxel dimensions in respective coordinate directions, the indicator function and one-point correlation function are 3D datasets with (*N*_{x}, *N*_{y}, *N*_{z}) dimensions. Figure 2(a) shows an REV from the tomogram presented in Fig. 1 (REV identification has been discussed in the authors’ previous work, e.g., Ref. 37). The rendering in Fig. 2(a) colors phases per corresponding one-point correlation function. The surface enveloping phase h is the $\varphi x\u2192|h=0.5$ isosurface. Note that the one-point correlation function is a scalar function of a vector variable, $x\u2192$.

The one-point correlation function implicitly contains the connectivity mapping, but such information is not readily apparent. The higher order correlation functions explicitly quantify the connectivity of different types. The two-point correlation function, $\psi r\u2192|h,g$ (shorthand notation $\psi r\u2192h,g$), characterizes^{45–49} the likelihood of phases h and g to be separated by a distance *r* in the direction $r^$. For any vector, $r\u2192$, $r\u2261magnituder\u2192=r\u2192$ and $r^\u2261directionr\u2192=r\u2192/r\u2192$. The two-point correlation function builds upon the information contained in the one-point correlation functions as

Here, Ω refers to the 3D structure under consideration (i.e., an REV for the present discussion). The following identities can be shown (*ε*^{h} is volume fraction of phase h):

The two-point correlation function also is a 3D matrix of scalars with voxel dimensions (*N*_{x}, *N*_{y}, *N*_{z}). The $\psi r\u2192pore,pore$ mapping for the REV in Fig. 2(a) is shown in Fig. 2(b). The calculation for $\psi r\u2192h,g$ [Eq. (3)] essentially checks every two voxels separated by $r\u2192$ in volume Ω if they have phases h and g and subsequently computes the fraction of positive occurrences, i.e., a frequentist probability. For a volume, Ω made up of *N*_{x} × *N*_{y} × *N*_{z} voxels, there are *N*_{x} × *N*_{y} × *N*_{z} pairs of two points separated by a given vector $r\u2192$. Additionally, in such a volume, there are *N*_{x} × *N*_{y} × *N*_{z} possible discrete separation vectors. Thus, the determination of $\psi r\u2192h,g$ involves $ONxNyNz2$ operations and can become prohibitive for realistic voxel dimensions. Also, for every separation vector, a few voxel pairs encounter the domain boundary, i.e., one of the voxels falls outside the voxel space {1, …, *N*_{x}} ∪ {1, …, *N*_{y}} ∪ {1, …, *N*_{z}}. Such exceptions are resolved by assuming periodic boundary conditions.^{45–49} The mathematical structure of the two-point correlation function allows for a faster calculation via Fast Fourier Transform (FFT)^{48} ( Appendix B summarizes algorithmic details).

Figure 2(b) pictorially presents the spatial structure of the two-point correlation function, $\psi r\u2192pore,pore$. Since the interest is in dissecting the pore network, $\psi r\u2192pore,pore$ is discussed subsequently. Equivalently, $\psi r\u2192active,active$ describes the active particle connectivity. On the other hand, $\psi r\u2192active,pore=\psi r\u2192pore,active$ embodies the interfacial geometry of active material–pore contact. Each point in Fig. 2(b) corresponds to a specific separation vector, with $r\u2192=0^$ at the centroid of the volume. Three perpendicular coordinate planes passing through this point as well as two-point correlation function distribution on these planes are highlighted. As a general trend, the two-point correlation function values decrease as separation grows. For small separations, the distribution appears to be nearly isotropic, while undulations are present at larger separations. $\psi r\u2192pore,pore$ approaches porosity in the limit $r\u2192\u21920^$ [as expected from Eq. (5)]. Such behavior indicates that the pore voxels are clustered together to form pores (small $r\u2192$ data). Nearly isotropic $\psi r\u2192pore,pore$ suggests that the pore morphology is approximately spherical. Intermediate separations show directionality in accordance with pore–to–pore connections. At even larger separations, $\psi r\u2192pore,pore$ drops, suggesting a gradually inferior correlation between two pores separated by a large distance. The concept of the two-point correlation function, $\psi r\u2192h,g$, is similar to the radial distribution function, $\chi rh,g$, commonly used in statistical mechanics treatments.^{27,56–59} It can be shown that the two are related [Eq. (7)]. The two-point correlation function contains directional dependence of the separation vector, and if averaging is carried out over all the directions at a fixed separation, one recovers the radial distribution function. The separation vector has equivalent representations based on the choice of coordinate systems. Such equivalence between the Cartesian and spherical coordinate systems can be expressed as $r\u2192=x,y,z=rcos\u2009\phi \u2009sin\u2009\theta ,sin\u2009\phi \u2009sin\u2009\theta ,cos\u2009\theta \u2261r,\phi ,\theta $, where *φ* ∈ [0, 2*π*) is the azimuthal angle and *θ* ∈ [0, *π*] is the polar angle,

## MAPPING PORE NETWORK TO TRANSPORT DISTORTION

The two-point correlations have been traditionally used to characterize ordering in complex structures given its probabilistic interpretation.^{60} As discussed in Fig. 2(b), essential information is contained in smaller separation vectors. In turn, the correlation functions provide a mathematical structure for dimension reduction, leading to applications^{45–49,61,62} in structural comparison, databasing, low dimensional visualization, and reconstruction. Many properties of multiphase composites depend on the structural arrangement of their constituents, and in turn, the connection between such effective properties and correlation functions is intuitively clear (the two-point correlation function is essentially a structural description). However, the explicit relationships are not so apparent given the dimensional disparity between the two. For example, the elastic response of composite material is strongly dependent on the internal arrangement. The mechanical properties defining such a response are a handful of elastic constants, while the correlation function contains considerably greater degrees of freedom (each separation vector represents a different degree of freedom). Such a dilemma of dimensional matching obfuscates the application of the aforementioned intuition, resulting in relatively fewer studies along similar lines.^{49,63,64}

The discrete two-point correlation, $\psi r\u2192pore,pore$, contains values of a scalar $\psi r\u2192pore,pore$ at different spatial locations. Each of these *N*_{x} × *N*_{y} × *N*_{z} locations is unique. In turn, $\psi r\u2192pore,pore$ is a sequence of numbers, say, *a*_{0}*a*_{1}*a*_{2}⋯*a*_{(N−1)}, where *N* = *N*_{x}*N*_{y}*N*_{z},

The index *p* refers to a unique position vector. As discussed earlier, the information is densely packed at smaller separations. Moreover, many directions are strongly correlated. Both these facts suggest that a lower dimensional representation is possible for the two-point correlation function.

For the structures of the same type, i.e., exhibiting qualitatively similar pore networks, one expects common features. If *M* structures of the same type are studied, one can identify class attributes which represent the group behavior,

Notation $apu$ refers to the two-point correlation function value at *p*th separation vector from the *u*th structure. In turn, the structure-to-structure variability is retained in $\delta \psi pu=\psi ppore,poreu\u2212\psi \xafppore,pore$ vectors,

where

The structure specific features vector *δψ*|_{u} is better suited for the dimension reduction since the class attributes are first filtered out. These data can be summarized in a matrix form,

Compute the corresponding covariance matrix, $CN\xd7N$, treating each row of Δ_{M×N} as a set of observations,

where each entry has the form

Physically speaking, the covariance matrix, $C$, quantifies the structure–to–structure variability in two-point correlation function, $\psi ppore,pore$. The diagonal entry $Cpp$ specifies the variance in structure specific features for the *p*th direction. The off-diagonal entry $Cpq,p\u2260q$ examines if the directions *p* and *q* contain correlated information, with a larger value signifying stronger correlation. These $Cpq,p\u2260q$ entries can be positive or negative. In the subsequent step, an eigen decomposition^{65} is carried out to identify hybrid coordinates that are unrelated. For each of these hybrid directions, the information content scales with the corresponding eigenvalues. These directions are constructed via a linear combination of the original spatial directions, and the specific combination rule is encoded in respective eigenvectors. Let the eigenvalues and directions be, respectively, denoted as *λ*_{j} and *e*^{j}’s.

The pore network is formed in electrodes during the electrode preparation stage.^{66} In addition to being dependent on the morphology of the solid phase,^{31,38,39} the pore shapes also depend on the presence of additional phases and processing factors.^{67,68} Thus, structurally speaking, electrodes with the same constituents, compositions, and fabricated under similar conditions belong to the same class. Multiple REVs from the electrode shown in Fig. 1 are studied to obtain corresponding two-point correlations as well as exact tortuosity values ( Appendix A). Successively larger learning sets are prepared with 8, 16, and 24 structures based on such calculations. For every set, the two-point correlations are analyzed to dissect class attributes and structure specific features [Eqs. (9) and (10)]. The analysis of the structure specific features [Eq. (11)] reveals the eigendirections containing maximum information. The first two directions (for this structure class) contained more than 99% information and are the independent (hybrid) variables containing pore network information,

Such independent variables, *y*^{k}, are more commonly referred to as principal components. These are subsequently used to build a formal relationship of the form,

The specific analytical nature of the relationship [Eq. (15)] depends on the underlying trends. Figure 3(a) presents the relationship *τ* = f(*y*^{1}, *y*^{2}) as a contour surface, along with data points. Figure 3(b) compares the predictions against data when three different learning sets are used and reveals that the functional relationship improves only slightly going from 16 to 24 structures. Henceforth, 16 structures are used to construct the expression (15) for tortuosity. It should be noted that the procedure outlined in Eqs. (9)–(14) is mathematically similar to Principal Component Analysis (PCA) or Singular Value Decomposition (SVD) as referenced in the specific literature.

## PORE NETWORK BASED TORTUOSITY PREDICTIONS

The mathematical steps discussed in Eqs. (9)–(14) facilitates a consistent approach to group the essential pore network features along two independent directions. It does not discard useful information but reexpresses data in a condensed manner. It is essentially a dimension reduction exercise. The electrode set used to build the relationship (15) contains $O100$ REVs, while the study of 16 (an order of magnitude less) provides a reasonable convergence. Figure 4(a) plots the (exact) tortuosity distribution for this electrode as a function of porosity. The anisotropy in the pore network transport is clearly visible. Moreover, for the same porosity, there are multiple tortuosity values, further debating the universality of porosity-based correlations. Figure 4(b) compares the tortuosity predictions from pore network analysis against exact values obtained through Direct Numerical Simulations (DNSs). The predictive accuracy, *ϵ*, is defined as

Table I reports an inaccuracy of about 5% (per prediction) for *x*–directional tortuosity prediction. Thus, the proposed approach sufficiently captures the pore network differences. Visually such an accuracy is reflected in scatter-points being tightly grouped around the 45° line [Fig. 4(b)]. Any point on this line represents *τ*_{learning} = *τ*_{exact}. Figure 4(b) conveys that the general nature of tortuosity distribution is fairly captured. An outlier point on such a plot would exhibit a negative correlation, i.e., small *τ*_{learning} for high *τ*_{exact}, or vice versa. The procedure for predicting tortuosity based on the formal relationship (15) is as follows:

- Step 1.
Compute one-point correlation function, $\varphi x\u2192|pore$, as per Eq. (2)

- Step 2.
Estimate two-point correlation function, $\psi r\u2192pore,pore$ ( Appendix B)

- Step 3.
Obtain structure specific feature vector, $\delta \psi =\psi r\u2192pore,pore\u2212\psi \xafr\u2192pore,pore$ [Eq. (10)]

- Step 4.
Compute principal components,

*y*^{k}=*δψ*·*e*^{k}[Eq. (14)] - Step 5.
Use analytical relation

*τ*=*f*(*y*^{k}) to estimate tortuosity [Eq. (15)]

Direction . | Accuracy, ϵ
. |
---|---|

X | 4.43 × 10^{−2} |

y | 5.67 × 10^{−2} |

z | 4.95 × 10^{−2} |

Direction . | Accuracy, ϵ
. |
---|---|

X | 4.43 × 10^{−2} |

y | 5.67 × 10^{−2} |

z | 4.95 × 10^{−2} |

Similarly, pore network → tortuosity mappings are constructed for tortuosity in other directions and subsequent predictions reported in Figs. 4(c) and 4(d) and Table I.

The discrete two-point correlation, $\psi ppore,pore$ [Eq. (8)], is an ordered sequence of *a*_{p} values, where sequence index *p* maps to Cartesian indices (*i*, *j*, *k*). Note that *p* ∈ [0, (*N* − 1)], *i* ∈ [0, (*N*_{x} − 1)], *j* ∈ [0, (*N*_{y} − 1)], *k* ∈ [0, (*N*_{z} − 1)] and *N* = *N*_{x}*N*_{y}*N*_{z}. The two-point correlation functions for transport in *y* and *z* directions are constructed by altering this mapping *p* ← (*i*, *j*, *k*) to *p* ← (*j*, *k*, *i*) and *p* ← (*k*, *i*, *j*), respectively (without having to reevaluate for every coordinate direction). Figures 4(c) and 4(d) plot the error in individual tortuosities (the average error is absolute 5% as documented in Table I). Individual errors are bound within ±10%, with most of them clustered much closer to the zero line, resulting in a smaller averaged inaccuracy. It is interesting to observe that the anisotropy is captured without losing the accuracy of predictions.

As alluded earlier, the pore morphology could vary with processing changes during electrode fabrication. Porous electrodes are often compressed to increase their energy density.^{69,70} The active material phase is typically brittle and prone to fracture. However, before the stresses reach high enough values to initiate particle-scale fracture, the particles likely rearrange. Such a particle relocation should not change the nature of the pore network and only alter the size distribution and connectivity. In other words, the pore network → tortuosity mapping [e.g., Eq. (15)] remains accurate as far as the qualitative nature of the pore network is invariant. Hence, the loss of accuracy of such a mapping is a direct evidence of nontrivial changes such as particle fracturing. Figure 5 presents tortuosity behavior for the same electrode composition (90 wt. % active material) at different degrees of compression (directly related to calendaring pressures^{44}). The mappings are constructed (Figs. 3 and 4) from the electrode compressed to 600 bars pressure. For smaller calendaring pressures, the same mapping provides fairly accurate predictions, in turn, justifying minimal changes in the pore morphology. On the other hand, at a much higher pressure (2000 bars), the predictions are poorer as evident from the data spread in Fig. 5. Such a trend proves the hypothesized regimes for microstructural changes upon calendaring:

For lower calendaring pressures, electrode compaction is due to particle relocation (no qualitative change in the pore network).

Above a threshold calendaring pressure, particle fracturing takes place that considerably alters the pore shapes and in turn pore network undergoes qualitative, i.e., changes in pore morphologies.

## PERTINENCE OF PORE NETWORK BASED TORTUOSITY ESTIMATION

Superficially, two competing approaches should be compared in terms of accuracy and cost (in the present context, cost translates to computational demand). As shown earlier, the accuracy of the proposed approach is quite reasonable. Figure 6(a) provides a cost comparison [refer to workflow in Fig. 6(b) for complementary details]. Consider an electrode with *M*_{total} subvolumes. The computations scale linearly with the number of subvolumes (a porosity based correlation does not guarantee accurate predictions, e.g., it does not resolve subvolumes having similar porosities but pore network differences, consequently, requiring explicit calculations). On the other hand, the pore network based approach only relies on explicit calculations to form the learning set. For these *M* structures (in the learning set), additional calculations are required for the two-point correlations. Hence, if only a small number of structures are to be analyzed, DNS is a straightforward approach. On the other hand, for a large dataset (with similar pore network characteristics) such that *M* ≪ *M*_{total}, the pore network approach is computationally superior, as, for the subsequent (*M*_{total} − *M*) structures, only two-point correlations are to be computed and then the pore network → tortuosity mapping provides an accurate estimation of the transport distortions. The pore network based approach provides certain advantages:

Unlike porosity correlations, the proposed approach deciphers pore network differences among structures with similar porosities

Once the learning step is carried out, it provides a concise treatment of more electrodes with similar pore networks, for example, multiple samples prepared under identical conditions

Since directional information is consistently characterized when the structure is deconstructed into its essential attributes [Eq. (14)], treatment of anisotropy is unambiguous

Given the sensitivity to pore morphology variations, it can serve as a measure to discern nontrivial network changes

On the other hand, porosity based correlations do retain useful aspects. First and foremost, they provide mean trends that help ascertain bounds on the physicochemical phenomenon under study.^{71–78} Second, they do not require explicit 3D structures for calculations. Despite advances in imaging, large electrode structures are luxury and not always available and/or necessary. The two approaches are rather complementary, with pore network philosophy being relevant for large 3D structures and spatial variabilities.

## BROADER CONSIDERATIONS

Fluid-filled pore networks host many different transport processes such as solute diffusion, ionic conduction, and fluid flow. The specific network arrangement essentially offers additional resistance to these processes and is characterized differently for different processes, leading to definitions of tortuosity, permeability, etc.^{19} Since each of these fundamentally relates to the pore network, the aforementioned approach can be equivalently extended to map other properties. For example, if the interest were to comprehend permeability, the pore network attributes (as quantified through the two-point correlation) are to be correlated with the permeabilities (as computed from Stokes flow in porous networks^{19,79,80}). Equivalently, Eq. (15) is replaced with a relationship, *κ* = f(*y*^{1}, *y*^{2}), where *κ* is the permeability in a given direction.

Alternatively, ion conduction is governed by ∇ · (*σ*∇*V*) = 0, where *σ* is the conductivity and *V* is the electric potential. When abstracting the role of pore network on ionic conduction, corresponding interface and boundary conditions are also mathematically identical to pore-scale solute diffusion ( Appendix A). In turn, the ratios *σ*^{eff}/*σ* and *D*^{eff}/*D* are related *σ*^{eff}/*σ* = *D*^{eff}/*D* = *ε*/*τ*. In other words, the tortuosity relation [Eq. (15)] simultaneously characterizes the distortion to diffusion and migration modes of species transport in fluids. Equivalently, the electrochemical literature employs identical tortuosity for both ion diffusion and ion migration.^{1,3,12–14,4–11}

Herein, we discuss the two-point correlation function, $\psi r\u2192pore,pore$, that describes the pore network given its direct relevance to electrochemical transport in electrodes. Other two-point correlation functions also provide useful information. For example, $\psi r\u2192solid,solid$ abstracts the electronic transport in the solid phase. On the other hand, if the solid phase were made up of two different materials, the REV scale conduction jointly depends on three functions $\psi r\u2192solid1,solid1$, $\psi r\u2192solid2,solid2$, and $\psi r\u2192solid1,solid2$. The first two account for conduction paths in each of the materials, while the latter describes the cross-conduction. Such extensions are trivial given that these different transport processes are defined by mathematically similar (elliptic partial differential equation) governing equations.

Most studies of coupled transport processes in porous networks^{71–76,81–85} still rely on simplified effective property relations, e.g., the Bruggeman^{43} relation for tortuosity and the Kozeny-Carman^{86,87} equation for permeability. Accounting for different pore morphologies or anisotropy is computationally expensive and prohibits more realistic predictions.^{20,40–42,88–96} The appeal of the proposed approach lies in its sensitivity to essential pore network features (size distribution, morphology, connectivity, anisotropy) with minimal computational overhead. Additionally, the framework is equally applicable to different transport modes (hence different effective properties) encountered for fluids in porous networks.

## CONCLUDING REMARKS

The functionality of high energy high power electrodes is limited by the transport inefficiency of convoluted pore networks formed due to their porous nature. Tortuosity is a measure of the increased resistance of such networks. It is conventionally expressed in terms of porosity (a measure of pore volume in electrodes) that fails to explicitly account for pore system features such as size distribution, connectivity, and morphology. Such attributes are particularly relevant to electrodes exhibiting spatial variabilities such as inhomogeneity and anisotropy.

We show that the two-point correlation function proffers the mathematical framework to quantitatively represent these aspects (i.e., deconstructing the pore network). Eigen decomposition groups these details into manageable descriptors that map to transport distortion (characterized by tortuosity). Unlike porosity based correlations which require a large dataset for correlation development, the present approach leverages the explicit interpretation of the networks to map tortuosity in terms of essential attributes (in turn requiring a smaller dataset). In addition to the computational efficiency, this elemental difference yields a naturalized treatment for the three-dimensional qualities and is sensitive to differences in equiporosity structures belonging to the same microstructure class. The tortuosity mapping also possesses investigative merits as it can discern pore morphology changes, for example, particle fracturing upon calendaring.

## ACKNOWLEDGMENTS

This work was supported in part by the National Science Foundation (Grant No. 1805656) and Purdue University College of Engineering Seed Funding on Data Science.

### APPENDIX A: DIRECT NUMERICAL SIMULATION (DNS) FOR ESTIMATING TORTUOSITY

The electrochemical devices employ flooded electrolytes, with a well-mixed (i.e., distributed throughout) salt. Hence, the entire pore network partakes in the transport. Such a transport response is very different from point injection of a tracer, e.g., during angiogram in medical imaging. The latter represents a situation where the affected pore network evolves in time.

Given such dynamical differences, transport distortion of electrode pore networks is conventionally interpreted from steady state diffusion^{1,6,61,97} (large time, *t* → ∞, limit), while the early time transients are more relevant for biological networks.^{62,64}

To characterize transport distortion (i.e., tortuosity) in a microstructure, concentration balance is solved in the pore phase with fixed boundary conditions to cause macroscopically unidirectional solute flux^{33,37,39} (Fig. 7). The discrete representation forms a system of linear governing equations and solved using matrix solvers.^{65,98–106} The simulated concentration field is analyzed to estimate the diffusive flux. The presence of impervious solid matter introduces additional transport resistance, and equivalently the diffusive flux reduces as compared to its clear medium limit. A comparison of the two fluxes quantifies the transport distortion. Consider pore network transport in the x-direction. A suitable set of boundary conditions fix two different concentrations at *x* = 0 and *x* = *l*_{x} faces, with no flux through the lateral faces. If the domain entirely consisted of electrolyte, the steady state solute flux is $Jclear=\u2212DCx=lx\u2212Cx=0/lx$ where *D* is solute diffusivity in the medium. In the presence of a solid skeleton, the flux reduces to *J*_{porous} = *J*_{clear} · *ε*/*τ* which offers a consistent interpretation of directional tortuosity. Note that given the steady state nature of the concentration balance, flux is constant at each cross-sectional plane perpendicular to the flow, for example,

Similar identities follow for characterization in other coordinate directions.

Often tortuosity intuition is developed as a ratio of diffusion path to the shortest distance between two points.^{107} Such an interpretation, however, loses its veracity in structures with highly convoluted pore network as is found in densely packed electrodes or the ones with aspherical particulate matter.

### APPENDIX B: COMPUTING TWO-POINT CORRELATION FUNCTION USING FAST FOURIER TRANSFORM

Based on the definitions of the correlation functions [Eqs. (2) and (3)], it is straight forward constructing $\varphi x\u2192|pore$ from the material phase distribution in a given microstructure. For voxelated data with *N*_{x} × *N*_{y} × *N*_{z} entries, it takes $ONx\u22c5Ny\u22c5Nz$ operations to compute the one-point correlation function [Eq. (2)]. However, for subsequent higher order correlations, the computing operations grow rapidly. In general, it takes $ONx\u22c5Ny\u22c5Nzn$ operations to compute an *n*–point correlation, directly from the definition. The Fast Fourier Transform (FFT^{108,109}) provides a faster means to evaluate the higher order correlations. The following discussion first sketches the necessary steps on a 1D space and subsequently provides relations for 3D computations.

Consider a one-dimensional microstructure with *N*_{x} voxels. Corresponding discrete one-point correlation, $\varphi x\u2192|h$, accepts discrete $x\u2192$ vectors. There are *N*_{x} discrete $x\u2192$ vectors, which can be assigned an index *p* ∈ [0, (*N*_{x} − 1)]. Equivalently, discrete one-point correlation can be identified as $\varphi ph\u2261\varphi x\u2192p|h$. In the same spirit, the one-dimensional discrete two-point correlation function becomes [upon conversion from continuous to discrete space, the integration operation in Eq. (3) changes to a summation],

where the discrete separation vector is referred to via index, *r* ∈ [0, (*N*_{x} − 1)]. The periodic boundary conditions imply that $\u2200p+r\u2265Nx:\varphi p+rg=\varphi p+r\u2212Nxg$.

For a discrete function *f*_{p}, *p* ∈ [0, (*N*_{x} − 1)], its FFT is defined as

where $F$ signifies the transform and *F*_{q}, *q* ∈ [0, (*N*_{x} − 1)] is the transformed function with *q* as the corresponding independent variable. Here, $i=\u22121$. If *p* represents time, then *q* identifies frequencies. For the present application, *p* refers to spatial locations, and equivalently *q* leads to wavelengths (*q* ≡ spatial frequency = 1/wavelength). The inverse FFT is defined as

Strictly speaking, relations (B2) and (B3) define Discrete Fourier Transform (DFT), and FFT is an algorithm for the fast estimation of DFT. However, these terms are often used interchangeably. For consistency, small case refers to functions in spatial coordinates, while title case denotes transformed function in the wavelength space.

The transformed correlation functions can be expressed as follows:

Here, *q*, *s* are wavelengths for *p*, *r*, respectively. Rearranging terms in $\Psi sh,g$,

Let *t* = *p* + *r* ⇒ *r* = *t* − *p* and the limits for inner summation change to *t* = *p* and *t* = (*N*_{x} − 1) + *p*,

The periodic nature of the structure simplifies the limits in the second summation,

Substituting for $\Phi qh$ and its complex conjugate $\Phi *qh$,

Using the definition for the inverse FFT [Eq. (B3)],

Expression (B8) is the identity helpful for a quicker estimation of two-point correlations. If the calculations are to be carried out in 3D, three-dimensional, FFT is to be employed. In turn, the frequency space contains 3D frequency vectors, rather than just frequency values. Each component of such a frequency vector identifies the frequency distribution in respective coordinate directions. The identity (B8) holds true and will involve a 3D inverse FFT instead [Eq. (B9)],

## REFERENCES

_{2}cathode

_{1/3}Mn

_{1/3}Co

_{1/3}]O

_{2}cathode

_{2}

_{2}sequestration