Quantitative Relationships Between Pore Tortuosity, Pore Topology, and Solid Particle Morphology Using a Novel Discrete Particle Size Algorithm

To sustain the continuous high-rate charge current required for fast charging of electric vehicle batteries, the ionic effective diffusion coef ﬁ cient of the electrodes must be high enough to avoid the electrode being transport limited. Tortuosity factor and porosity are the two microstructure parameters that control this effective diffusion coef ﬁ cient. While different methods exist to experimentally measure or calculate the tortuosity factor, no generic relationship between tortuosity and microstructure presently exists that is applicable across a large variety of electrode microstructures and porosities. Indeed, most relationships are microstructure speci ﬁ c. In this work, generic relationships are established using only geometrically de ﬁ ned metrics that can thus be used to design thick electrodes suitable for fast charging. To achieve this objective, an original, discrete particle-size algorithm is introduced and used to identify and segment particles across a set of 19 various electrode microstructures (nickel-manganese-cobalt [NMC] and graphite) obtained from X-ray computed tomography (CT) to quantify parameters such as porosity, particle elongation, sinuosity, and constriction, which in ﬂ uence the effective diffusion coef ﬁ cient. Compared to the widely used watershed method, the new algorithm shows less over-segmentation. Particle size obtained with different numerical methods is also compared. Lastly, microstructure-tortuosity relationship and particle size and morphology analysis methods are reviewed.

Lithium-ion battery (LIB) electrodes have complex porous microstructures that are linked to their macroscopic performance.Typical LIB macroscale electrochemical models [1][2][3][4][5][6][7][8][9][10] are based on porous electrode theory [11][12][13][14][15] and thus abstract microstructural heterogeneity of composite electrodes using effective macroscopic properties.Four microstructure parameters are usually considered: 2 (i) the porosity e; (ii) the specific surface area S , p used to scale the reaction current; (iii) the particle diameter d , 50 which controls the maximum diffusion length of the solid phase diffusion; and (iv) the tortuosity factor t , i which denotes the effect of the convoluted, tortuous, path of the pores that hinders the Li-ion diffusion along the direction i. Parameters and acronyms are referenced in Table I.Achieving an accurate estimation of these parameters is essential to provide accurate predictions in agreement with experimental data. 16,17Fast ion transport along the electrode thickness (i.e., low tortuosity factor for a given porosity, cf.8][19][20][21] Tortuosity factor is deduced from the normalized effective diffusion coefficient (defined as the ratio between the effective diffusion coefficient D eff i , along the direction i and the bulk, dense, diffusion coefficient D bulk ) and the porosity (cf.Eq. 1). 22Previous work has demonstrated tortuosity factor is a key parameter to enable fast charging for thick electrodes. 17It is then relevant for automotive applications to design LIB electrodes with the lower tortuosity factors to facilitate fast charging.
While tortuosity factor can be experimentally determined, 16,23 its value can also be obtained through a homogenization calculation performed on a three-dimensional volume reconstruction of the electrode microstructure obtained either by X-ray computed tomography (CT) or focused ion beam scanning electron microscopy (FIB-SEM).7][28][29] However, the homogenization calculation is an indirect determination of the tortuosity factor, as it computes the normalized effective diffusion coefficient and then only deduces the tortuosity factor using Eq. 1.This method does not give any insights about the intrinsic relationship between the microstructure morphology and the tortuosity factor, which is valuable information for designing tailored lowtortuosity electrodes.In the battery community, the standard microstructure-tortuosity relationship is the Bruggeman's analytical law 30,31 (cf.Eq. 2) and its empirical counterpart, the generalized Archie's relationship 32,33 (cf.Eq. 3).The limits of the Bruggeman's analytical law have been extensively discussed in the literature.Its validity is restricted to unisize isotropic (i.e., spherical) particles in an electrode with moderate to high porosity. 16,31,33,34The Bruggeman law fails to capture the tortuosity factor of low-porosity and/or non-spherical particles such as anisotropic flake-like graphite. 16The empirical Archie's relationship is more flexible, as it introduces two microstructure-dependent parameters, g i and a , i with i being the investigated direction.Although it has demonstrated accurate correlations, each new microstructure type requires a new set of parameters g i and a , i which drastically limits its predictive capability for new microstructures. 16,26Furthermore, g i and a i do not have a geometric definition that could be used to design a lowtortuosity factor electrode.6][37] Many of these relationships tend to disagree for the low porosity range. 36,37There is then an identified need to establish a microstructure-tortuosity relationship based not only on porosity, but also on additional geometrically defined parameters, in order to provide a more generic relationship that will apply on a wider range of different microstructures, eventually producing valuable design recommendations for low-tortuosity and low-porosity (hence high-rate) electrode materials.
In this work, two issues of microstructure-tortuosity correlation are addressed: lack of genericity, which limits prediction solely to microstructures for which correlation parameters have been fitted, and use of nongeometrically defined parameters, which prevents design recommendations.To achieve this, pore topology and particle morphology of various electrodes are quantified through a set of geometrically defined metrics using a novel particle identification algorithm, and then correlated with the tortuosity factor and the normalized effective diffusion coefficient.The new established correlations stand for a wide range of microstructures (NMC and graphite) with different porosities and particle shape and can be then considered generic, in contrast with the empirical Archie's relationship, which is microstructure specific.The tortuosity factor is then defined directly (i.e., as a function of defined geometric parameters) rather than indirectly (cf.Eq. 1), which enables design recommendations.An open-source electrode library is used (cf.paragraph section Electrode library, over-segmentation, and particle identification) so that other groups can compare identification results with their own methods.The next two sections review microstructure-tortuosity relationships and particle size and morphology analysis methods reported in the literature, as well as their limits.

Review of Microstructure-Tortuosity Relationship
Tortuosity factor encompasses all the contributions of the microstructure, not only the porosity, in the normalized effective diffusion coefficient.Several authors [38][39][40][41][42][43] have attempted to give a geometrical definition of the factor of tortuosity, thus a direct definition instead of the indirect definition given by Eq. 1. Two geometric properties are commonly considered in addition to the porosity and are illustrated in Fig. 1 for different geometries.
First, the sinuosity of the microstructure is quantified through the geometric tortuosity t , geo i , usually defined as the average shortest distance within the pores from the face normal to the direction i of the microstructure domain to the opposite face á ñ l i normalized with the domain's length along the same direction L i (cf.Eq. 4 and Fig. 1a).Geometric tortuosity is an oriented parameter with a direction: t + geo i , and t - geo i , are the geometric tortuosity along direction i calculated from bottom to top and from top to bottom, respectively.Domains with a small field of view, or case-study geometries, may have a different number of paths between opposite faces from bottom to top and from top to bottom, which can result in t t ¹ + geo i geo i , , (cf.Fig. 1a).5][46] Although, the skeletonization algorithm will induce a bias in the analysis as the extracted topology may be different depending on what algorithm has been used. 40urthermore, other criteria than the shortest path can be used, such as highest permeability path for maximal flow. 47= á ñ l L 4 Second, the variation of the section area along the pathways [38][39][40][41][42][43] (cf.Fig. 1b) is described with the constriction factor b, the constrictivity d, and the geometric constrictivity d .
geo The constriction factor, initially introduced by Petersen, 38 is defined as the ratio between the large section A max (the bulge) of radius r max and the small section A min (the bottleneck) of radius r min of the pore domain embedded in periodically ordered spherical particles (cf.Eq. 5).Because the particles are aligned, his case study corresponds to a unit geometric tortuosity.Petersen 38 analytically demonstrates that the diffusion within this simplified pore medium, which exhibits periodic constrictions, is controlled by the constriction factor.Michaels 39 demonstrated later that the variation of the section area has a much higher impact on the effective diffusion than the distance between the bottleneck and bulges.Holzer et al. 40 have fitted the effective normalized diffusion coefficient according to Petersen's result (cf.Eq. 6).However, the simple geometry used in Petersen's Journal of The Electrochemical Society, 2020 167 100513 study limits the use of this relationship for actual, more complex electrode microstructure geometries.Holzer et al. 40 have then proposed to characterize the constriction factor directly from actual microstructure tomography through a particle-size distribution analysis that combines a continuous particle-size distribution (c-PSD) algorithm to get the average largest section á ñ A max and a mercury intrusion porosimetry (MIP) simulation to calculate the average smallest section á ñ A min (cf.Eq. 7).Indeed, size distributions calculated with MIP simulation tend to overrepresent the bottlenecks, thus providing an estimation of the minimum radius. 40owever, the choice of using c-PSD to determine the maximum section is disputable, as c-PSDs are known to underestimate the actual particle size, especially for non-spherical particles. 48,49is used to decompose the tortuosity factor between the geometric tortuosity and the constrictivity (cf.Eq. 8). 29,42As defined, the constrictivity determination is indirect: it is deduced from the geometrical tortuosity, the porosity, and the normalized effective diffusion coefficient.To achieve a direct determination (i.e., based only on geometrically defined parameters), the constrictivity is correlated with the geometrically defined constriction factor (cf. Eq. 9).Because the constrictivity is an oriented parameter, it implies that the definition of the constriction factor should also be anisotropic.The constrictivity obtained from such a correlation is named geometric constrictivity and noted d .
geo For instance, Holzer et al. 40 have established an empirical relationship for porous olivine diaphragms (cf.Eq. 10).The normalized effective diffusion coefficient is then a function of geometrically defined parameters: the porosity, the geometric tortuosity, and the constriction factor through its correlation with the geometric constrictivity.Another interest of this approach consists in its orthogonality: while the tortuosity factor is a combination of geometric tortuosity and geometric constrictivity for complex electrode microstructure, some simple geometries exhibit either no geometric tortuosity or no geometric constrictivity (cf.In addition to these geometric parameters, several authors have pointed out the role of the particle alignment for non-spherical particles on the tortuosity factor. 16,50,51Pore embedded between ellipsoid-like particles exhibits a much lower tortuosity factor along the direction aligned with the ellipsoids long diameter compared with the direction orthogonal with it. 16,50,51Indeed, particle alignment impacts both geometric tortuosity and constriction factor, as illustrated in Fig. 2.This illustration suggests that constriction factor should be defined as an oriented parameter.These three geometric parameters are not independent, which implies that a microstructuretortuosity relationship can include all of them or only a subset of them.The definition of the constriction factor, geometric tortuosity, and particle elongation used in this work is given in the section Particle morphology and pore topology parameter determination. Whatever the microstructure-tortuosity factor relationship is, the tortuosity factor should always converge towards 1 as the porosity is increasing.This implies that the contributions or weights of these additional geometric parameters (geometric tortuosity, constriction factor, and particle alignment/elongation) is decreasing with the porosity, i.e., their values are less important.Contrariwise, for the low porosity range, not taking them into account can lead to inaccurate tortuosity predictions.This explains the wide prediction spread in the low porosity range obtained with porosity-tortuosity factor relationships reported in the literature. 36,37Evaluating particle morphology (to deduce particle elongation and alignment) and pore topology (to deduce geometric tortuosity and constriction factor) is then particularly essential for establishing microstructure-tortuosity factor relationships valid in the low-porosity range for a wide type of electrode materials that could be used to provide practical design recommendations for LIB electrode design for electric vehicle applications.Review on particle size and morphology analysis methods.-Particleshape and size can be quantified through optical observation of the powder.However, particle spatial distribution and alignment within the electrode microstructure, which controls both the constriction factor and the geometric tortuosity (cf.Fig. 2), requires more in-depth observation techniques such as CT and FIB-SEM.Metrics to characterize particle size and morphology are well defined. 52However, three-dimensional volume reconstructions obtained with these techniques usually do not show distinct particles, but rather a unique connected cluster (cf.Fig. 3 and Refs.16, 49).While some methods allow for calculating particle size without identifying distinctively particles (c-PSD, 48,53,54 MIP simulation, 53 covariance, 16,25,49,55 chord length, 56 and intercept method 57,58 ), they provide no or limited (e.g., size anisotropy 16 ) information on the particle morphology.Therefore, a proper particle identification or segmentation step is required prior to the morphology analysis.Particle identification is challenging as no unique definition exists for what a particle is for a connected cluster.Furthermore, LIB electrode material particles exhibit different geometric features (e.g., cracks) and shapes (e.g., quasi-spherical, flake-like) that prevent the use of a simple geometric definition (cf.Fig. 3).In addition, particles with significant cracks could be considered as distinct particles.Also, the bottleneck critical size used to visually distinguish (respectively, merge) adjacent geometric features into two particles (respectively, one particle) is highly subjective, especially for microstructures that exhibit a wide particle size distribution (cf.Fig. 3).Therefore, a particle identification numerical method should be parameter-free to avoid subjectivity issues.
][61][62][63][64][65][66][67] This family of algorithms is based on the same approach: simulating a user-defined phenomenon (e.g., erosion, 64 flooding or immersion, [59][60][61][62]67 drop of water principle, 65 or seed growing 66 ) within the microstructure segmented image, and then identifying the particles based on the simulation result interpretation. Thse algorithms identify particle bulk (called seeds 66 or catchment basins 65 depending on the employed method) and particle connections (called dividing, or watershed, lines 65 or concave bottle necks 48 ) to label distinct particles.52 Most d-PSD methods are based on the Euclidean distance map (EDM), which is the minimum distance to the domain's boundary 64,[68][69][70]. A isual representation is given in a previous work.49 These algorithms often consider the EDM as a topographical relief 60 and are based upon the mathematical morphology notion of watershed line.Watershed line is an essential tool developed by Beucher and Lantuejoul, 71 initially for segmenting images, which delimit the zone of influence of each phase or particle.A large number of methods have been developed to calculate them. 65 Onof the most common watershed algorithms consists in marking the opposite EDM minima as water sources and then progressively flooding the domain, considering the opposite EDM as a topography map, each source point being the bottom of a distinct lake (i.e., the immersion of flooding approach).[59][60][61][62] When two lakes, or catchment basins, meet each other, an imaginary dam is placed on their boundary to prevent distinct lakes from merging.60 The union of all dams are the Figure 3. Two-dimensional slices extracted from CT volumes for (a) graphitic anode and (b) nickel-manganese-cobalt-oxide (NMC) cathode. (a Solid red lines represent hand-made particle identification based on the assumption that particle external contour is roughly spherical, while dash lines are alternative plausible particle identification without any shape presupposition. Inthis example, the pore 3D domain is fully connected.(b) NMC quasi-spherical particles appear to be connected, with bottlenecks of various size (blue lines) that make their identification subjective.
watershed lines, which delimit the influence zone of each catchment basin. 59,60The iterative process stops when all voxels of the phase have been attributed to a lake, which are then interpreted as distinct particles.Illustrations and results of this method can be found in a previous work. 49Immersion-type algorithms are still developed, with an emphasis on parallel computation in order to increase the size of analyzable domain in a decent time. 67,72Alternatively, Cousty and Bertrand 65 have formalized the "drop of water" principle, for which a drop of water falling on a topographic surface (represented by the EDM opposite) follows a descending path until it reaches a stable position, the local minimum of the EDM opposite.In this approach, the watershed lines are the separating lines of the domain of attraction of each catchment basin.Saadatfar et al. 66 have developed another EDM-based method, which considers the local EDM maxima as distinct seeds, one for each individual particle.A region growing algorithm is then used until all voxels of the phase have been attached to the initial seeds' voxels. 66Lastly, Münch et al. 64 have developed a splitting recognition algorithm, based upon a stepwise erosion process of the binary image, that will progressively disconnect touching microstructure features and then identify them.This algorithm also uses the EDM as eroded regions are identified based on their distance from the domain's boundary.
Due to the absence of any particle morphology assumption, d-PSD methods are expected to provide particle sizes closer to the actual sizes, compared with those obtained with a spherical assumption (continuous particle-size distribution, noted c-PSD), 48,49,73 especially if particles are actually non-spherical.However, because watershed methods are gradient based, they usually suffer from over-segmentation or over-identification due to the presence of concavities on the particle surfaces (e.g.non-smooth interfaces), which induce a noisy gradient image and generate a lot of irrelevant, oddly shaped small particles. 59,60,62,64,74Correct watersheds lines can be then lost in a mass of irrelevant ones, even if the gradient image has been previously filtered. 59Several authors have proposed methods to limit over-segmentation. 49,64,73Adjacent catchment basins can be merged according to some criteria, or preexisting knowledge of the geometric features to be identified can be used to perform an initial marking step. 59Münch et al. 64 use an additional parameter to control the magnitude of the erosion process, which can be tuned to reduce oversegmentation.Usseglio et al. 49 proposed to discard particles assigned with a diameter lower than the one obtained with a c-PSD algorithm and reassign their voxels to adjacent particles, assuming that the smallest possible geometric feature that can be considered as a particle is the largest sphere and everything below is an artifact (i.e., assuming c-PSD gives the particle size lower bound).Except for case study geometries (e.g., spheres 66 ) with no surface roughness, the outcome of a d-PSD algorithm is difficult to predict due to the absence of a morphology constraint. 64Therefore, visual inspection is required to validate the identification results.Furthermore, some methods use parameters either to control the magnitude of the physical mechanism simulated 64 or to filter the gradient image, such as by merging nearby local minima based on a user-defined critical distance. 70However, the choice of these user-defined parameters reduces the result's unicity.Lastly, d-PSD methods suffer from a border effect that underestimates the size of the truncated particles located at the domain's boundaries.Münch et al. 64 proposed two methods to correct this edge effect: (i) an unspecific particle volume correction factor applied on all particles based on the fact that the probability of larger particles being truncated is higher than for smaller particles, and (ii) a specific particle volume correction factor of 2 n applied only on particles truncated by n adjacent planes, assuming particles are randomly oriented.
Aim and organization of the article.-Tortuosityfactors have been determined in previous work 16,49 with an indirect approach (cf.Eq. 1), thus without establishing a direct microstructuretortuosity correlation valid for a wide range of different microstructures and involving only geometrically defined parameters.Such relationships are essential to provide design recommendations for low-tortuosity and low-porosity electrode materials.Particle morphology and pore topology analysis are required to provide the parameters of the desired microstructure-tortuosity relationships, and a particle identification step is needed prior to the particle morphology analysis of typical LIB electrode reconstructed volumes.Section Particle size and morphology original algorithms introduces two novel algorithms used to calculate particle size and to identify particles.Section Particle morphology and pore topology parameter determination defines the morphology and topology parameters that are extracted from the particle identification.Section Method applied with reference geometries illustrates particle identification performed on case-study geometries, with a reference watershed (immersion approach) method and the novel algorithm.Section Methods applied with electrode library shows results obtained on tomography-based LIB electrode reconstructed volumes.The electrode open-source library is described in paragraph section Electrode library, over-segmentation, and particle identification.Particle sizes obtained with different methods are compared in paragraph section Particle size comparison, particle and pore morphology results are shown in paragraph section Particle morphology, and pore topology (i.e., geometric tortuosity and constriction factor) in paragraph section Pore topology.Finally, paragraph section Pore tortuosity factor correlation with microstructure geometry shows microstructure-tortuosity correlation established based on the metrics results discussed in the two previous paragraphs section.Discussion and conclusion sections conclude the article.

Particle Size and Morphology Original Algorithms
In this work, particle sizes are evaluated through a set of different numerical methods: covariance, from specific surface area, c-PSD, and watershed (immersion), with algorithms detailed in previous work. 49In addition, two new methods are introduced and are detailed in the paragraph below: Euclidean distance map fitting (EDMF) and pseudo-Coulomb repulsive field (PCRF).Table II summarizes the pros, cons, and assumptions of each method.
Euclidean distance map fitting method.-Euclideandistance map (EDM), also called Euclidean distance transform, is defined for each point x of the investigated phase V as the shortest distance from this point to the phase boundary ¶V (cf.Eq. 11).Previous work shows an illustration of such a map. 49The EDM is calculated with the MATLAB built-in function bwdist based on a linear time algorithm. 68The EDM cumulative function C R d , EDM sphere ( ) of a sphere of radius R is defined as the volume sum of all the infinitesimally thin layers distant from the sphere boundary with a distance superior or equal to d, normalized with the sphere volume (cf.Eq. 12).Applied to an arbitrary microstructure phase V discretized in N voxels x i of volume W, the EDM cumulative function C d EDM ( ) is the volume sum of all the voxels distant from the phase boundary with a distance superior or equal to d, normalized with the phase volume (cf.Eq. 13).
The Euclidean distance map fitting (EDMF) original method consists in finding R so that the difference between C R d , EDM sphere ( ) and C d EDM ( ) is minimal.An example is shown in Fig. 17c.If the microstructure is constituted of nonoverlapping spheres with a unique radius R, then this method provides the exact radius.Even though the Journal of The Electrochemical Society, 2020 167 100513 fitting step relies on a spherical assumption, the calculation of C EDM is assumption-free.Furthermore, this method is parameter-free.
The general approach used in this original method is to first choose a function F known as a function of a parameter P for a reference geometry G, then to fit F with P to match the function calculated on an arbitrary domain so that the fitted value of P would provide an approximation of P for this domain.Here, F is the Euclidean distance map cumulative function, P the characteristic size, and G a sphere.This approach has been used in a previous work, 49 with F being the covariance function and G an infinite medium constituted of spherical particles with a unique diameter.It is believed by the authors that other doublets (F, G) exist to determine the characteristic size of an arbitrary domain.
Discrete particle size distribution: pseudo-Coulomb repulsive field method.-Theoriginal method presented in this section is a discrete particle size-distribution (d-PSD) algorithm that provides particle size but also identifies distinctly particles that is valuable to quantify particle morphology and pore topology metrics defined in section Particle morphology and pore topology parameter determination used to correlate tortuosity factor.The new method is expected to be more robust than the reference watershed method for reasons detailed in this section.The novel method is described below.
Boundaries of the investigated phase are considered as fixed walls, positively charged.Imaginary positively charged particles are then dropped or generated between these fixed walls (i.e., within the phase volume).We will call these imaginary charged particles "charged points" to avoid confusion with active material or lithiumion particles.Because both points and walls are charged with the same polarity, they repulse each other according to a pseudo-Coulomb's law (cf.Eq. 14).Since walls are fixed, only the charged points will move from their initial position.Points are pushed away from the phase boundary and converge toward particle centers.Eventually, the sum of the trajectories that converge to the same Table I.Nomenclature. 1 As defined in this work.d d d : : ( ) This is valuable to evaluate spatial heterogeneities, edge effect, and/or graded electrodes.b) Diameter is known along the three spatial directions.c) Particles are identified distinctly from the analysis of the fully connected domain.It can be used to quantify particle morphology and pore topology through a domain's skeletonization.d) c-PSD, covariance, and watershed methods used in this work are detailed in previous work. 49e) Spherical particle assumption is used only to fit the algorithm result.The assumption is then less strong than for the c-PSD method where spherical assumption is used to directly calculate the algorithm result.
Journal of The Electrochemical Society, 2020 167 100513 stable position defines a unique particle.The method is illustrated in Fig. 4. We call this method in this work the pseudo-Coulomb repulsive field (PCRF) method.
The pseudo-Coulomb field F ¯is calculated as described here (cf.Fig. 4a).First, voxels that belong to the complementary phase but are adjacent to the investigated phase are identified and noted w, ¯and represent the fixed wall.Then, for each voxel x i that belong to the investigated phase, the local value of F x i ¯( ) is calculated (cf.Eq. 14), with n as the number of wall voxels and l ij a line of sight binary value: 1 if a straight path exists between x i and w j and 0 otherwise (cf.Eq. 15).In the example shown in Fig. 4, the stable local minima of F   are not the two spheres' centers.The value of k is set to 2 by default to keep the analogy with Coulomb's inverse-square law, although different values can be tested (higher values decrease the contribution of distant walls).This is the only parameter of the method and its impact on the particle identification is discussed in sections Method applied with reference geometries and Methods applied with electrode library.
The second step consists of simulating the displacement of charged points within the microstructure according to the field F ¯to identify particles.A 3D array P is first initialized (zero) and will keep track of the particle identification process, with p as the number of particles identified (p is initialized at 0).An initial position x ¯is randomly chosen within the microstructure and assigned to The charged points are then moving stepwise according to the sign of F, ¯i.e., either positive or negative (cf.Fig. 4b), and not according directly with the gradient of F x ¯( ¯) (cf.Fig. 4c).Stable points between the two approaches, however, are identical.Three cases are then possible: x x ¯is calculated (cf.Fig. 5a).The complete trajectory t of the particles is saved in memory, as explained later.2. If the new position brings back the point to a previous location of the current trajectory or if the point displacement stalls (cf.Fig. 5b), i.e., respectively ¢ = + P x p 1 ( ) and ¢ = x x, ¯then the particle number is incremented = + p p 1. A complete stop corresponds to a stable or unstable position of F, ¯while a backtracking can happen if the local minimum actual position is between two voxel positions.The last two positions x ¯and ¢ x are also saved and correspond to the particle center of attraction.A new initial position x ¯is randomly chosen among the points not already assigned to P, and the algorithm continues until all voxels of the phase have been assigned to a particle number p.

The last case corresponds to
which means the point intersects another trajectory that belongs to a particle previously identified (cf.Fig. 5c).In such a case, all voxels of the current trajectory t are assigned to ¢ P x , ( ) i.e., = ¢ P t P x ( ) ( ) and the particle number is not incremented.Note that if we calculate the next displacements, the trajectory will overlap with the one just intersected.A particle is then defined as the sum of all the trajectory that converge to the same final location.A new initial position x ¯is randomly chosen among the points not already assigned to P, and the algorithm continues until all voxels of the phase have been assigned to a particle number p.
The method is domain agnostic, as it can be applied separately on the solid (cf.Fig. 4) or pore (cf.Fig. 6) phase identically depending on what information is wanted: solid particle morphology or pore topology, respectively (cf.section Particle morphology and pore topology parameter determination).To avoid confusion, the identified regions are named "particles" and "pore segments" when referring, respectively, to the solid domain and the pore domain.
Four post-processing steps are performed to improve the particle identification.First, particles with adjacent centers of attraction are merged.Second, it has been observed that in some cases, a chessboard pattern emerges at the intersection between two particles (cf.Fig. 6a).Such numerical artifact is derived from the discrete displacement step.Voxels that belong to such a pattern are identified and assigned to the particle with the nearest center of mass (center of mass, or centroids, are calculated for each particle, ignoring the voxels that belong to chessboard patterns).The operation is repeated until no chessboard pattern is detected.Third, particle continuity is checked.If they are found to be noncontiguous (in rare occasion, the chess-pattern removal algorithm induces such noncontinuity), the largest cluster is kept while all the other smaller clusters are assigned to the particle for which they share the largest interface.Lastly, onevoxel-sized particles are removed, and their unique voxel is assigned to the particle for which they share the largest interface.Some particle identification algorithms suffer from over-identification or over-segmentation, which leads to small irrelevant features. 59,60,62,64,74In the case of the watershed method, such artifact issue is exacerbated as the basin source locations are derived from the minimum distance from the boundary, thus a unique distance value.Because of this, multiple irrelevant basin sources may be incorrectly generated if the phase surface is rough.In this work, the local field value used to simulate the point displacement is based on multiple distance values (cf.Eq. 14).Because of this, surface roughness is expected to cause less over-segmentation and then the new method is expected to be more robust than the reference watershed method.Over-segmentation of both methods are compared in section section Method applied with reference geometries and section Methods applied with electrode library.In order to remove any over-segmentation, the same technique previously developed for immersion watershed 49 is applied.Briefly, it consists in reassigning voxels for which the diameter is lower than the diameter of the largest sphere that contains it without overlapping the complementary phase (cf.Fig. 6b).This over-segmentation correction method is not specific to the pseudo-Coulomb repulsive field method and can be used for any other d-PSD algorithms.In this work, it is also used to correct watershed oversegmentation results, as done previously. 49he general approach used in this original method is to simulate a mechanism M within an arbitrary domain, for which the result interpretation I allows to determine a parameter P. Here, M is dropping charged points between fixed walls charged with the same polarity, I is the points trajectory draw the particle, and P is the particle identification.Similarly, for watershed methods, M is flooding the domain with a liquid, I is the catchment basins linkages correspond to particle-to-particle connections, and P is also the particle identification.It is likely that other doublets (M, I) that could provide a good particle identification exist.

Particle Morphology and Pore Topology Parameter Determination
Once particles and pore segments are identified, whatever the method is, particle morphology and pore topology can be quantified based on a set of metrics. 52Indeed, d-PSD algorithms are used on the solid domain (identified regions are named "particles") to deduce solid particle equivalent diameter, elongation, and sphericity, and on the pore domain (identified regions are named "pore segments") to deduce pore constriction factor and pore geometric tortuosity, as described in this section.Note that all metrics defined below are domain agnostic.For instance, pore sphericity or diameters can be calculated using the same description (simply replacing "particle" with "pore segments" in the definition).
[( ) ] / / with V being the particle volume defined as the volume sum of the voxels assigned to the particle.Sizes of truncated particles at the field of view boundary are not corrected.Mean particle diameter d 50 or á ñ d is calculated according to a rule of mixture weighted by the particle volumes: The particle average volume á ñ V is then deduced: Particle elongation or aspect ratio.-Particleaspect ratio d d d : : / thus as normalized aspect ratio and not as ellipsoids long and short diameters.Average particle elongation is calculating according to a rule of mixture weighted by the particle volumes: Sphericity.-Sphericity y of a particle with a volume V and a surface area ¶V is defined as the surface area of an ideal sphere of volume V divided with the surface area of the particle ¶V: 16.Sphericity quantifies the roundness of a particle: for a sphere y = 1, otherwise y < 1 and is decreasing as the particle shape is getting different from this ideal shape.However, the surface area of Journal of The Electrochemical Society, 2020 167 100513 a sphere of volume V discretized in cubic voxels is six times the disc area, 49,75 then the corrected sphericity is defined as simplified in Eq. 17 / / representative of the average pore segment (cf.Fig. 7, top right).The area of the disc normal to the direction i and that intersects with the ellipsoid centroid, á ñ A i max , provides an estimation of the average maximum section area (the bulge) normal to this direction (cf.Eq. 18).The interface between adjacent pore segments are fitted with the equation of a three-dimensional plane, to deduce the interface normal n , ] N being the number of interfaces, and = n 1.
j   The interface area A j is also calculated (cf.Fig. 7, bottom right).Since the interface area is the sum of voxel faces, its value is overestimated and must be corrected with a factor of 2/3 using a spherical assumption. 49,75The average minimum section area (the bottleneck) normal to the direction i, á ñ A , i min, is then defined as the average interface area weighted with n e ., [ ] (cf.Fig. 7 bottom right and Eq.19).For instance, if the interface normal is almost parallel with e 1 and perpendicular with the two other directions on average (i.e., á ñ ñ 1 0 0 ¯[ ] ), then the bottleneck section area along direction 1 á ñ A min,1 is almost equal to the average interface area á ñ A 2 3 , / while the two other directions á ñ A min,2 and á ñ A min,3 are almost 0. lastly, the constriction factor b i along direction i is defined as the ratio of á ñ A i min, over á ñ A i max, (cf.Eq. 20).As defined, the construction factor is an oriented property, as suggested by Fig. 2.
. , Geometric tortuosity.-Geometrictortuosity t geo i , along direction i is also deduced from the pore domain identification.Once pore segments are identified, the pore domain is skeletonized using a graph analogy: graph nodes are the pore segments' centers of mass, while graph edges connect nodes that correspond to adjacent pore segments in contact (cf.Fig. 8).Edges are bidirectional (undirected graph) and weighted by the edge distance.To calculate t + geo i , or t -, geo i , sources nodes are all the nodes that correspond to pore segments in contact with the first (or last, respectively) slice of the 3D field of view normal with the direction i, while target nodes correspond to pore segments in contact with the last (or first, respectively) slice of the 3D field of view normal with the direction i.Then, the shortest distance + l i k , (respectively - l i k , ) from source node k to any target nodes are calculated with the Dijkstra algorithm 76,77 for all source node k (the MATLAB built-in function distances is used).Source nodes that do not have a path to a target node are ignored.They correspond to isolated (i.e., non-connected) pore domains.Distances from field-of-view boundaries to source and target nodes are added.The average shortest distance á ñ + l i (respectively á ñ l i ) is then deduced and normalized with the domain's length L i along the direction i to provide the geometric tortuosity t + [ ] [ ]

Method Applied with Reference Geometries
Reference geometries without surface roughness.-Predictingthe outcome of a discrete particle-size algorithm for complex topologies such as LIB electrode materials is a difficult exercise due to the absence of morphology constraint.It is then recommended to first test such algorithm on simple geometries to get an understanding of how the method works, and to evaluate the identification relevance, as the ground truth (i.e., the expected identification) is known most of the time.Identifications achieved with watershed (reference method) and PCRF methods are compared for a set of ten geometries (cf.Table III).These geometries have no surface roughness, except for the domain discretization in pixels.To compare their relative over-segmentation, results are shown without the over-segmentation correction method discussed in paragraph section Discrete particle size distribution: pseudo-Coulomb repulsive field method.Since watershed identification is based on the opposite EDM, the latter is also shown in the table.The pseudo-Coulomb field F ¯is also shown, as PCRF identification relies on it.
Both algorithms performed similarly for the first four cases #1-4 (except for a minor edge effect in #4): a unique sphere (#1) and its complementary volume (#2) and overlapping spheres (#3) and their Table III.Particle identification illustrated for case-study two-dimensional geometries compared between watershed immersion method (reference) and PCRF method (new method).Cases 1-4 provide similar identification, while Cases 5-10 are better with the new method.Over-segmentation is not corrected on purpose for the sake of the comparison.The number of identified particles is noted n.

Watershed method
Pseudo-Coulomb repulsive field method (PCRF) Journal of The Electrochemical Society, 2020 167 100513 complementary volume (#4).Cases 1 and 3 correspond typically to solid particle analysis, while Cases 2 and 4 better suit the pore topology analysis.Note that the bottleneck identified in Case 4 would provide the exact section area A min used to calculate the constriction factor b.For the remaining Cases 5-10, the new PCRF method achieved a better identification compared with the reference watershed method.Case 5 tests algorithms for spheres connected with a long channel: watershed failed to distinguish spheres from the channel, while PCRF identified the distinct geometries with minor errors.Note that the best algorithm to use for such spheres and cylindrical channels is actually the c-PSD method as it will provide two sizes: sphere diameter and cylindrical diameter.Cases 6-8 test algorithm responses for spheres with inclusions.Case 6, a sphere with a centered spherical inclusion, is a particularly tough case to resolve as the local maximum of the distance map is theoretically a continuous ring.There are two possible solutions to this problem: either considering the domain as a unique particle, or as the union of an infinity of one-dimensional particles radially oriented.While both algorithms fail to achieve either solution, over-segmentation is halved for the PCRF method.One reason that explains this numerical result is the numerical error: the absolute maximum of the calculated EDM is not a ring, but isolated pixels far from each other that prevent merging.Additionally, since the geometry is theoretically symmetric, any tiny asymmetry due to the pixel discretization would shift the identification process of both methods in one direction.Case 7 relaxes the problem difficulty by giving an orientation to the inclusion.PCRF correctly identified the four distinct regions, while the watershed method failed.Case 8 relaxes the difficulty by breaking the inclusion symmetry.Even though it can be subjective to decide what should be the perfect identification for such a case, the PCRF result is better than the watershed result based on the large odd (i.e., with a 90°angle) shapes present in watershed identification but absent from the PCRF result.The odd shapes identified with the watershed algorithm are similar to those found in the literature. 64Lastly, Cases 9 and 10 test algorithm performance for cracked particles.For both cases (an open crack in #9 and a closed crack in #10), the identification is better for the PCRF method.
The choice of the parameter k (cf.Eq. 14) did not change the identification of Cases 1-5 and 7, slightly impacts Cases 6 and 8, and significantly changes Cases 9 and 10, with a better identification obtained with k = 3. Identifications obtained with k = 1 are not shown as they were systematically less adequate, which disqualifies this value.Over-segmentation for these geometries can be reduced, and in some instances completely removed, using the correction discussed in paragraph section Discrete particle size distribution: pseudo-Coulomb repulsive field method (cf.Figs. 6 and 9).However, it is preferable that the initial over-segmentation is the lowest possible so that the result keeps a strong link with the identified catchment basins (watershed method) or center of attraction (PCRF method).For further validation, particle identifications obtained with the watershed method using the algorithm developed in a previous work 49 has been compared with the MATLAB built-in function watershed (cf.Fig. 9).Our in-house watershed algorithm shows less over-segmentation, but it is also slower.
Reference geometries with surface roughness.-Theorigin of over-segmentation is often explained in the literature 64 by noisy Euclidean distance map induced by small concavities at the particle surface (i.e., surface roughness).Geometries investigated in the previous paragraph have no surface roughness, expect for the domain's discretization in pixels.Therefore, their EDMs are not noisy.However, they still exhibit significant over-segmentation with the watershed method (Cases 6-10) and in one case (#6) for the PCRF method.The source of their over-segmentation lies in the geometry's theoretical symmetry and the theoretical absence of preferential direction (e.g., Case 6).Any asymmetry or preferential direction induced by the pixel discretization actually controls the identification process.Indeed, adding directions (Case 7) or breaking the symmetry (Case 8) reduces the over-segmentation.These ideal geometries are not representative of actual LIB electrode materials, for which particles and geometric features are never completely symmetrical.Therefore, testing or illustrating identification    algorithms on symmetric smooth geometries tests the algorithm over-segmentation sensitivity with symmetry, while over-segmentation for actual microstructure electrodes comes from surface roughness and not from symmetry, as there is no ideal symmetry in these materials.Evaluating identification algorithms on symmetric features without surface roughness is then a biased approach, as it is not representative of their performance on actual microstructures.Geometries are modified to investigate the impact of surface roughness on identification algorithms.A simple approach consists in randomly assigning the pixels at the domain boundary, either to the domain or to the background.This method is enough to test geometries with initially no surface roughness at all (cf.Fig. 10a).For such geometries, the watershed method is extremely sensitive with the surface roughness, generating a lot of irrelevant structures (i.e., over-segmentation), while far less artifacts are generated with the PCRF method, for the reason explained in paragraph section Discrete particle size distribution: pseudo-Coulomb repulsive field method.For rounded geometries (e.g., Case 3 in Table III), which already start with a non-smooth boundary due to the pixel discretization (i.e., with a 1-pixel-thick surface roughness), a thicker surface roughness is generated by upscaling the geometry obtained with the random pixel assignment method used previously (cf.Fig. 10b).A stronger surface roughness is obtained by upscaling an initial coarse geometry several times (cf.Fig. 10c).No over-segmentation is induced with either method.Surface roughness is then more an issue for channel-type geometries, for which the EDM local maximum is a continuous line, rather than for particle-type geometries, for which the EDM local maximum is a point at the particle center.

Methods Applied with Electrode Library
Electrode library, over-segmentation, and particle identification.-Electrodelibrary.-Athree-dimensional reconstructed volumes open-source data set from a previous work 16   volume fraction of the pore and the additive phase (i.e., carbon binder).Therefore, the binary domain contains: 0 (background) = pore + additives, and 1 = active material.This choice results from the inability to distinguish the additives from the pore phase with X-ray CT.The reader is invited to read 16 for details about the segmentation.In addition, five other graphite electrodes have been manufactured by CAMP: SLC1506T, SLC1506T2, SLC1520P (natural coated graphite), MCMB (artificial mesocarbon microbeads standard type G15), and BTR BFC-10 (artificial graphite Targray SPGPT805), with porosities of 0.374, 0.345, 0.351, 0.381, and 0.321, respectively, and thicknesses of 44, 70, 45, 48, and 49 μm, respectively.All five samples share the same weight ratio of the anodes previously investigated. 16Imaging and segmentation for these new electrodes have been done in a similar way compared with the previous data set. 16Open-source tomographic data are also available. 78Porosity of the total 19 electrodes investigated in this work are listed in Table V.
Over-segmentation.-Watershed and PCRF identifications are shown in Fig. 11.To provide a fair comparison, both results are displayed before over-segmentation correction.The watershed method Figure 12.Particle identification before and after over-segmentation correction illustrated for a A12 graphite microstructure (solid domain).While some artifacts are completely removed (white circles), some odd shapes are not corrected because of their large size (red circle).In addition, odd shapes may be corrected but still result in an over-partition of the particle (orange circle).In contrast, PCRF generates initially far less over-segmentation, with most of the artifacts located at particle interconnection, which results in a better final identification (also visible in Fig. 6).
Figure 13.Identification achieved for the solid (left) and pore (right) domains, with PCRF method and after over-segmentation correction, illustrated for a A12 graphite, SLC1506T graphite, and NMC microstructures.
generates a significant number of irrelevant artifacts with odd shapes, similar to those reported in the literature. 64In contrast, identification obtained with the PCRF is much less noisy.Little difference is observed between the identification obtained with different values of k (2, 3, and 4).There is then no real reason to choose one value over another.For the rest of the article, k = 3 has been chosen.
Over-segmentation correction and final identification.-Most of the small irrelevant artifacts are removed by using the oversegmentation correction algorithm discussed in paragraph section Discrete particle size distribution: pseudo-Coulomb repulsive field method, even for the watershed images.However, the final identification still appears visually better with the PCRF method rather than with the watershed method due to a better initial identification (cf.Fig. 12). Figure 13 shows an identification obtained with the PCRF method, after over-segmentation correction.No or very few artifacts are visible for the three different microstructures.
Until this point, all identification results have been calculated from two-dimensional slices.This allowed for visual evaluation and comparison of the two identification methods, as all of the domain topology is visible.However, quantifying particle morphology and pore topology requires work on three-dimensional volumes.From this point, all morphology and topology results are obtained from identification performed on three-dimensional volumes using the PCRF method with k = 3.To achieve 3D identification in a reasonable time, volumes have been cropped and image resolution has been downscaled (2×).Figure 14 shows such an identification.Visual interpretation of the identification results for the quasi-spherical NMC particles is possible in 3D.However, due to the complex particle 3D morphology of the two graphite electrodes, it is not possible to visually interpret identification results obtained from a 3D volume plotted in 2D slices.
Particle size comparison.-Particlesizes calculated with the numerical methods listed in Table II are shown in Fig. 15 for all the electrodes.NMC and SLC1520P electrodes have the largest particle sizes, while A12 and SLC1506T have the smallest and SLC1506T2, MCMB, and BFC-10 have intermediate sizes.The c-PSD systematically provides the lower bound due to three reasons.First, the spherical assumption induces a constraint that limits the particle size as discussed in Refs.48, 49, 73.For instance, the diameter of the largest sphere contained within an ellipsoid is the smallest ellipsoid's diameter.Second, the particle volume not contained within the largest sphere is assigned to smaller spheres, which induces an even smaller volume average particle size.Figures 16a and 16b    these two sources of underestimation.Particle elongation induces a significant particle size underestimation.For instance, for an ellipsoid of diameter ratio 1:1:2, the ellipsoid equivalent (respectively, largest sphere) diameter is ∼1.58 (respectively, ∼1.24) times higher compared with the mean diameter obtained with the c-PSD method.Because of these tiny sphere diameters, particle-size distributions calculated with c-PSD exhibit a peak for the very small diameters, absent for d-PSD algorithms such as watershed and PCRF (cf.Figs.17a and 17b).Particle elongation decreases sphericity as defined in Eq. 17, which controls the difference between the equivalent ellipsoid diameter and the mean c-PSD diameter (cf.Figs.16c and 16d).Flake-like particle (i.e., elongation 1:1:x with x < 1) induces a stronger decrease of sphericity and size underestimation compared to an oblong ellipsoid (i.e., elongation 1:1:x with x > 1).Third, inclusions or cracks of any size within a particle will drastically reduce the calculated mean c-PSD, as illustrated in Fig. 16b.Contrariwise, discrete particle-size algorithms such as PCRF are likely to consider that a relatively small crack or inclusion will not partition a particle into smaller particles (cf.Fig. 16b).Therefore, c-PSD strongly underestimates actual particle size, as already reported in the literature. 48,49,73Differences between c-PSD and d-PSD results are stronger for graphite than for NMC electrodes and will be discussed in the paragraph section c-PSD underestimation induced by the spherical assumption.
The largest diameters are calculated with the watershed and PCRF methods.Indeed, as opposed to c-PSD, these two d-PSD algorithms do not rely on any morphology assumptions.Similar diameters are obtained between the two methods for NMC, while higher values are reported for graphite electrodes with PCRF (cf.Fig. 15).This indicates that particle identification is roughly independent of the d-PSD algorithm for quasi-spherical particles such as NMC (cf.Fig. 13) after over-segmentation is corrected, while particle identification of more complex morphology is more dependent of the choice of the d-PSD algorithm, even after oversegmentation correction.Some particles are still over-partitioned after over-segmentation correction with the watershed method, indicating that their size is underestimated (cf.Fig. 12).This suggests that PCRF provides a better estimation of the particle diameters.Because of over-partitioning, watershed size distributions exhibit smaller highest diameters compared with PCRF (cf.Figs.17a  and 17b).Particle-size distribution indicates that both active material particles and pores have a wide unimodal diameter distribution (cf.Figs.17a and 17b), indicating that modeling phases with a unique mean diameter is an oversimplification.
The EDMF method calculated diameters similar with those obtained with the watershed method (cf.Fig. 15).Even though EDMF provides less information (no diameter spatial distribution and no particle morphology analysis) compared with the watershed method, it is also simple and much faster.Therefore, it is a very interesting method to characterize particle size of quasi-spherical particles (such as NMC), as EDMF, watershed, and PCRF provide similar diameters for such particle morphology.In addition, the distance-to-surface distributions, which is the definition of Euclidean distance map (cf.Eq. 11), are also plotted in Figs.17c and 17d.Such an alternative metric provides a relevant characteristic diffusion distance.
Covariance provides intermediate values among the different methods (cf.Fig. 15).However, to obtain a relevant fit of the equivalent diameter as detailed in Ref. 49, a large volume is required, which has prevented getting a diameter for some of the electrode volumes.The diameters reported in the Fig. 15 are the average of the diameters calculated along the three directions.
Lastly, the diameters have been also determined by first calculating the specific surface area of the active material particles (with the "direct" method as detailed in Ref. 49 with a corrective factor of 2/3), and then applying the relationship, diameter equal to six times the specific surface area inverse, which only stands for nonoverlapping spherical particles with no surface roughness.Even though the results are quite coherent with EDMF and watershed methods (cf.Fig. 15), and that such simplification is often used in the literature, 16 the authors strongly recommend not relying on this method as surface roughness can alter the result.For example, spherical particles in Fig. 10b have very similar diameters while having very different specific surface areas, due to a difference in surface roughness.To illustrate this, the graphite particle mean diameter of SLC1506T2 has been calculated at different voxel sizes, and thus different surface roughness, by applying a downscaling algorithm on the whole volume (cf.Fig. 18).As voxel size gets coarser, the surface details vanish while the particle volume is still preserved (as demonstrated by the nearly constant c-PSD diameters), and thus the specific surface area decreases and the diameter deduced from it increases.While diameters obtained from c-PSD are stable with voxel size (as particle volume is quite insensitive with image resolution), diameter obtained from the specific surface area is not (as surface is very sensitive to image resolution).Even worse, specific surface area is expected to diverge as imaging performed with a higher image resolution will reveal new, finer surface details not visible at a coarser image resolution.The analysis of the specific surface area fractal behavior is beyond the scope of this paper and questions the definition and dimension of specific surface area used in battery models for real microstructures.Indeed, should the battery community keep the current definition of this parameter, which results in a fractal behavior (thus a dimension between 2 and 3 for the surface, while using a dimension of 2 in the battery model equation), or redefine this parameter so that a finite converging value could be determined?To the authors' knowledge, this issue has not yet been Journal of The Electrochemical Society, 2020 167 100513 investigated thoroughly in the lithium-ion battery field, although some work has been done for fuel cell microstructures. 79Therefore, it is not recommended to rely on surface area to deduce particle size unless the particle surface is numerically smoothed prior to the calculation.
Particle morphology.-Activematerial particles and pores morphology metrics, corrected sphericity y c (cf. Eq. 17), and elongation are reported in Table IV.Both phases for most electrodes have very similar in-plane diameters, i.e., d d / with in-plane over through-plane diameters roughly identical and higher than one, i.e., ~>  d d  d d 1.
/ / These indicate that the active material and pore particles' typical morphology is an ellipsoid with a short diameter along the electrode thickness and identical in-plane diameters.However, pores and active materials particles do not share the same morphology, especially for the NMC, as their respective sphericities are different (discussed more in detail in the next paragraph).NMC active material particles are quasi-spherical with a relatively small elongation, while A12 graphite particles are significantly elongated, with their sphericity and elongation linearly correlated (cf.Fig. 19).One NMC data point (4-NMC-C) shows a significant offset, attributed to truncated particles that bias the particle morphology analysis.Indeed, 4-NMC-C is significantly thinner (34 μm) compared with the other NMC electrodes (88-160 μm) (cf.Table II of Ref. 16).Calendared NMC electrodes (noted NMC C) are slightly more elongated and less spherical than their uncalendared counterparts, (noted NMC UC) while having similar particle sizes.This indicates that calendaring slightly modified the particle morphology.SLC1520P, SLC1506T, SLC1506T2, and MCMB have intermediate sphericity and elongation compared with A12 and NMC electrodes.
Sphericity and elongation distribution.-Sphericityand elongation averaged values only provide a limited insight into the particle morphology.In-depth analysis has been performed on three electrodes, 1-NMC-C, SLC1506T2, and 5-A12-C, which respectively covers the high, intermediate, and low sphericity (cf.Fig. 19).Sphericity as a function of equivalent diameter is plotted for each individual active material particle and pore in Figs.20a and 20b.The region covered by the points {diameter, sphericity} is different for each microstructure, like a fingerprint.Indeed, these microstructures are significantly different: quasi-spherical for 1-NMC-C, channel-like for SLC1506T2, and cracked-potato-shape for 5-A12-C (cf.Fig. 20a inserts).Large active material particles tend to be less spherical for SLC1506T2 and 5-A12-C, while such a trend is not visible for 1-NMC-C.Point clouds and average metrics for 5-A12-C and SLC1506T2 are similar between active materials and pores (cf.Table IV and Figs.20a and 20b).Contrariwise, pores of 1-NMC-C are significantly less spherical.This result is derived from the active material particle morphology: the complementary volume of a channel-like phase (e.g., SLC1506T2) is also a channel-like structure, while the complementary volume of spheres (e.g., 1-NMC-C) is a star-shape, as depicted in the insert of Fig. 20b, which is much less spherical.Lastly, active material particles and pores of 5-A12-C have equally complex shapes.Differences in pore morphologies implies differences in permeability, i.e., the ability of the porous electrode to be filled with the electrolyte.Sphericity distributions are shown in Fig. 20c, with NMC and 5-A12-C having the narrower and the wider distributions, respectively, while SLC1506T2 exhibits an intermediate dispersion.
Similarly, particle elongation distributions have also been calculated for these three electrodes (cf.Fig. 21).Elongations d1/d2 are centered around one for both electrodes, elongations d1/d3 and d2/d3 are roughly identical for each microstructure (but different between electrodes), and through-plane over in-plane elongation are significantly different between electrodes (see Table IV).As with sphericity, distributions NMC and 5-A12-C exhibit the narrower and the wider distribution, respectively, while SLC1506T2 has an intermediate dispersion.This result indicates that active material and pore particle morphology heterogeneity is different between electrodes, with heterogeneities of 1-NMC-C, SLC1506T2, and 5-A12-C in ascending order.is plotted as a function of active material sphericity in Fig. 22.The difference increases as the sphericity decreases.Indeed, as sphericity is decreasing the spherical assumption of the c-PSD method is getting less relevant and the c-PSD underestimation is worsened, while the PCRF method is expected to provide a value closer to the actual diameter as it does not rely on any morphology assumption.This result clearly demonstrates that the c-PSD method does not provide a relevant lower bound for diameter but can be considered as an incorrect value when investigating nonspherical Note that a high ratio (2.8) has been also reported in the literature. 73The trend plotted in Fig. 22 is expected, as it measures the relevance of the spherical assumption and therefore provides an indirect validation of the particle identification step, as both sphericity and d PCRF 50 are derived from the particle identification.Discrete particlesize algorithms are not as available as the c-PSD method due to their higher complexities; therefore, one could use the ratios from Fig. 22 to provide a better estimation of the particle diameter based on a simple c-PSD calculation and a prior knowledge of the particle sphericity.For comparison, results obtained with ideal ellipsoids are also shown in Fig. 22.While the same trend is obtained, differences arise as sphericity does not uniquely characterize particle shape.and tortuosity factor t i along direction i are reported in Table V. Volume fraction from recipes for both pore and additives (carbon-black domain, noted CBD) are given for information.As stated in paragraph section Electrode library, over-segmentation, and particle identification, volumes are segmented so that the "pore" domain is the union of the void and the additive phases.If not specified otherwise, the term e refers to the volume fraction of the segmented combined phase pore and CBD.Constriction factor as geometrically defined in section Particle morphology and pore topology parameter determination (cf.Eq. 20) ranges from 0 (worst) to 1 (ideal).Geometric tortuosity as geometrically defined in the same paragraph section (cf.Eqs.21 and 22) ranges from 1 (ideal) to plus infinity (worst).Tortuosity factor as indirectly defined in the introduction section (cf.Eq. 1) ranges from 1 (ideal) to plus infinity (worst) and has been calculated using the open-source MATLAB software package TauFactor, 24 as done in a previous work. 16If the additives were to be considered, higher tortuosity factors would be obtained as detailed in previous work. 16Among the investigated electrodes, porosity ranges from 0.387 to 0.617.For most electrodes, in-planes metrics are similar (i.  .Correlation between these parameters is specifically investigated in paragraph section Pore tortuosity factor correlation with microstructure geometry.Geometric tortuosity distribution and asymmetry.-Contraryto tortuosity factor, for which a unique value is obtained per volume and direction, geometric tortuosity offers more granularity as its value is the average of all the shortest paths as illustrated in Fig. 1a.Furthermore, geometric tortuosity may be asymmetric, as discussed in paragraph section Review of microstructure-tortuosity relationship and section Particle morphology and pore topology parameter determination (cf.Eq.21 and Fig. 1a).Geometric distributions are plotted for three  / and (d) electrodes in Fig. 23.Very little asymmetry is obtained, as expected for such large volumes.Narrow distributions are calculated for all geometric tortuosity except for the through-plane direction of 5-A12-C.This result indicates that 5-A12-C graphite microstructure exhibits a geometric tortuosity heterogeneity, i.e., depending on their entry point, lithium ions have a significantly different shortest path length to reach the current collector from the separator, and vice versa.In other words, submillimeter in-plane regions of the electrode have access to significantly different shortest path to reach the opposite in-plane region.It can be assumed that diffusion path length will also be heterogeneous for this electrode.
Pore tortuosity factor correlation with microstructure geometry.-Previoussections have been focused on analyzing the particle morphology and pore topology metrics independently.In this section, they will be correlated with the tortuosity factor to deconvolute each microstructure geometric parameter that contributes to its value.
Tortuosity-porosity correlation.-Tortuosityfactors are plotted as a function of porosity in Fig. 24 for all electrodes and directions.The generalized Archie's relationship (cf.Eq. 3) is fitted for NMC and A12 electrodes separately.As expected, the correlation stands only when limited to a single type of microstructure.However, it clearly appears that porosity alone cannot explain the variation of tortuosity factors when all different microstructures are considered, especially along the electrode thickness.Furthermore, while both inplane tortuosity factors are rather similar, through-plane tortuosity factors are much higher.Limitations of tortuosity-porosity correlation are detailed in the introduction section.
Tortuosity-microstructure correlation.-Correlationwith porosity, geometric tortuosity, constriction factor, and/or geometric constrictivity.Tortuosity factors are plotted as a function of geometric tortuosity and constriction factor along the through-plane direction ( i.e., along electrode thickness) in Figs.25a and 25b.As expected, as geometric tortuosity (respectively, constriction factor) is increasing (respectively, decreasing) tortuosity factor is increasing.Linear and inverse correlations, respectively, are established between t 3 and t geo,3 (cf.Eq. 23) and t 3 and b 3 (cf.Eq. 24).Geometric tortuosity provides an underestimation of tortuosity factor since the two other geometric contributions to the tortuosity factor (porosity and constrictivity) are not yet considered together.While some points show significant offset from the fitted curves, as expected since only one geometric contribution is considered, it is noticeable that geometric tortuosity or constriction factor alone is enough to fit the tortuosity factor with a monotonic function.This result suggests that all three parameters (porosity, constriction factor, and geometric tortuosity) are not independent.Parameter orthogonality is investigated later in this section.Through-plane over in-plane tortuosity factor anisotropy can also be fitted with geometric tortuosity and constriction factor   , (cf.Eq. 8) and plotted as a function of the constriction factor b i (cf.Fig. 26).The aim is to find a relationship between d i (a nongeometric parameter, as t i ) and b i (a geometric parameter) so that a new parameter, the geometric constrictivity d , geo i , is defined as a function of b i (cf.Eq. 9).Therefore, the constriction contribution to the tortuosity factor is estimated with , (a geometric parameter as defined as a function of b i ) and no longer with d .
i Once this empirical relationship is established (cf.Eq. 26), the normalized effective diffusion coefficient can be predicted knowing the geometric tortuosity, geometric constrictivity (through the constriction factor), and the porosity (cf.Eq. 9), i.e., only geometrically defined parameters.The correlation d b = f geo i i , ( ) is performed considering all the electrodes and all the directions together (cf.Fig. 26).
´-´+ 9.237 0.924 0.593 26 [ ] Figure 27 shows different tortuosity-microstructure correlations for all the electrodes and all directions.If only the porosity is considered, a generic correlation does not exist (cf.Fig. 27a).We define generic correlation as a relationship that stands for a wide range of microstructures and for all orientations.If the porosity, geometric tortuosity, and constriction factor are considered, a generic relationship exists (cf.Fig. 27b and Eq.27).To achieve a one-to-one prediction, the porosity, geometric tortuosity, and Correlation with particle elongation.Through-plane tortuosity factors t 3 are plotted as a function of the active material particle elongation 1:1:x (x being the through-plane diameter) in the insert of Fig. 28.While a decent correlation is achieved, which suggests that particle elongation is correlated with geometric tortuosity and/or geometric constrictivity (as discussed later in this section), particle elongation alone is insufficient to accurately predict the throughplane tortuosity factor.Indeed, the porosity difference between the calendared and uncalendared electrode induces a variation (cf.Fig. 28a).To improve the prediction, tortuosity factor is fitted with the product of the particle elongation with the square root of the porosity (simply multiplying by the porosity provided a less correct fit) (cf.Eq. 28 and Fig. 28a).Differences from the fitted line arise from the truncated particles (edges effect) that bias the particle elongation calculation.For instance, the NMC-UC offset point (of elongation ∼0.7) corresponds to the 4-NMC-C electrode, which is significantly thinner (34 μm) compared with the other NMC electrodes (88-160 μm).SLC1520P has the largest particle size (cf.Fig. 15), which induces a stronger edge effect compared with the other electrodes with a similar field of view and explains the deviation from the fitted line.In addition, particle elongation and porosity predict tortuosity factor anisotropy well (cf.Fig. 28b).Parameter orthogonality.Ideally, to deconvolute the impact of each microstructure contribution to the effective diffusion coefficient, independent parameters should be used.Results from Figs. 25  and 28 suggest that geometric tortuosity, geometric constrictivity, and particle elongation are correlated.Particle elongation and porosity are not correlated in the design space investigated.Geometric tortuosity and porosity are not correlated when all directions are considered.Geometric constrictivity generally increases with porosity, but dispersion is too high to establish a proper relationship.Contrariwise, geometric tortuosity is correlated with particle elongation (cf.Fig. 29a) and geometric constrictivity is correlated with geometric tortuosity (cf.Fig. 29b).Since porosity is not correlated with the three other metrics, any tortuosity-microstructure relationship must consider porosity.Geometric tortuosity can be deduced from particle elongation (cf.Eq. 30), and geometric  constrictivity can be deduced from geometric tortuosity (cf.Eq. 31).Therefore, knowing one of these three parameters is enough to determine the other two.As a consequence, to predict tortuosity factor and thus the normalized effective diffusion coefficient, either {porosity, particle elongation}, {porosity, geometric tortuosity}, or {porosity, geometric constrictivity} is enough, as illustrated in Fig. 28a (Eqs.28 and 29) and 30 (Eq. 32).Journal of The Electrochemical Society, 2020 167 100513 on two-dimensional images, as it allows visual inspection.Furthermore, the original method (PCRF) used in this work to provide all the particle morphology and pore topology metrics has been compared with a reference method (watershed) and proved to be better, based on over-segmentation for both test case and CT electrode geometries.The physics-based PCRF identification method is more robust, as the field F ¯it relies on is deduced locally from multiple distances, while the watershed method relies on a field EDM deduced locally from a unique distance (the shortest distance) and is then more sensitive to surface roughness, as thoroughly explained in paragraph section Discrete particle size distribution: pseudo-Coulomb repulsive field method.The absence of morphology constraint on d-PSD algorithms, such as watershed and PCRF, is both an advantage, as it prevents biased results due to incorrect assumptions, and an hindrance, as it is difficult to visually judge the identification quality for complex (i.e., non-spherical particles) 3D microstructures such as the various graphite electrodes investigated here.Since there is not a unique definition of a particle within a fully connected cluster, particle identification is by definition very difficult to validate, as the ground truth for real microstructures is not defined.Nevertheless, the various relevant correlations established from metrics obtained from particle identification indirectly support the adequacy of the identification results.Namely, the expected variations between particle elongation and sphericity (cf.Fig. 19); c-PSD underestimation and sphericity (cf.Fig. 22); tortuosity factor, geometric tortuosity, and constriction factor, including their anisotropy (cf.Fig. 25); and tortuosity factor and particle elongation, including anisotropy (cf.Fig. 28).Furthermore, the few offset points are explained by edged effects (truncated particles), (cf.Fig. 28a).Lastly, the fact that simple (i.e., linear or quadratic) monotonic (i.e., strictly increasing or decreasing) correlations of porosity, constriction factor (or geometric constrictivity), and geometric tortuosity were enough to match the normalized effective diffusion coefficient (cf.Figs.27b and 27c) supports the identification results.Indeed, while it is possible to fit a given function with any non-defined parameter, the resulting correlation may be extremely sensitive, indicating overfitting (i.e., meaningless correlation) has been done. 80Thus, the combination of simple monotonic relationships with meaningful (i.e., defined) variables is a strong indicator that no or few overfittings have been performed in this work.More complex relationships have also been established (cf.Fig. 30), but only when these metrics were correlated between one other, i.e., reducing the number of parameters using their nonorthogonality but at the price of an increase of the correlation complexity.
The identification algorithm could be improved by distinguishing stable and unstable end displacement points (cf.Fig. 4b) by applying small perturbations.Particles would be then defined as the sum of trajectories that converge towards stable points, while unstable end points would be assigned to the nearest identified particles.Such a modification would improve voxel assignments near particle connections.
The identification step could be greatly simplified upstream by coating particles with a thin surface layer distinctly visible with CT prior to the electrode fabrication and specifically for the imaging, as particle boundary would be clearly visible.Furthermore, such an experiment could be used to validate the identification algorithm, as it essentially provides the ground truth.
Tortuosity-microstructure relationships relevance and possible improvements.-Beingable to establish simple monotonic relationships supports the fact that porosity, sinuosity (through the geometric tortuosity), and constriction (through the geometric constrictivity and the constriction factor) are relevant parameters to geometrically define and predict the normalized effective diffusion coefficient for a wide variety of microstructures (cf.Figs 27b, 27c,  and 30a).However, they are not the only geometrically defined metrics that can predict it.For instance, particle elongation in combination with porosity can be used instead (cf.Figs. 28, and  30b).This indicates that there is not a unique tortuosity-microstructure generic relationship, but several of them.This is also supported by the large variety of tortuosity-microstructure relationships reported in the literature (cf.introduction section).Therefore, it is likely that geometrically defined metrics other than those considered in this work could be used to correlate microstructure topology with transport properties.Our work suggests that a minimum of two independent (i.e., orthogonal) parameters, one being the porosity, is enough to characterize the normalized effective diffusion coefficient and the tortuosity factor.However, additional parameters are required to account for the role of lower-scale solid phase, if any, as done in previous work 16 for the carbon-black additives.It is still an open question if sinuosity and constriction  / All directions are considered (i.e., i = 1, 2, and 3).
could be defined so that their values are not correlated once calculated (i.e., orthogonal).All microstructures investigated in this work roughly share the same attribute: active material particles can be characterized (on average) as ellipsoids with their two inplane diameters similar and their through-plane diameter slightly or significantly smaller.While methods and parameters defined in this work are expected to correctly characterize microstructures for which particles have three distinct diameters, as both geometric tortuosity and constriction factor are oriented parameters, the relationships established in this work may be adjusted for such different microstructures.
The geometric constrictivity approach relies on an assumption which is not guaranteed to be true: the existence of a monotonic relationship between the constriction factor and the constrictivity (cf.Fig. 26 and Eq.26).A monotonic relationship is a desired feature, as it guarantees a relative stability for the prediction (i. ), and therefore a lower requirement on the constriction factor precision.The relative dispersion obtained (cf.Fig. 26) indicates that either the definition, the method to calculate it, both the definition and the method, or the field of view should be refined, improved, or enlarged to reduce this scattering.However, the geometric constrictivity approach is not mandatory to provide accurate prediction (cf.Figs.27b and 28).
Geometric tortuosity alone provides a decent correlation with the tortuosity factor (cf. Figs.25a and 25c).However, it must be strongly corrected (scaled) in order to provide a better estimation (cf.Eq. 23).A possible approach to lower this underestimation would be to define a graph weighted not only by the length between nodes, but also by the local section area or local section area variation.Note that this modification would essentially consist in merging both geometric tortuosity and constriction in one parameter.As these two parameters are correlated (cf.Fig. 29b), this parameter reduction is justified.
Tortuosity-microstructure relationships practicability.-Tortuosity-microstructurerelationships implying geometric tortuosity and constriction factor or geometric tortuosity (cf.Eqs. 23, and 27) are essentially valuable from an understanding point of view, i.e., deconvoluting the different microstructure contributions on the normalized effective diffusion coefficient or the tortuosity factor.However, they have little value from a design guidance point of view, as recommending decreasing the geometric tortuosity or increasing the constriction factor is as pointless as recommending decreasing the tortuosity factor for improving the effective ionic diffusion for the fast charging of thick electrodes: without additional information, it is unclear how to control these parameters, and their ideal values are already known.
From a calculation point of view, relationships involving porosity and geometric tortuosity (cf.Fig. 30a and Eq.31) can be interesting.Indeed, while calculating geometric tortuosity through pore identification is a complex task, another path consists simply in calculating the shortest paths considering all voxels, at the price of a much larger graph and thus larger CPU and memory requirements.Lengths of identified shortest paths should also be corrected to remove the staircase bias induced by voxels.Such an approach (convert pore domain in a graph, calculate shortest paths, determine geometric tortuosity, and then deduce normalized effective diffusion coefficient) could replace the Laplace homogenization approach (convert pore domain in a mesh, solve Laplace equation, and then deduce normalized effective diffusion coefficient).
From a design guidance point of view, the most interesting relationships are those involving porosity and particle elongation, as both are controllable to a certain extent.One can use Eqs.28 and 29, which correlates directly with tortuosity factor (as well as tortuosity factor anisotropy), with particle elongation and porosity, or Eq.32, which takes advantage of the non-orthogonality between particle elongation, geometric constrictivity, and geometric tortuosity.In both cases, it provides a much better prediction of the tortuosity factor and of the normalized effective diffusion coefficient compared with relying only on porosity.It also provides a strong incentive on controlling the particle alignment during the manufacturing process and a simple explanation of why electrodes with the same porosity exhibit significant different effective ionic transport properties.In this work, particle elongation required the complex particle identification step to be determined.
Particle size relevance.-Particlediameters have been thoroughly compared in paragraph section Particle size comparison and c-PSD underestimation induced by the spherical assumption.
Because "particle" is not uniquely defined within a fully connected cluster, a wide range of numerical methods exist to calculate its diameter.However, these methods are not equivalent (cf.Fig. 15) and especially two should be avoided: (i) c-PSD is strongly biased due to an incorrect assumption and should be disregarded for nonspherical particles (cf.Fig. 22), and (ii) diameter deduced from specific surface area suffers from a fractal behavior (cf.Fig. 18) that disqualified it.The EDMF original method presented in paragraph section Euclidean distance map fitting method provides a valuable alternative because it is easy to implement, fast, and provides a similar mean diameter with the more complex watershed method.Discrete particle size algorithms, such as the new PCRF method, are mostly valuable when the particle morphology is far from spherical (e.g., graphite A12) or for specifically investigating the microstructure topology as done in this work.

Conclusions
Particle identification has been performed on a set of 19 various LIB electrode (NMC and graphite) three-dimensional microstructures obtained from CT images with a wide range of porosity, using an novel algorithm, the PCRF method, in order to quantify particle size and shape as well as pore topology to ultimately correlate tortuosity factor and normalized effective electrolyte diffusion coefficient with geometrically defined metrics.Established correlations improve understanding of the different contributions of the microstructure that control the effective ionic diffusion and provide ways to improve it, which is valuable for fast charging of thick LIB electrodes.Main elements are recapitulated below.
Particle identification.-Thenew particle identification algorithm (PCRF) introduced in this work consists of dropping charged points between a fixed domain's boundary also charged with the same polarity and identifying particles based on the trajectory of these charged points according to a pseudo-Coulomb's law.This new method is much less sensitive to surface roughness and geometry symmetry in regards with over-segmentation compared with the watershed reference method for both test-case geometries (cf.Fig. 10 and Table II) and CT-imaged battery electrodes with complex microstructures (cf.Figs.11 and 12).We demonstrate that surface roughness is more of an issue for particle identification in the case of channel-type geometries than for particle-type geometries (cf.Fig. 10).While over-segmentation is typically attributed to surface roughness, we illustrate that symmetry of simple test geometries is another source of over-segmentation, as breaking such symmetry strongly reduces the over-segmentation (cf.Table II).Testing identification algorithms on only symmetric geometries may bias analysis, as symmetry is rarely an attribute for LIB electrode microstructures.
Particle size and shape.-Inaddition to the new PCRF method, a second original numerical algorithm has been also introduced in this work, although limited to calculating particle size: the EDMF method.EDMF, numerically simple and fast, provides a lower diameter bound with values similar to those obtained with the more CPU-expensive watershed method, while PCRF, numerically complex and slow, provides a higher diameter bound, with values expected to be closer with the actual diameters as watershed exhibits more particle over-partitioning in comparison (cf.Fig. 12).
Journal of The Electrochemical Society, 2020 167 100513 Diameters have been calculated with six different methods: specific surface area, c-PSD, covariance, EDMF, watershed, and PCRF (cf.Fig. 15).Particle size obtained through specific surface area proved to be unreliable due to fractal behavior of the surface area (cf.Fig. 18).Particle size underestimation of the c-PSD method has been correlated with particle sphericity, as it quantifies the relevance of the spherical assumption (cf.Fig. 22).The covariance method offers intermediate values.EDMF and watershed provide similar larger diameters.Lastly, PCRF gives the largest diameters, expected to be the closest to the actual diameters among the other methods.
Particle morphology has been deduced from the particle identification step through two metrics: sphericity and elongation or aspect ratio.On average, the particle aspect ratio is anisotropic with similar in-plane dimensions and a slightly (e.g., NMC) or significantly (e.g., A12 graphite) smaller through-plane dimension (i.e., along the electrode thickness) depending on the electrodes.Particle sphericity and diameter provide a signature of the medium investigated, as each different type of microstructure exhibits a different distribution (cf.Fig. 20a).Besides, active materials and pores may or may not share the same distribution (cf.Fig. 20b).Results show significant particle size and shape differences between the different microstructures; for instance, those in NMC are larger and more spherical compared with A12 graphite.Furthermore, morphology heterogeneity is different across all the electrodes microstructures, with more heterogeneity calculated for the A12 graphite (cf.Figs.20c and 21).
Tortuosity-microstructure relationships.-Aset of oriented and geometrically defined metrics has been determined from the particle morphology analysis (particle elongation) and from the pore topology analysis, applying particle identification algorithm on the pore domain (constriction factor, geometric constrictivity, and geometric tortuosity).Generic tortuosity-microstructure relationships that stand for all microstructures and all directions (unlike the generalized Archie's relationship that is microstructure specific; cf.Figs.24 and 27a) were then established: normalized effective diffusion coefficient with porosity, constriction factor or geometric constrictivity, and geometric tortuosity (cf.Figs.27b and 27c); porosity and geometric tortuosity or particle elongation (cf.Fig. 30); and tortuosity factor with porosity and particle elongation (cf.Fig. 28).Geometric tortuosity, geometric constrictivity, and particle elongation have been found to be correlated (cf.Fig. 29).Providing such a generic relationship was the initial aim of this work, and their relevance and interest, e.g., in terms of design guidance, is debated in the discussion section.
In summary.-Anovel method to identify particles or pores embedded in a connected cluster was developed and proved to be more robust than the reference watershed method.Particle morphology and pore topology analysis enabled by the aforementioned new identification algorithm was carried out to establish generic correlations between geometrically defined metrics and the tortuosity factor and normalized effective diffusion coefficient on various LIB electrode three-dimensional microstructures from an open-source library.It is expected that the new algorithm will be useful for future microstructure characterization of heterogenous materials, and that the tortuosity-microstructure relationships will provide relevant design guidance to improve ionic transport of LIB electrodes.

Figure 1 .
Figure 1.(a) Tortuosity factor is a combination of both the geometric tortuosity t geo and the constriction factor b. In this example, there are three shortest paths from bottom to top noted + - l i ,1 3 and four from top to bottom noted -- l ; i ,1 4 (b) Tortuosity factor due to the constriction factor; (c) Tortuosity factor due to the geometric tortuosity; (d) Unit tortuosity factor.

1 ,while b  1 Holzer(
Note that Holzer et al. have chosen the opposite area ratio for the constriction factor compared with the initial Petersen's definition: b  Petersen cf.Eqs. 5 and 7).The constrictivity d introduced by Brakel and Heerjes

Fig. 1
). Parameter orthogonality for actual microstructures is discussed in detail in the paragraph section Tortuosity-microstructure correlation.

Figure 2 .
Figure 2. Tortuosity illustrated for periodically aligned ellipsoidal particles for diffusion path (a) normal with the ellipsoid long diameter and (b) aligned with it.
Particle elongation or aspect ratio: d 1 and d 2 in-plane dimensions, d 3 through-plane dimension >0

Figure 4 .
Figure 4. (a) Illustration of the calculation of the pseudo-Coulomb field F ¯for a case-study geometry.The domain's wall is positively charged (+).(b) Actual stepwise displacements of charged points (in magenta) as used in the algorithm according to the sign of F. ¯Stable and unstable positions (end of displacement) coincide with the local minimum of F .  (c) Displacements of charged points (magenta, orange, and grey) according to the gradient of F ¯(not used in the algorithm).
Particle elongation is illustrated in Fig.7.

Figure 5 ./
Figure 5. Voxel assignment to a particle numerated p illustrated during charged-point displacement.(a) Standard displacement, (b) displacement reaches F   local minimum, and (c) current trajectory (in orange) intersects with a previous trajectory (in magenta).Trajectories are merged: the particle is the sum of all the trajectories which intersect themselves.

Figure 6 .Figure 7 .
Figure 6.Example of post-processing steps illustrated on the pore domain of an NMC electrode CT image.(a) From left to right: (i) binary image, (ii) pore segment identification (PCRF method) with chessboard patterns (circled in red) localized at the bottlenecks between some neighboring segments, and (iii) pore segment identification after chessboard pattern removal.Please note chess pattern is not systematic.(b) From left to right: (i) binary image; (ii) pore segment identification (PCRF) with small irrelevant segments identified at the bottlenecks.Some voxels of these segments could be assigned to larger pore segments circled in red; (iii) Voxels for which - - d PSD x c PSD x ( ) ( ) are marked; and (iv) reassigned to neighbored pore segments iteratively until no voxels verify -- d PSD x c PSD x .( ) ( )

,
) (cf.Eq. 21).As already discussed in paragraph 1.1 and illustrated in Fig.1, t + geo i , and t - geo i , may not be equal.Their difference is discussed in paragraph section Pore topology.In this work, the geometric tortuosity along the direction is the mean value of t +

Figure 8 .
Figure 8.(a) Pore segment identification (PCRF method) illustrated for a pore domain two-dimensional slice of an NMC electrode (white background is the solid particles).(b) Undirected graph of the pore domain weighted with the edge distance.Pore segment equivalent diameter is represented by a disc.(c) Graph visualization of a pore three-dimensional domain.

2 .
field F ¯(illustrated with k = 2) Particle identification c,d) F   a,b) a) Local minima (dark blue) indicate the centers of attraction.b) Iso-lines are shown.Iso-values are following a geometric progression for visibility's sake.c) Results are shown without the over-segmentation correction discussed in paragraph section Discrete particle size distribution: pseudo-Coulomb repulsive field method and illustrated in Fig. 6.If over-segmentation occurs, the particle identification obtained after its correction is shown in Fig. 9. d) PCRF is used with parameter = k If particle identification obtained with = k 3 differs, then it is also displayed.

Figure 9 .
Figure 9.Comparison between the particle identification obtained with watershed method using (a) our in-house algorithm and (b) the MATLAB built-in watershed function.Over-segmentation correction allows for removal of most of the over-segmentation for the chosen example (residual tiny artifacts are circled in red).The number of identified particles is noted n.

Figure 10 .
Figure 10.(a) Over-segmentation induced by surface roughness applied to perfectly smooth geometries.(b)-(c) No over-segmentation is induced by surface roughness applied to rounded geometry, except for some tiny features located at the particle surface which are insignificant for volume averaged analysis.The number of identified particles is noted n.

Figure 11 .
Figure 11.Particle identification without over-segmentation correction between watershed (reference) and PCRF methods for (top left) A12 graphite, (top right) SLC1506T and (bottom) NMC microstructures, applied for solid (top rows) and pore (bottom rows) domains.Red circles indicate some over-segmentation and odd particle shapes. illustrate

Figure 14 .
Figure 14.Identification obtained on three-dimensional volumes (black is pore domain).

Figure 15 .
Figure 15.(Top) Active material mean particle sizes calculated on various NMC and graphite electrodes, using c-PSD, specific surface area S p ( = d S 6 p 50 / ), EDMF, watershed, PCRF, and covariance; (bottom) Active material mean particle size normalized with c-PSD.

Figure 16 .
Figure 16.(a) c-PSD size underestimation illustrated with an ellipsoidal particle; (b) c-PSD size underestimation and corrected sphericity calculated as a function of ellipsoid elongation 1:1:x for one ellipsoid; (d) c-PSD size underestimation as a function of corrected sphericity; (e) Diameter map comparison between c-PSD and PCRF methods for a SLC15020P particle with a tiny void inside.

Figure 17 .
Figure 17.SLC1506T2 size distribution of (a) active material particles and (b) pores calculated with c-PSD, watershed, and PCRF methods; (c) Active material particle and (d) pore distance to boundary distribution of SLC1506T2 (not to be confused with diameter).Distributions have been smoothed using a moving average filter.(c) Insert shows the Euclidean distance map cumulative function (C , EDM cf.Eq. 12) and the fitted cumulative function for a sphere (C , EDM sphere

Figure 18 .
Figure 18.(a) Active material mean particle diameter and (b) specific surface area calculated as a function of voxel size for SLC1506T2.Volume fractions are conserved within the range of voxel size investigated: porosity varies from 0.400 to 0.407.

Figure 20 .
Figure 20.Corrected sphericity as a function of equivalent diameter for (a) active material particles and (b) pores.(a) Inserts show a 2D slice extracted from a subvolume of the electrode segmented volume.(c) Distribution of corrected sphericity (sphericity is weighted with the particle volume).

Figure 22 .
Figure 22.Active material particle mean particle size calculated with PCRF, normalized with diameter obtained with c-PSD, plotted as a function of the active material particle mean corrected sphericity.Dashed line is a linear correlation, while solid blue lines are copied from Fig.16dto provide comparison between results from ideal ellipsoid and actual microstructures, with x the ellipsoid elongation 1:1:x.

Figure 23 .
Figure 23.Geometric tortuosity distribution for (a) 1-NMC-C, (b) SLC1506T2, and (c) 5-A12-C.Solid lines are for t - geo i , (i.e., paths calculated from top to bottom) and dashed lines are for t + geo i , (i.e., paths calculated from bottom to top).

Figure 24 .
Figure 24.(a) First in-plane, (b) second in-plane and (c) through-plane tortuosity factors as a function of porosity for all electrodes.

Figure 25 .
Figure 25.Through-plane tortuosity factor t 3 as a function of (a) through-plane geometric tortuosity t geo,3 and (b) through-plane constriction factor b ;3

Figure 26
Figure 26.Constrictivity d t t = i g e oi i , / as a function of the constriction factor b .iDashed line is the geometric constrictivity d geo i , fitted so that it matches the constrictivity, i.e., d d .geoi i, ˜All directions are considered (i.e., i = 1, 2, and 3).

Figure 27 ./
Figure 27.(a) Tortuosity factors t i as a function of porosity e; (b) normalized effective diffusion coefficient e t i / as a function of e b t ´; i geoi , / and (c) normalized effective diffusion coefficient e t i / as a function of e d t ´. geo i geo i , , / All directions are considered (i.e., i = 1, 2, and 3).

Discussion
Particle identification relevance and possible improvements.-Asignificant effort has been done to assess the identification results

Figure 28 .
Figure 28.(a) Through-plane tortuosity factor t 3 and (b) tortuosity factor through-plane in-plane anisotropy as a function of the product of the active material particle elongation 1:1:x, with = + d d d x 2 3 1 2 (( ) ) / / with the square root of the porosity e.(a) (Insert) Less accurate correlation is achieved when only the particle elongation is considered.

Figure 29 .
Figure 29.(a) Geometric tortuosity t geo i , plotted as a function of particle elongation x i and (b) geometric constrictivity d geo i , plotted as a function of geometric tortuosity.All directions are considered (i.e., i = 1, 2, and 3).

Figure 30 .
Figure 30.Normalized effective diffusion coefficient e t i / as a function of (a) porosity e and geometric tortuosity t geo i , and (b) porosity e and particle elongation x i with = + Journal of The Electrochemical Society, 2020 167 100513 TableII.Summary of the different particle size methods employed in this work with pros (+) and cons (−).Not only is average diameter d 50 known, but also the local diameter d x .
is reused.It corresponds to a variety of LIB 7 cathodes Li(Ni 0.5 Mn 0.3 Co 0.2 )O 2 from TODA America Inc. (noted NMC in this work) and 7 anodes CGP-A12 graphite from ConocoPhillips Inc. (noted A12 graphite in this work).Calendered (noted C) and uncalendered (noted UC) electrodes were fabricated by the Cell Analysis, Modeling, and Prototyping (CAMP) facility at Argonne National Laboratory and imaged by University College London with micro-and nano-CT.Details of the preparation, imaging, and segmentation can be found in Ref. 16. Segmented volumes investigated in this work correspond to the case where the background domain volume fraction is the combined

Table IV .
Particle morphology characterized with corrected sphericity and elongation for active materials and pores for all electrodes.Mean values have been weighted with the particle or pore volume.

Table V .
Pore volume fraction and tortuosity-related metrics.Journal of The Electrochemical Society, 2020 167 100513 anisotropy (cf.Fig.25candEq. 25).This indicates that the oriented definition of the geometric tortuosity and constriction factor matches well with the diffusion anisotropy of the porous electrodes.Constrictivity d i is deduced from the tortuosity factor t i and the geometric tortuosity t geo i