Tortuosity in electrochemical devices: a review of calculation approaches

ABSTRACT The tortuosity of a structure plays a vital role in the transport of mass and charge in electrochemical devices. Concentration polarisation losses at high current densities are caused by mass transport limitations and are thus a function of microstructural characteristics. As tortuosity is notoriously difficult to ascertain, a wide and diverse range of methods have been developed to extract the tortuosity of a structure. These methods differ significantly in terms of calculation approach and data preparation techniques. Here, a review of tortuosity calculation procedures applied in the field of electrochemical devices is presented to better understand the resulting values presented in the literature. Visible differences between calculation methods are observed, especially when using porosity–tortuosity relationships and when comparing geometric and flux-based tortuosity calculation approaches.

Effective diffusion coefficient of a species D bulk Bulk diffusion coefficient Binary diffusion coefficient e i aPDF Speed of a species in the particle distribution function Effective heat flux R Ideal gas constant r 2 (t) Mean square displacement r P Mean pore radius T Temperature t Time V phase Volume fraction of analysed phase w Mass fraction x Mole fraction x PDF Location of a species in the particle distribution function Symbols α Bruggeman exponent α PDF Direction of movement of a species in the particle distribution function γ Scaling factor δ Constrictivity ε Porosity 1/t 2 Diffusibility or effective relative diffusivity κ Tortuosity factor κ geo Geometric-based characteristic tortuosity factors κ CFD Flux-based characteristic tortuosity factors λ bulk Bulk thermal conductivity μ Dynamic viscosity σ bulk Bulk conductivity σ eff Effective conductivity τ Tortuosity τ C Characteristic tortuosity V i aPDF Collision term of a species in the particle distribution function Introduction Electrochemical devices, including fuel cells and batteries, will play an increasing role in our lives, particularly as we transition to a low-carbon economy. However, in order to accelerate their commercialisation across a range of applications, an improved understanding of the underlying material characteristics is required. The importance of the effect of microstructure on the performance of electrochemical devices has been widely demonstrated [1], which is why studies of microstructural analysis techniques [2] are crucial for optimising vital parameters. Among these parameters, tortuosity plays an essential role in mass transport and concentration polarisation resistance [3,4]. Yet, calculating the tortuosity is not trivial, which is why a wealth of tortuosity calculation methods have been developed, not only in the electrochemical community, but also across many fields of research (optics, magnetism, geology, medicine, etc.), each with associated definitions and areas of application [5,6].
The microstructure of porous electrode and support layers in electrochemical devices is a major contributor to performance losses, especially in mass transport limiting operating regimes. This is valid for batteries, fuel cells, super-capacitors and oxygen transport membranes alike, where tortuosity is used to relate the effective transport properties of diffusion and electric or ionic flux, to its respective bulk property. As such, tortuosity is an indispensable parameter in modelling and quantifying fuel cell [7] and battery behaviour [8]. In addition, tortuosity serves as an input parameter in Newman-type models of battery performance [9] and the Adler-Lane-Steele model for electrode kinetics [10].
Owing to the importance of tortuosity for electrochemical devices, and the multitude of calculation approaches, this article reviews the use of tortuosity and various tortuosity calculation methods in the field of electrochemistry. These methods can differ considerably from each other in terms of calculation approach and data preparation techniques. Here, we deal with each in turn.

Definition of tortuosity
In geometrical terms, tortuosity τ is defined as the fraction of the shortest pathway through a porous structure Δl and the Euclidean distance between the starting and end point of that pathway Δx, illustrated in Figure 1 and equation (1). As such, τ always amounts to a value equal to or greater than unity.
In general, when analysing a porous structure, there exists only one shortest pathway and one tortuosity value. From this geometric perspective, constrictions or bottlenecks of the pore structure are not considered. The concept of tortuosity has been adopted in a variety of sciences [5,6], such as gaseous mass transport and electronic and ionic conductivity through porous, functional layers. In these fields, tortuosity is applied in a broader way than just a simple geometric measure of the shortest path length; tortuosity is also used to quantify and describe the resistance of a structure to a flux. In this respect, the difference between 'tortuosity' and 'tortuosity factor' was coined by Epstein in 1989 [11], who used a capillary model to show that the tortuosity τ is the square root of the tortuosity factor κ, as presented in Eq. (2).
The tortuosity factor accounts for both the additional path length and its change in the velocity of a species when migrating through a porous structure. Epstein then applied this derivation in the field of diffusion, where the tortuosity factor is used to calculate the effective diffusion coefficient D eff based on the bulk diffusion coefficient D bulk , shown in Eq. (3), which is also valid for ionic and electronic conductivities. The nomenclature distinguishing between κ and τ is adopted in this review, and values are converted accordingly, where necessary.
Yet, the theory behind the tortuosity factor is controversial, especially in the field of diffusive mass transport: van Brakel and Heertjes [12], for example, defined a constrictivity factor δ to account for the variation in pore diameter along the diffusion pathway, which is included in calculating the effective transport property via Eq. (4).
This constrictivity factor was later adopted by Holzer et al. [13] who stated that the implementation of τ 2 was used to explain high values of experimentally derived tortuosities. Consequently, the authors differentiated between two types of tortuosity: [13,14] (1) That, which is acquired by indirect calculations based on experimental data τ exp . (2) And that, which is determined via geometric algorithms from reconstructed 3D volumes τ geo .
Additionally, when analysing diffusive mass transport problems, depending on the diffusion mechanism taking place through a porous medium (ordinary, Knudsen and/or viscous diffusion) [15] and on the gases involved [16], different tortuosity values may dominate.
The geometric definition of tortuosity clearly suggests that there exists only one shortest pathway through a porous membrane. Yet, this pathway might not be the predominant diffusion pathway of gases and does not account for constrictive pores; not all molecules will be affected by the microstructure to the same extent when migrating through such a layer. The inherent difference of the mean free path between each gaseous species leads to different Knudsen numbers and thus, different diffusion pathways for different species at different temperatures, gas compositions and transport regimes. It can thus be inferred that, for each species and each transport regime, a different tortuosity value is dominating.
Moreover, in experimental approaches, tortuosity is not always presented explicitly, but is rather combined with porosity into a 'diffusibility' [17,18] or 'effective relative diffusivity' [19,20] value expressed as 1/t 2 . Additionally, in the field of battery research, tortuosity is contained in the MacMullin number N M , which relates the bulk conductivity of the electrolyte σ bulk to the effective conductivity of the porous electrolyte σ eff [21][22][23][24].
These different definitions and applications of tortuosity cause differences in its interpretation and calculation approach. For example, geometric-based tortuosity takes only the shortest path length into account, while flux-based values rather account for the path of least resistance. Hence, the resulting values differ appreciably. These discrepancies are reflected by the vast number of different tortuosity calculation approaches shown here.

Porosity-tortuosity relationships
Employing a porosity-tortuosity relationship is one of the most fundamental and straightforward approaches to derive a tortuosity (or effective medium property) of a porous structure. Such relationships, of theoretical or empirical origin, directly calculate a tortuosity value solely based on a porosity of a sample.
In the comprehensive work by Shen and Chen [25], a review of past and present correlations is provided, among which the Bruggeman equation is the most well-known and most widespread relation in the field of electrochemistry [26]. Eq. (6) presents the generally used form of the Bruggeman relationship, where α is the Bruggeman exponent which, in its standard form, is considered to be 1.5. Recently, the authors have provided a translation and explanation of the mathematical formulation of Bruggeman which is used to derive the widely used model [27].
While the history of the Bruggeman correlation can be traced back to the 1930s, its proliferation was not notable until the 1950s: Hoogschagen was one of the first to use the Bruggeman and Maxwell relation (cf. Eq. (7)) to validate experiments, where gas diffusion through glass spheres was measured. He observed that values for the labyrinth factor (1/t 2 ) lay between the Maxwell and Bruggeman correlation, but slightly closer to the latter [17].
De La Rue and Tobias achieved similar results when measuring the effective conductivity values of liquid ZnBr 2 electrolyte solution. A variety of non-conducting glass spheres of different sizes were embedded into the electrolyte to achieve different volume fractions. The conductivity as a function of volume fraction of the embedded phase was evaluated. As was the case in Hoogschagen's publication [17], results lay between the Maxwell [28,29] and Bruggeman relation [30]. Since then, the Bruggeman equation has become a commonly used method to derive effective medium properties of porous structures in batteries [31][32][33][34][35] and proton exchange membrane (PEM) fuel cells [36][37][38][39][40][41][42][43][44][45]. Moreover, it has been implemented as a standard addition to predicting microstructures in electrochemistry models, such as in the COMSOL Multiphysics modelling software (COMSOL, Inc.) [24].
However, the predictions given by the Bruggeman correlation are not always consistent with experimental results [24,46]. As a consequence, researchers have adjusted the Bruggeman equation by altering the exponent α to fit experimental values. Thorat et al. [47]. even included an additional scaling factor γ to correlate the Bruggeman model with their experiments, resulting in Eq. (6) to be extended to the following form: Thorat et al. used AC impedance spectroscopy and the polarisation-interrupt method (cf. section 'Electrochemical experiments') to extract the tortuosity of a battery separator (Celgard 2400) and cathode samples (LiFePO 4 and LiCoO 2 ). Tortuosity values of the battery cathode samples were plotted as a function of porosity, and an exponential fitting curve was superimposed. The exponent of the fitting curve amounted to −0.53, which is equivalent to a Bruggeman exponent of 1.53 and thus, very close to its derived value. However, achieved tortuosities were almost twice as high as predicted by the standard Bruggeman relationship, which is why a scaling parameter γ amounting to 1.8 was introduced. This approach of adjusting α and γ was widely adopted, showing that depending on the analysed structure, both parameters can deviate from the ideal values of 1 and 1.5, respectively [32,[47][48][49][50][51][52][53][54][55].
A further refinement of this approach was realised by Zacharias et al. [53], who made α and γ a function of their battery electrode composition. For this, the dry weight fractions of graphite, carbon black and polyvinylidene fluoride were considered, resulting in higher γ values (2.5 and 2.6) and lower α values (1.27 and 1.28) compared to values from Thorat et al. [47]. Table 1 and Figure 2 compare several derived Bruggeman exponents and scaling parameters for different porous materials for battery applications. These were each extracted as a function of several experimental measurement points and used to extrapolate the presented curves as function of porosity. It is notable that even for this small class of materials, values for α and γ differ significantly from each other. The differences in manufacturing techniques, and also the differences of composition, pore size distribution and other microstructural characteristics of each battery layer, contribute to such a large spread of values. Some of these derivations, however, achieve tortuosity values below unity when extrapolated to high porosity values, which is in contradiction to the definition and physical significance of τ. Moreover, a porosity of 100% necessitates a tortuosity of unity, yet, this is not achieved by all correlations. Both of these findings cast doubts on the usefulness of this method. As a consequence, the application and interpretation of α and γ values have to be analysed with caution.
Hence, evaluating the validity of the Bruggeman correlation is still an ongoing field of research. Chung et al. [56] used X-ray computed tomography and simulation techniques for an extensive study to evaluate the effect of battery membrane fabrication and processing methods on the tortuosity. In total, 16 LiNi 1/3 Mn 1/ extracted by simulating mass transport according to Fick's law across the sample volume (see section 'Flux based algorithms'). It was shown that calculated tortuosity values always lie slightly above the Bruggeman correlation. For further investigation, samples based on the particle size distributions of the imaged samples were computer generated, for which the orientation and particle packing was varied. It was discovered that perfectly ordered particle distributions result in tortuosities close to the Bruggeman relationship throughout the range of porosity values [56].
Continuing work in the field of battery research from Wood and co-workers (cf. [54,56,57]) culminated in the development of an open source programme called BruggemanEstimator [58]. This programme allows the extraction of the Bruggeman exponent α in each dimension of a 3D sample volume using two 2D images, namely one top view and one cross-sectional view. The Bruggeman exponent of the sample is achieved by applying the differential effective medium approximation method introduced by Bruggeman. In comparison to previously obtained values, the results calculated by the BruggemanEstimator software agreed well with numerical tortuosity calculation methods [58] and have been recently applied in practice [59]. This approach is similar to stereological methods which quantify 3D properties based in 2D image slices [60]. The advantage of stereology is the reduced experimental efforts necessary to extract results. However, Taiwo et al. [2] recently concluded that values based on stereological approaches may deviate appreciably from 3D measurements.
Moreover, a wide range of recent studies report conflicting results on the validity of the Bruggeman correlation when compared to the calculations done using tomography techniques. Conclusions vary substantially as in some instances, simulations agree well with the Bruggeman correlation [34,56], while considerable disagreement was observed in other cases [45,[61][62][63]. The reason for this seems to be sample specific, as heterogeneity and geometry are characteristics of porous materials that are not accounted for by the Bruggeman correlation. The aforementioned studies have shown that the characteristic shape of the analysed microstructure has considerable effects on the validity of the Bruggeman relation: spherical structures, which follow Bruggeman's initial hypothesis very closely, adhere to the correlation. The correlation, however, is less suitable for connected solid phases and complex porous networks. This is further complicated by the distinctions (or lack thereof) between geometrical and transport limiting tortuosity [64]. Moreover, porosity-tortuosity relationships provide limited information in areas, where the analysed sample consists of several layers with different microstructural features, such as multilayer battery separators [52]. These combine different properties into a single separator; i.e. each individual layer exhibits distinct structural properties, and for this reason, the simplified assumption of a homogenous sample volume made by the Bruggeman correlation is no longer valid. As a conclusion, it can be stated that porosity-tortuosity relationships are only applicable and reliable when executed across homogeneous microstructures which are similar to the microstructure used to derive the respective relationship.

Experimentally derived tortuosity
Historically, the lack of detailed geometrical information on complex porous media in 3D has limited the ability of researchers to extract meaningful data on the tortuosity of a porous body. In the absence of this information, effective transport properties of porous structures have been derived experimentally by means of diffusion cell experiments [17,18,[65][66][67][68][69][70][71] and electrochemical measurements [16,20,47,72].

Diffusion cell experiments
As reviewed by He et al. [73], diffusion measurement methods in the field of fuel cell research aim at extracting effective diffusion coefficients for distinct gas mixtures. Typically, a porous sample is mounted between an upper and a lower gas channel where two different gases enter the upper and lower chamber. Owing to the concentration gradient across the porous material, diffusion of either gas to the opposite channel is induced.
Measuring the concentration of either gas in both streams allows the calculation of the diffusion fluxes across the membrane via a mass balance over the cell, as illustrated in Figure 3. The effective binary diffusion coefficient and, in turn, the tortuosity of the sample, are subsequently derived by applying a suitable diffusion model. The applicability of such models is dependent on the diffusion mechanism taking place through the porous medium (ordinary, Knudsen and/or viscous diffusion cf. [15,74,75]). The most direct diffusion model considering ordinary and Knudsen diffusion is the Fick model. For this, Fick's law is extended by combining the Knudsen diffusion coefficient D i,K , (Eq. (9)) with the effective binary diffusion coefficient D eff ij via the Bosanquet equation (Eq. (10)). Moreover, Fick's model is capable of incorporating viscous flux via Darcy's law, as presented by Eq. (11). In this case, the model is referred to as advective-diffusion model [76].
In the following equations, r P is the mean pore radius, R is the ideal gas constant, T is the temperature, M is the molar mass, D eff i is the effective diffusion coefficient of a species, B O is the viscous flow parameter and μ is the dynamic viscosity and p is the pressure  Table 1. and w is the mass fraction.
While Fick's law assumes equimolar diffusion in a binary gas mixture, Mills [77] suggested that diffusion follows equimass principles. By converting the molar concentration gradient of Fick's first law into a gradient of mass fraction, the governing equation for equimass diffusion was achieved: Again, the Bosanquet equation and Darcy's law can be used to extend the formulation to cater for additional diffusion mechanisms besides ordinary diffusion. The effective diffusion coefficients are directly achieved via the above models which make their application straightforward. More complex models, such as the Dusty Gas Model (DGM) [74] or Maxwell-Stefan Model (MSM), shown in Eq. (13) and Eq. (14), respectively, combine several diffusion modes a priori. The DGM, for example, includes expressions for ordinary, Knudsen and viscous flux, where x is the molar fraction: The MSM uses the same correlation as the DGM, but neglects Knudsen diffusion effects. As a result, coefficients related to Knudsen diffusion drop out of equation Eq. (13) and result in the following formulation: The accuracy of these models has been discussed [78] and evaluated in literature, predominantly by comparing them to measured concentration polarisation losses in SOFC anodes [79,80]. In these experiments, the DGM achieved highest accuracy among the analysed models, which might be one reason for its widespread use in literature [80][81][82][83][84], while the simplicity but lower accuracy of the Fick model was frequently highlighted. However, in these cases, tortuosity is usually used as a fitting parameter to tailor calculation results to measured data. Consequently, the extracted tortuosity values are highly dependent on the accuracy of the applied model.
One recent application of extracting the tortuosity of a porous sample via diffusion cell experiments was presented by Tjaden et al. [85]. In their analysis, a variety of binary gas mixtures were tested on a planar YSZ porous support layer of an oxygen transport membrane. It was shown that the tortuosity for different binary gas mixtures is not a constant value, but depends on the involved gaseous species. Owing to the difference in mean free path between the different constituents, different diffusion pathways dominate and thus, different tortuosity values are achieved. For example, at ambient temperature, tortuosity based on the diffusive flux of CH 4 in the CH 4 -N 2 binary gas mixture amounted to approximately 2.3, while tortuosity of the N 2 diffusion flux amounted to approximately 2.5. In addition, the authors showed that tortuosity increased with increasing temperature: the average tortuosity of all binary gas mixtures increased from 2.36 to 2.73.
The contrary temperature-dependent behaviour of measured diffusion coefficients was observed by Zamel et al. [18]. A Loschmidt cell was used to measure the effective diffusion coefficient of a O 2 -N 2 gas mixture migrating through carbon paper, which is commonly applied as the gas diffusion layer (GDL) in polymer electrolyte membrane (PEM) fuel cells. When increasing the temperature from 25°C to 80°C, the bulk diffusion coefficient of the gas mixture, achieved via a resistance network model based on Fick's law, increased from approximately 0.2 cm 2 s −1 to 0.275 cm 2 s −1 , while the effective diffusion coefficient increased from approximately 0.05 cm 2 s −1 to 0.075 cm 2 s −1 . This causes the factor 1/t 2 to increase by approximately 11.5% from 0.252 to 0.281. Thus, when considering a constant porosity value, the tortuosity decreases to the same extent. In addition, the authors compared the calculated diffusibility values to a set of porosity-tortuosity relationships, among others, the Bruggeman relation. In all cases, these relationships overestimate the effective diffusion coefficient (and thus, underestimate tortuosity), which was also observed by Tjaden et al. [85].
The discrepancies in temperature dependence might be caused by the fundamental difference in microstructural aspects: while YSZ-based porous structures are tailored to feature connected solid phases with specific porosities, carbon paper-based gas diffusion layers are an accumulation of randomly oriented, fine fibres, typically with much higher porosity. The resulting differences in pore size distribution, mean pore diameter and porosity cause the observed diffusion mechanisms to differ significantly among such samples.

Electrochemical experiments
Mass transport limitations play a vital role in electrochemical devices as they are responsible for concentration polarisation at high current densities. For example, as current densities increase, the fuel demand in a fuel cell increases linearly, as shown in Eq. (15), whereṅ fuel is the molar flow rate of fuel gas, i is the current density, A is the area, n is the equivalent electrons per mole of reactant and F is the Faraday constant: [7] The fuel consumption rates at the active sites of a fuel cell are limited by the maximum diffusion rate of reactant through the porous structures. As introduced in previous sections, diffusive mass transport and, as such, mass transport limitations are a function of the complex microstructure of the involved porous membrane layers. Hence, microstructural parameters, such as tortuosity, are achievable by measuring concentration losses of fuel cells and applying gas diffusion theory.
In this respect, SOFCs offer the possibility to investigate the effect of fuel gas compositions on the performance owing to their wide fuel flexibility. A thorough study of this topic was presented by Jiang and Virkar [16]. As the effects of mass transport limitations are dominating under high current density operations, Jiang and Virkar modified Fick's law to express the effective diffusion coefficient as a function of the limiting current density of the fuel cell under specific operating conditions. The resulting equation is presented thereafter, where i lim is the limiting current density, p 0 fuel is the partial pressure of fuel at the gas inlet and d is the thickness: The limiting current density was measured experimentally from polarisation curves for a set of binary and ternary fuel gas mixtures, including H 2 73. Moreover, in direct comparison between the two binary gas mixtures, it was revealed that fuel cell performance was higher using H 2 as fuel gas rather than CO which, besides the lower electrochemical activity of CO, was due to the significantly faster diffusion rate of H 2 . These results confirm the findings of different tortuosity values for different binary gas mixtures presented in the previous chapter.
Brus et al. [72] adopted the same methodology to compare electrochemically derived tortuosity values with an image-based tortuosity calculation method, namely the random walk method (cf. section 'Voxel-based calculation methods'). For their experiments, a button-type SOFC sample was manufactured to measure impedance spectra and polarisation characteristics at 700°C and 800°C. This way, the limiting current densities were extracted for H 2 concentrations between 2.5 and 90% in N 2 and inserted into Jiang and Virkar's model. After these experiments, the 3D microstructure of the anode was reconstructed using focused ion beam-scanning electron microscope (FIB-SEM) tomography, and the random walk method was executed. For each hydrogen concentration and for both operating temperatures, a distinct tortuosity value is calculated, whereas the random walk method results only in a single value, as shown in Figure 4. Here, only the tortuosity values calculated for low hydrogen concentrations and, as such, high concentration polarisation were considered as accurate representative values. In these cases, the experimentally derived tortuosities agreed well with the random walk value. Hence, under standard fuel cell operating regimes, where activation and Ohmic losses dominate, concentration losses, and thus the tortuosity of the porous layers affect the performance only slightly.
However, experimental-based tortuosity values are valid only for the specific experiment at hand. While the results between image-and experimental-based tortuosity values in the above case are close, this agreement might not be reproducible when the fuel gas composition changes. Figure 4 also shows that higher temperatures have a positive effect on tortuosity: for each fuel gas composition, the tortuosity is lower at higher temperature, which can be explained by the higher catalytic activity and faster diffusion rate. Yet, aside of the effect of temperature, the influence of structural parameters, such as the layer thickness on the tortuosity of SOFC anodes, is of interest. This was investigated by Tsai and Schmidt [86][87][88], who, again, applied Jiang and Virkar's approach for this purpose. While they observed the same dependency of tortuosity on H 2 concentration as Brus et al. [72], Tsai and Schmidt [88] showed that electrode thickness had no effect on the achieved tortuosity values, which is expected for steady state operation.
Electrochemical experiments have also been applied to study microstructures of lithium-ion battery materials. Thorat et al. [47] used polarisation-interrupt (or restricted diffusion) experiments [89][90][91] and impedance spectroscopy to measure the tortuosity in electrode and separator layers. Using the polarisation interrupt technique, Thorat et al. derived the tortuosity of two distinct active material films consisting of LiFePO 4 and LiCoO 2 , respectively. On the other hand, AC impedance spectroscopy was carried out to determine the effective conductivity of the electrolyte in the separator and ultimately, the MacMullin number or the tortuosity of the separator itself. While the authors used the AC impedance experiments to validate the polarisation interrupt experiments, the effect of porosity on the tortuosity of the active material films was in the centre of their research and led to the adjustment of the scaling factor to 1.8 and Bruggeman exponent α to 1.53 [47] as discussed in section 'Porosity-tortuosity relationships'.
With the development of advanced manufacturing techniques, lithium ion battery electrode microstructures can be tailored and optimised to meet userand application-specific demands. Bae et al. [92], for example, applied a two-pronged approach to improve electrode design: first, using a modified model by Doyle and Newman [93], the tortuosity of different electrode microstructures with periodically spaced flow channels was calculated. Based on these results, LiCoO 2 electrodes mimicking the modelled microstructures were manufactured using a co-extrusion procedure. In their model, electrodes with flow channel spacing equal to or smaller than the electrode thickness offered lowest tortuosity values. To validate these findings, charge and discharge curves of the manufactured samples with large, medium and small channel spacing were measured. As predicted, the sample with finest and most closely spaced channels yielded highest specific capacity of approximately 8 mAhcm −2 at C-rates of one and two. The authors attributed this improved capacity to the lower tortuosity of their manufactured electrode, validating their model.
In general, experimental set-ups can be adjusted to fit the operating conditions of the analysed specimen. However, as the derived results are fitting parameters, the tortuosity values are highly dependent on the applied model. Moreover, while fuel cell experiments can be highly versatile in terms of operating temperature and applied fuel gas, batteries are not subject to such variations. Hence, it appears to be easier to extract an overall valid tortuosity value for a battery layer than a fuel cell layer.

Tortuosity calculation in 3D volumes
The advent of sophisticated and easily accessible tomography methods has increased the amount of obtainable data of porous samples which fundamentally changed the perception of microstructural characterisation in 3D [94]. FIB-SEM slice and view tomography [95], and X-ray computed tomography (X-ray CT) [96,97] are among the most prominent methods of reconstructing a sample in three dimensions. Even though the operation and image acquisition of both methods are radically different, comparative studies showed that acquired data are identical when the resolution is the same [85,98,99].
FIB-SEM tomography is commonly used in SOFC research to achieve solid phase segmentation between the YSZ and Ni-phase in the anode [100]. Moreover, this technique offers the possibly of using additional analysis methods, such as energy dispersive spectroscopy, electron backscatter diffraction and secondary ion mass spectroscopy [1]. Hence, 3D chemical and crystallographic data can be acquired in parallel with microstructural parameters. The advantage of X-ray tomography lies in the non-destructive nature of operation which allows in-situ and in-operando analyses. Changes in microstructure as function of continuing and varying operating conditions can be tracked [101]. This allows a four-dimensional analysis of microstructural evolution during processing, redox reactions and degradation [1].
In recent years, tomographic reconstruction of microstructures in electrochemical devices, illustrated in Figure 5, has become increasingly widespread, offering the possibility to evaluate vital parameters, such as triple phase boundary length in SOFCs [102,103], connectivity [104], phase distribution [105] and tortuosity [81,102] at different length scales [106]. Additionally, the effect of microstructural parameters on the performance of electrochemical devices has been evaluated by generating synthetic 3D volumes in silico [34,56,107,108]. The purpose for this process is to directly evaluate the effect of specific microstructural variations, such as porosity, pore size distribution, shape or packing orientation of particles on mass transport.
There remains some confusion in the literature regarding the different definitions of tortuosity for the purpose of image-based modelling: here, we distinguish between two main approaches in extracting tortuosity: (1) Geometric-based algorithms, which aim to determine the shortest path length through a porous structure by purely considering geometric aspects. (2) Flux-based algorithms, which mimic mass transport and diffusion behaviour, which is not taken into consideration in geometric-based algorithms. These methods are further divided into the following two subsections: (a) Voxel-based algorithms that take the extracted dataset and directly execute tortuosity extraction techniques across the voxel domain of the analysed phase. (b) Mesh-based approaches which rely on generating a volume mesh of the analysed phase to prepare the sample for computational fluid dynamics (CFD) programmes.
It is evident that the increase in the development of such techniques correlates with the increasing accessibility of tomography equipment and high-performance computers.

Geometric-based algorithms
Geometric algorithms are commonly used to find the shortest pathway through a porous structure and thus, its tortuosity. The pore centroid method [63, Figure 5. Illustration demonstrating that high-resolution tomography is necessary to extract microstructural features which affect diffusive mass transport. 106,[109][110][111], the fast marching method (FMM) [2,112,113], the distance propagation method [114], as well as or other shortest path search methods [115,116] achieve this by being executed on the voxel domain of the analysed phase. These methods are straightforward in their application, as mesh preparation and refinement is not required. In addition, the results directly follow the initial definition of tortuosity, making them conceptually easier to interpret. Furthermore, apart from the pore centroid method, these algorithms create a distance map, which incorporates the distance of each pixel to the starting plane of the algorithm. Using the resulting distance map allows not only the identification of the shortest pathway, but also the generation of a tortuosity histogram.
The FMM achieves this by simulating an advancing front, starting from one plane of the sample towards the opposite plane in the considered phase. The algorithm measures the time it takes for the front to reach each pixel on its way. By knowing the speed of advance of the front and the time it takes to arrive at a pixel, the distance between each pixel and the starting plane is achieved and tabulated in a distance map. Finally, tortuosity is calculated by dividing the shortest path length between two opposing planes by the Euclidean distance of the two endpoints of that path.
Jørgensen et al. [112] exploited the FMM-based tortuosity histograms of a strontium-substituted lanthanum cobaltite (LSC) and gadolinium-substituted ceria (CGO) SOFC cathode to understand microstructural characteristics of each phase. In accordance with each phase's volume fraction, LSC features higher tortuosity values than CGO. The distinct shapes and specifics of each phase's tortuosity achieved by the FMM-based histograms are able to provide more insight into the microstructural build-up of a sample compared to a single, mean tortuosity value.
Yet, tortuosity histograms do not show where the specific high or low tortuosity values are located within the sample. This was realised by Chen-Wiegart et al. [114], who combined different tomography methods and distance propagation-based tortuosity calculation approaches on various samples. Specimens included, among others, a LiCoO 2 battery cathode, which was reconstructed using an X-ray tomography. Geometric tortuosity values were then achieved by pixel counting and distance-measuring techniques. The resulting values were presented as 3D distribution across the battery cathode sample, as shown in Figure 6. The local variation in the image slices ranges from 1 to 2.5, which can also be ascertained from the tortuosity histogram. However, as tortuosity poses a resistance to mass and charge transport, the local tortuosity distribution is capable of pinpointing areas of low reactivity. It can be used to explain regions of increased charge transfer, areas of low fuel conversion, uneven charging or catalyst utilisation and degradation. Shearing et al. [115] extended the approach of spatial distribution of geometric tortuosity to include additional characteristics such as volume specific surface area and porosity. A reconstructed graphite Li-ion battery electrode was segmented into a mosaic of equally sized volumes and for each tile, the aforementioned parameters were calculated and visualised to highlight the relation between them. While in most cases, tiles with high porosity featured low tortuosity, some sub-volumes exhibited low tortuosity coupled with low porosity. Even though this combination seems counterintuitive, it emphasises the complex interrelation between different microstructural parameters which are not always as clear as expected.
For comparative purposes, Chen-Wiegart et al. executed a diffusion simulation analogous to the one used in [4] (cf. section 'Mesh-based calculation methods') across the same sample volumes. It was shown that the results between the distance propagation and diffusion method of the pore phase in the LiCoO 2 sample agreed well. However, when applying the same calculation approaches to two SOFC samples, the geometrically derived tortuosity values for the pore and YSZ phases were consistently below the diffusion-based tortuosity methods. The difference might stem from the inherent difference between geometric-and diffusionbased considerations: the geometrically shortest path through a structure is not always the path of least resistance for a flux, owing to the presence of constrictions and pore necks. Further discussion on the differences of these considerations is presented in section 'Flux-based algorithms'.
In contrast to the aforementioned algorithms, the pore centroid method does not provide a histogram of tortuosity values or spatial distribution of tortuosity, but rather arrives at one specific value of tortuosity along each dimension of a sample. The calculation algorithm follows the centre of mass of a phase of a 2D plane along the third axis of the volume. The length of the pathway going through each centroid is then calculated and used to determine the tortuosity as depicted in Figure 7. Despite its shortcomings in comparison with the previous algorithms, the pore centroid is a standard option in image and volume processing programmes such as Amira and Avizo (both FEI). As such, it is easily applied for comparative studies and capable of giving a quick tortuosity estimate.
Cooper et al. used the pore centroid method for comparison reasons when studying an SOFC [111] and a battery electrode (cf. section 'Mesh-based calculation methods') [63]. In [111], the tortuosity of the solid and pore phase of an LSCF SOFC cathode was determined by a variety of calculation algorithms, namely heat flux simulation (cf. section 'Mesh-based calculation methods'), Avizo XLab plugin, diffusion simulation, random walk method (cf. section 'Voxelbased calculation methods') and pore centroid method. These algorithms were executed across the same sample after imaging at 14°C and 695°C using synchrotron X-ray nano CT which have previously been extracted by Shearing et al. [117]. The pore centroid method produced the lowest tortuosity values for both phases at both temperatures and closely followed the Bruggeman relationship. Yet, the flux-based calculation algorithms agreed well with each other as values lay between the heat flux simulation and the random walk method. The average tortuosity for the pore phase amounted to approximately 1.21 in all three dimensions at both temperatures and lay visibly below the values reported by Gostovic et al. [109] using the same method. Large variability in homogeneity of a sample significantly affects the results achieved by the pore centroid method, causing visible fluctuations. In this respect, Cooper [118] pointed out that if the analysed characteristic feature becomes small compared to the control volume, the centroid of each 2D plane will tend towards the centre, resulting in a tortuosity of unity which casts doubt on the applicability of this approach.

Flux-based algorithms
Even though geometric-based tortuosity calculation algorithms can extract useful data concerning the distribution of geometric tortuosity across a sample, these algorithms do not mimic the flux-like behaviour of transport phenomena. For example, small connections consisting only of one voxel would contribute only a negligible amount to the overall flux of transported species, while they are fully included in the above calculation methods. As a result, flux-based algorithms focus on simulating the transport mechanism at hand to extract the tortuosity of a sample. Here, this method is separated into two parts, namely voxeland mesh-based calculation approaches.

Voxel-based calculation methods
Voxel-based algorithms are directly executed across the voxel domain of the reconstructed volume. This means that for the methods introduced below, no additional re-tessellation or re-meshing steps are necessary after the sample has been segmented. In most cases, a binarised 2D image sequence is sufficient to operate the calculation procedure.
One of the first applications of combining X-ray nano tomography with image-based tortuosity calculation was presented by Izzo et al. [81], where X-ray CT was used to gather microstructural parameters of a porous SOFC anode, including porosity, tortuosity and pore size distribution. The authors solved the Laplace equation of diffusive mass transport through the pore phase of the electrode as explained in a different publication from the group [119]. Grew et al. [120] applied the same methodology, but extended its application to the solid phases of a Ni-YSZ SOFC anode. As effective ionic conductivity and electronic conductivity are affected by the tortuous nature of fuel cell electrode layers (cf. equation 3), tortuosities of solid phases are equally as important as of pore phases. Yet, they were at least a factor of 1.2 higher.
Their work was further refined in [121] by calculating the representative volume element (RVE) of the pore phase tortuosity by solving the Laplace equation using the same method. Cooper [118] programmed a MATLAB (Mathworks) Laplace solver called TauFactor [122][123][124] to extract the tortuosity of a two-phase segmented 3D tiff stack as shown in Figure 8. The solver then determines the tortuosity in each dimension for both phases. In [118], Cooper compared the results of the TauFactor solver to previous work presented in [111], revealing that the solver gives similar results as the Avizo package XLab Thermo and the heat flux simulation.
Aside from solving the Laplace equation to arrive at the tortuosity of their sample, Izzo et al. [81] included the lattice Boltzmann method (LBM) [125] to model multi-component gas transport coupled with an electrochemical model to visualise the H 2 distribution in the anode. Owing to the capability to model gaseous, ionic and electronic transport, the LBM became widely applied in fuel cell research, also with the focus of extracting tortuosity in different phases of a functional layer [103,[126][127][128][129][130][131]. For this, the LBM uses the particle distribution function (PDF) f i a PDF , which is a function describing the probability of encountering a particle of a species i at a certain location x PDF with a certain speed e i a PDF at a certain point of time t moving in a certain direction α PDF .
The LBM consists of two steps, namely streaming and collision, which are carried out on each point of a lattice: during streaming, the particles migrate to adjacent lattice points while during collision, the interactions between particles at each lattice point governed by the collision term V i a PDF are computed. Both steps are collectively expressed by the lattice Boltzmann equation: [119,132] f i a PDF (x PDF + e i a PDF , t + 1) − f i a PDF (x PDF , t) = V i a PDF (17) Using this approach, Iwai et al. [103] arrived at tortuosity values for each phase in the porous Ni-YSZ anode by calculating the effective diffusion coefficient and effective ionic as well as electronic conductivities of the respective phases. The anode sample was reconstructed based on FIB-SEM tomography, where the Ni and YSZ phases were identified via EDX mapping, to correlate the correct phase to the respective electron image. Table 2 compares the achieved tortuosity values for all three phases along each dimension using the LBM as well as the random walk method, which is introduced thereafter. It is evident that the tortuosity values of the solid phases are higher compared to those of the pore phase, which is identical to the findings presented by Chen-Wiegart et al. [114] using a distance mapping approach. Nevertheless, values for the pore phase tortuosity are lower than but comparable to the values found by Izzo et al. [81] However, owing to the observed directional anisotropy of the solid phase tortuosities, Iwai et al. concluded that the sample volume was not sufficiently large to present effective ionic and electronic conductivity values. Vivet et al. [133] achieved similarly high Ni-phase tortuosity values using a finite difference method. However, owing to the higher YSZ fraction in their sample, achieved YSZ tortuosities lay below the values reported by Iwai et al. [103].
The aforementioned random walk method [29,103,[134][135][136][137] mimics a diffusion process by distributing a number of non-sorbing particles, so-called walkers, across the segmented voxel phase. The algorithm then starts a time-step sequence, where at each step, every walker choses one neighbouring voxel as its next location. If that neighbouring voxel is of the same phase (e.g. pore phase), the walker migrates to that new location. However, if the chosen neighbouring voxel is of a different phase (e.g. solid phase), the walker remains at its current location and chooses a Figure 8. Results of the TauFactor solver by Cooper [118] running across the pore phase of a porous sample showing the binary image map, the initial, linear concentration distribution and the concentration distribution at steady state. different neighbouring voxel at the following time step. By repeating this sequence, the mean square displacement r 2 (t) of the walkers in the analysed phase is calculated which, in turn, is used to achieve an effective diffusion coefficient D eff , where V phase is the volume fraction of the analysed phase: Tortuosity is then calculated by comparing the effective diffusion coefficient to the bulk diffusion coefficient through an empty volume of equal dimensions. The random walk approach was first formulated in the 1990s [5,138,139] and found its way into electrochemistry via Kishimoto et al. [134], after having been used to extract the tortuosity of porous rocks [140]. However, the obtained tortuosity is affected by the number of walkers and by the number of time steps chosen for the calculation. This is why, in [103], 100,000 walkers and 10,000,000 time steps are chosen to ensure high accuracy of the results (cf. Table 2).
A similar comparative study using the random walk method was carried out by Tariq et al. [137]. The tortuosity values of a Li-ion battery anode, calculated by the random walk method, were compared to the results based on a sub-grid scale finite volume method explained by Kishimoto et al. [141]. As shown in Table 3, results for both methods agree excellently, revealing a higher tortuosity along the z-axis of the pore phase. The authors noted that a RVE analysis would reveal, if this anisotropy was persistent or if the high value was caused by a local heterogeneity. Yet, it was explained that the computation time needed for the random walk method is only a fraction compared to the finite volume method.
Mesh-based calculation methods By applying the same tomography methods mentioned in the previous section, extracted datasets can be represented as volume meshes for additional analysis algorithms enabled, for example, by CFD or finite element software packages. These programmes allow the simulation of heat, mass and/or charge transport through the generated mesh of the investigated structure to subsequently evaluate the tortuosity. In the data preparation process, parameters chosen for sample smoothing, surface repair and mesh generation affect mesh quality and thus the simulation results. Hence, care must be taken when choosing these parameters [85] and sensitivity analyses should be carried out to verify the consistency of the chosen values.
Pioneering work in this field was realised by Wilson et al. [102], who reconstructed an SOFC anode using FIB-SEM tomography. The tortuosity of the pore phase was then extracted to assess the mass transport limitations at high current densities. For this, the sample volume was converted into a finite element mesh to solve the Laplace equation in FEMLAB (now COMSOL Multiphysics).
Extensive simulation work in the field of electrochemical devices, using a similar approach as presented above, has been carried out by Ivers-Tiffée and coworkers: initially based on COMSOL Multiphysics, the group developed the 3D finite element tool Par-Cell3D to model the behaviour of fuel cells [142][143][144][145] and batteries [146]. Joos et al. [147] used this tool to investigate the RVE of tortuosity of an SOFC cathode for both phases, namely the pore and the mixed ionic-electronic conducting LSCF phase. In total, the RVE of porosity, volume-specific surface area and tortuosity were calculated for three separate volumes, of which the latter one is presented in Figure 9. The results for both phases in sample volumes 1 and 3 agree excellently with each other, achieving a flat development for electrode thicknesses of l cat > 10 μm.  Table 3. Tortuosity values for graphite and pore phase using the random walk method and finite volume method [137]. Random Figure 9. RVE analysis of the tortuosity factor for the pore and LSCF phase of an SOFC cathode as function of electrode thickness reproduced with permission from Elsevier [147].
However, the tortuosity of the LSCF phase in sample two takes an electrode thickness almost twice as long as for the other sample volumes to produce a flat curve. To follow the nomenclature of this review, it has to be pointed out that τ in Figure 9 ought to be replaced by τ [2]. Besides COMSOL Multiphysics [148,149], researchers have calculated tortuosity using programmes, such as Cast3M [150] or custom made models, which focus on a specific electrochemical device, such as Batts3d [34,56,151].
In addition to simulating mass and charge transport, the tortuosity is also computable by exploiting the mathematical similarity between Fourier's law of heat conduction and Fick's law of diffusion shown in Eq. (19) and Eq. (20), where J eff is the effective diffusion flux, c is the molar concentration,q eff is the effective heat flux, λ bulk is the bulk thermal conductivity [85,106,[152][153][154]. By rearranging the temperature gradient, the heat flux and thus, the tortuosity can be calculated along each axis of a sample. Cooper et al. [63] scanned a commercially available LiFePO 4 battery cathode using an X-ray synchrotron nano CT and investigated the tortuosity of the pore phase using heat flux simulation. A cube of 8.8 μm side length was cropped and meshed using an adaptive polyhedral volume mesh. The heat flux across the porous phase of the sample was simulated in StarCCM+ (CD-adapco), resulting in a temperature distribution across the analysed volume (see Figure 10), where the temperature of each mesh element can be understood as a concentration value of a migrating species.
For further analysis, the sample was divided into eight non-overlapping sub-samples to compare the pore centroid and heat flux-based tortuosity factors. To extract a single tortuosity value for each subvolume, the concept of the characteristic tortuosity τ C was introduced: [63,85] The geometric-based characteristic tortuosity factors (κ geo ) were then plotted as a function of the fluxbased characteristic tortuosity factors (κ flux ). The resulting graph revealed a high degree of correlation following Eq. (22). This equation reveals that the simulation-based tortuosity factor (and as such, tortuosity) is always higher than the geometric-based value of the same sample in cases where either value is >1. As mentioned previously, geometric tortuosity algorithms do not take the effect of pore constrictions into account, which would affect a transport flux. This means that any connection consisting of only one voxel in diameter is fully considered in the calculation process, while from a flux point of view, such a pore would not allow a significant amount of mass/charge to pass through. As a result, geometric-based tortuosity values tend to be visibly lower compared to simulationbased values [63,85]. k geo = 0.5 ln (k flux ) + 1 (22) It is common practice to subdivide a given sample volume into an array of smaller sub-samples and extract the tortuosity for each individually [51,147,155]. Although non-trivial, this approach allows researchers to extract similar results as tortuosity histograms and tortuosity distribution maps (cf. section 'Geometric-based algorithms') using flux-based methods. This approach reveals the homogeneity or heterogeneity of a sample, comparable to tortuosity histograms, and pinpoints the locations of high or low tortuosity. Kehrwald et al. [51] were among the first ones to apply this methodology on a battery electrode by solving Fick's law using the programme Star-CD (CD-adapco) on a total of 12 sub-volumes. Local tortuosities showed differences of a factor of three, which might lead to inefficiencies during charging and discharging of the battery: Li + ions will avoid areas with higher tortuosity, but seek areas with low tortuosity, highlighting the need for homogeneous microstructures in this field to prevent spatial distribution of performance and degradation [155]. In addition, microstructural inhomogeneities might be the cause of failure mechanisms and material fractures [51].  Table 4 and Table 5 list pore phase tortuosity factors and tortuosity values for different image-based calculation methods along all three axes of porous samples. Calculation approaches, which take flux-like behaviour into account (cf. Table 4), arrive at higher tortuosity values compared to purely geometric-based algorithms (cf. Table 5). It is thus imperative to distinguish between these two approaches as otherwise, misinterpretation may ensue. Furthermore, differences in tortuosity values are observed even when analysing the same type of samples. This can be explained by the chosen imaging resolution; higher resolution uncovers smaller pore structures and improves pore connectivity. However, the effect of varying the pixel size can differ substantially, depending on the analysed microstructural parameter. Shearing et al. [106] imaged the same commercially available Li-ion battery cathode using different X-ray systems: a laboratory-based micro CT system and a laboratory-based nano CT system. The pixel size for each system amounted to 597 and 65 nm and the extracted volumes amounted to 10,434,731 μm 3 and 75,164 μm 3 , respectively. The authors calculated the porosity, geometric tortuosity using the pore centroid method and the volume-specific surface area for both volumes. It was shown that the porosity (36.3 and 38.0%, respectively) and geometric tortuosity values (2.0 and 1.9, respectively) agreed well between the two datasets despite the significant difference in image resolution. However, the volume-specific surface area of the nano CT scan was more than twice as high compared to that of the micro CT scan (1.3 μm 2 μm −3 and 0.6 μm 2 μm −3 , respectively). This finding is comparable to the analogy by Mandelbrot [156] who stated that the length of a coastline depends on the resolution of the map.

Summary
Moreover, the size of the sample volume has to be sufficiently large so that extracted values are representative of the sample bulk [157]. Hence, the higher the resolution, and the larger the extracted volume, the more likely the extracted values are accurate, and representative, and not affected by microscopic heterogeneities. Hence, the selection of resolution and the sample volume are crucial to ensure high accuracy in calculated microstructural data and an appropriate imaging technique has to be employed. At the same time, these factors become limiting when analysing structures where different length scales are important such as in battery materials [158]. When comparing the work by Wilson et al. [102] and Iwai et al. [103], who analysed the same type of sample using the same imaging technique achieving similar pixel sizes, no difference in tortuosities is observed even though Iwai et al. reconstructed a nine times larger sample volume. Similar findings are revealed when comparing Laurencin et al. [150] and Tjaden et al. [85]. Yet, this does not contradict the previous statement as homogeneous samples will yield representative values even for small sample volumes. Also, analysing the same sample type does not necessarily mean that these are structurally similar.
In addition, the above comparison revealed that different flux-based tortuosity calculation algorithms yield comparable results, which was also affirmed when executing different algorithms on the exact same sample [103,137]. This suggests that the choice of a flux-based computation algorithm has a smaller effect on the results than sample preparation technique, imaging parameters and the structure of the sample itself, which includes pore size distribution and volume fractions of the constituent phases. The interplay between these additional parameters and the tortuosity is visible when inspecting the work by Wilson et al. [102] and Izzo et al. [81]: while Izzo et al. presented higher tortuosity values, the porosity of their sample is a factor of 1.5 higher. Hence, the tortuosity itself does not give a full picture of the microstructure and the performance of the analysed sample, but has to be evaluated with respect to other microstructural characteristics [115,159]. Also, care must be taken when applying purely continuumbased models which do not account for Knudsen diffusion effects, such as the heat flux simulation. Such simplifications might cause visible differences between experimental-and simulation-based results [85].

Concluding remarks
The large number of tortuosity calculation methods is testimony of the significance of tortuosity in the field of electrochemistry. Here, we have reviewed different tortuosity calculation approaches which span from porosity-tortuosity correlations and image-based techniques to experimental methods.
Among these, a certain trend is revealed: porositytortuosity relationships, such as the Bruggeman equation, are more common in battery and PEM research, while flux-based algorithms are popular in SOFC research. Yet, each approach features distinct advantages and disadvantages. While easily applied, the Bruggeman relationship is valid only for spherical structures and is generally unfit to predict accurate values for complex porous networks. Similar caution must be exercised when applying image-based tortuosity calculation algorithms, and one must be aware of the difference and significance of geometric-and flux-based tortuosity. Results of either calculation procedure differ visibly, where geometric values lie below flux-based algorithms.
Moreover, tortuosity values determined using experimental techniques are valid only for the specific experiment at hand, as changes in temperature, set-up and gas composition affect the results. Also, tortuosity is usually used as a fitting parameter in these cases. The results are thus highly dependent on the applied calculation model. Furthermore, when comparing fluxbased algorithms across similar sample types, it is shown that tortuosity is a complex function of microstructural parameters and has to be interpreted while taking pore size distribution and volume fraction of constituents into consideration. Also, purely continuum-based models, which do not consider Knudsen effects, have also to be applied with caution as they might be the reason for discrepancies between simulation and experimental values.

Disclosure statement
No potential conflict of interest was reported by the authors.