Review of parameterisation and a novel database (LiionDB) for continuum Li-ion battery models

The layer thicknesses of components within the cell, as depicted in figure 1, are straightforward to understand and relatively easy to measure, e.g. using calipers or a micrometre gauge. Electrode layer thicknesses are closely linked to the areal capacities of active material coatings. Reported electrode thicknesses either measured or used in DFN models typically fall in the 10–100 µm range, and commercial polymer separator thicknesses in the 15–30 µm range. Layer thicknesses may be affected by the fabrication process, so these properties should be understood as their values after cell formation. This distinction is particularly important for the separator, which is usually made of a porous polymer, whose thickness prior to cell construction can be substantially different from the thickness after being wetted with electrolyte and wound or stacked in a cell, incurring a static mechanical load.

It is also possible that layer thicknesses can change with cycling of a battery, a process that requires further characterization and modelling. These temporal variations are almost always neglected in DFN models, where thicknesses for the electrodes and separator are usually taken to be those measured individually, prior to assembly [17], or following teardown after use [19]. Volume expansion or contraction of the active materials accompanies the lithiation/delithiation processes [53, 54], which is causes concomitant thickening/thinning of the electrodes. For silicon-based anodes, the volume expansion associated with lithiation is so large (up to ∼400% for pure Si) that it cannot be justifiably neglected. In that case some modelling work has aimed to extend the DFN framework to couple volumetric changes with cell electrochemistry, see for example McDowell et al [55] and Zhang et al [56].

The two electrodes usually contain four distinct phases: active-material particles, a conductive additive, a polymeric binder, and the pore-filling phase. The conductive additive, most commonly carbon black, is present to increase the electronic conductivity of the composite electrode, aiding electron transfer between the particles and current collector. It is often extremely fine (much finer than the electrode particles) and is typically well mixed with polymeric binder, such that the two phases are hard to distinguish using most imaging techniques. For the purposes of DFN modelling these two phases are usually lumped together and referred to collectively as the carbon–binder domain (CBD). Note that incomplete wetting by electrolyte during cell formation or side reactions during cycling can lead to the formation of a fifth, gaseous pore-filling phase [57], which is typically neglected in DFN models for Li-ion cells, although it has been considered in some cases [58].

Many models, including the original DFN models, set the volume occupied by the CBD to zero [5, 6, 19]. This assumes that cells are manufactured to maximise energy density by minimising the required CBD volume. As such, the electrolyte and active material volume fractions are related via:

$$\begin{align} \varepsilon_k + \varepsilon_{\textrm{act},k} = 1. \end{align} \tag{ 29 }$$

Here, ε_k is the electrolyte volume fraction (or solid-phase porosity) and $\varepsilon_{\textrm{act},k}$ is the active-material volume fraction. Figure 3 illustrates the distribution of porosities for positive and negative electrode materials, along with polymer separators sampled from the LiionDB parameter database. Porous electrodes in LiBs have electrolyte volume fractions of approximately 0.4–0.5, while separators skew towards higher porosities.

Zoom In Zoom Out Reset image size

Neglecting the CBD volume fraction via the assumption of spherically symmetric electrode particles ignores the fact that parts of an electrode particle's surface will be in contact with the CBD, instead of the electrolyte, resulting in some portion of the particle's surface not being available for (de)intercalation. This sort of consideration could dramatically impact model assumptions, particularly the spherically symmetric diffusion approximation in the microscopic equations and boundary conditions.

The only role that the CBD plays in the model equations presented here is indirect: the CBD volume fraction determines the effective electronic conductivity of the electrodes, $\hat{\sigma}_k$, a quantity that is typically measured directly using electrochemical impedance spectroscopy (EIS). The active-material volume fractions, in combination with the layer thicknesses, are needed to determine the theoretical areal capacity of the cell. They are also used to calculate the maximum Li concentration in active materials, as well as for anode–cathode stoichiometric balancing, as discussed in section 4.1.2. The electrolyte volume fraction has a significant role in determining the rate performance of the cell, because it impacts the transport of ions in the electrolyte required for charge transfer. Further discussion of porosity and its effect on transport efficiency is provided in section 3.1.5. On a practical note, incomplete electrolyte infiltration within pores is known to reduce cell performance. Gas generation from degradation during operation may also cause additional inactive regions [57, 59].

The surface area of active material that is available for intercalation per unit volume of electrode, b_k , depends on the size and shape of the active particles as well as their arrangement with respect to the CBD [35]. It comes into the liquid-phase charge balance in (3) when relating the local interfacial current densities on pore surfaces to the current transferred per unit of porous-electrode volume. As such, the surface-area to volume ratio impacts the rate at which lithium can transfer between the liquid electrolyte and the solid active materials.

Naively, one might expect that b_k should be maximised for the optimal performance (to provide lots of surface area for intercalation); however, in reality, the active material surfaces tend to form a layer called the solid electrolyte interphase (SEI), which consumes lithium inventory, thereby decreasing device capacity, and so a balance must be struck. In addition, the practicalities of manufacturing with very fine powders, which are hard to produce and process, can also be a consideration.

As highlighted in section 3.2, and plotted in figure 4, the measured surface area varies dramatically depending on the measurement technique used. There is also little consensus on the choice of surface area, with literature measurements choosing between definitions such as: (a) the whole solid-pore interface (i.e. including CBD surface), (b) the active particle surface area, and (c) the active particle-pore interface. As stated in section 4.1.3, the surface area defined must be paired carefully with measured exchange current densities to parameterise electrode kinetics consistently.

Zoom In Zoom Out Reset image size

The radii of the active-material particles establish the distance that intercalated lithium need to diffuse to reach the core, as well as determining the characteristic relaxation time of solid-phase diffusion, which in turn controls rate performance. Large particles do not perform well at high rates, because the depletion of or saturation with lithium at particle surfaces limits the accessible range of reaction current, thereby inhibiting performance by lowering power efficiency.

Real electrodes contain a distribution of particle sizes and shapes that confounds the determination of a single representative radius for each electrode, as is required by the basic DFN model [60]. For polycrystalline materials, the choice of a characteristic length based on primary or secondary agglomerated structures can change the radius parameter by a factor of 20 [61]. Particle radii tend to fall below 20 µm according to figure 5. The spread in the literature owes to the choice of primary or secondary structures, along with the variation among active materials.

Zoom In Zoom Out Reset image size

The particle radius plays a large role in the way intrinsic parameters of an electrode material are extracted, as discussed in section 4. While average radii are adequate for most models, newly engineered electrodes with graded radii may require models that treat particle radii as a variable with position [35, 62]. Different approaches towards averaging can also alter the effective radii [35]. In addition, some materials, such as graphite platelets, tend to form highly non-spherical particles [40, 63]. Anisotropic lithium diffusivities also result in different observed behaviours than that expected for an isotropic sphere (particularly in the case of graphite) [64]. For these reasons, defining a particle radius to characterise and model electrode geometries is not straightforward, even if all of the true morphological information is known.

The original DFN formulation captures the impact of pore-network geometry on electrolyte transport properties by using effective transport parameters $D_{\textrm{eff},k}$ and $\sigma_{\textrm{eff},k}$. When a diffusion medium can be suitably homogenized, the domain-specific properties of a pore-filling electrolyte relate to their values in an unobstructed electrolyte through a domain-specific transport efficiency, $\mathcal{B}_k$, as mentioned above. Figure 6 shows that transport through the porous layers in LiBs is typically constricted by more than 60%.

Zoom In Zoom Out Reset image size

It is also common to see MacMullin numbers ($N_{\textrm{M},k}$) being used to quantify the hindrance of transport caused by the electrode geometry [65]. As $\mathcal{B}_k$—the inverse of $N_{\textrm{M},k}$—is bounded between 0 and 1, it is comparatively more convenient for manipulation and interpretation. Transport efficiency is often further decomposed into two distinct contributions: the reduction due to the pore volume fraction, captured by ε_k , and the reduction due to the pore tortuosity factor, $\unicode{x03C4}_k$. The four macroscopic factors that characterize pore networks relate through

$$\begin{align} \mathcal{B}_k = \frac{1}{N_{\textrm{M}, k}} = \frac{\varepsilon_k}{\unicode{x03C4}_k}. \end{align} \tag{ 30 }$$

Measurement and interpretation of tortuosity factors has been the subject of extensive discussion [66].

Fuller et al [6] make use of the well-established Bruggeman correlation (i.e. $\unicode{x03C4}_k = \varepsilon_k^{1-b}$, where b is 1.5 for effective transport through spherical packed beds) to approximate the tortuosity factor. We note that there are some subtleties in how Doyle et al [5] implement the correction; for example, they take a factor of ε_k out from $D_{\textrm{eff},k}$, since their variable $N_\mathrm{e}$ describes anion flux through a whole volume element, whereas their $c_\mathrm{e}$ is the concentration of anions in the pores.

Bruggeman's tortuosity correlation is commonly deployed in situations where experimental geometrical information is limited. Although it is a convenient means to approximate the effect of microstructure on transport, it has been shown that the assumptions in its derivation do not necessarily align with the conditions established by the microstructures of real battery electrodes [67]. For example, use of a Bruggeman factor suggests that all electrodes with the same average particle radii and porosity would have the same pore-phase transport properties, but this is clearly not the case. Although correlations can be extended to account for more complex scenarios [68], many additional simplifications and assumptions are still required.

Unlike porosity, which has a meaningful definition within a planar differential cross section (ratio of cut-through pore area to total plane area), the tortuosity factor (and therefore ${\cal{B}}_k$ and $N_{\textrm{M},k}$) does not. This is because the tortuosity factor is an emergent property of transport through a microstructure; therefore sufficient path length (i.e. representative volume) is required before a value can be meaningfully measured. The multiscale approach used to construct the DFN model assumes that the characteristic lengthscale of the electrode particles is negligible compared to that of the electrode thickness, but this is not the case, with typical scale differences of only a single order of magnitude. The error that this incurs is likely to be related to the ratio of particle size to electrode thickness, i.e. of the order of 10%.

More recently, alternative concepts, including the 'electrode tortuosity factor' [69] have been reported that are able to measure the effect of graded microstructures. Many literature examples report 'geometrical tortuosities' or 'path-length tortuosities', calculated using a wide variety of algorithms applied to image data. This is distinct from the tortuosity factor which is found from using the multiscale approach. However, as made clear by Epstein [70], although path-length tortuosities are relevant to systems with well-defined paths, e.g. rivers and veins, they are not relevant to the types of microstructure found in real battery electrodes, unless they happen to consist of 'capillary' pores with relatively constant cross section.

In DFN models, the effective electronic conductivity of the electrodes is often assumed to be constant within each electrode, taking values of $\hat{\sigma}_\mathrm{n}$ and $\hat{\sigma}_\mathrm{p}$) in the negative and positive sides respectively. Probe methods capture the effective conductivities directly, and are therefore usually adopted in practice, although two-electrode measurements of direct-current resistance are sometimes used [71, 72]. Although the active particles are electronically conductive, they have a much lower conductivity than the CBD, and so it is often relatively accurate to approximate the conductivity of the electrode as a whole using only a contribution from the porous CBD. In practice, electronically-isolated regions of active material can become disconnected from the CBD network [73]. In reality such particles are then inert, and cannot participate in battery operation, but these local effects are not captured by the effective quantities.

Scanning electron microscopy (SEM) and transmission electron microscopy (TEM) have been used extensively for 2D image acquisition. SEM scans a beam of electrons across the target in a raster pattern and then detects the scattered electrons, with some devices achieving a resolution of below 1 nm [74]. By contrast, TEM transmits electrons through the sample before being detected, offering a way of capturing the internal detail, such as crystal structure and morphology. Additionally, TEM can achieve finer resolution than SEM, with some reports of 50 pm. However, samples need to be less than ∼10⁻⁷ m thick (much thinner than practical electrode coatings) to allow full electron penetration.

2D images can then be processed to extract measurements for layer thicknesses, particle radii and particle size distributions. It is important to note that because sample sizes are often very small compared to the entire electrode, a given image may not be statistically representative of the whole structure. It is generally difficult to obtain sufficiently large datasets that are truly representative for analysis.

Both focused ion beam-scanning electron microscopy (FIB-SEM) and x-ray computed tomography (XCT) have been used extensively in 3D studies; many examples can be found in fuel cell [75, 76], and battery research [77–79]. The FIB-SEM method images the 2D surface using an SEM, with the same resolution capability as SEM, and then mills the surface away using a controlled beam of charged positive particles, typically gallium, to reveal another layer. These 2D images are then reconstructed to produce a 3D microstructural volume. The milling process means FIB-SEM is destructive, i.e. in situ or operando imaging is not possible [80]. XCT systems fire an x-ray source at the sample and measure the resulting attenuation pattern, producing a 2D projection image. The sample is rotated to obtain images from many angles, which are then reconstructed, typically through a back-projection algorithm, to produce a 3D volume [81]. Image analysis methods are used to directly characterise the 3D microstructure, helping to understand particle distribution and anisotropy [82–84].

Image-based simulations can calculate geometric parameters, such as effective conductivities and volume fraction. Several open-source tools and methods for data processing are available [85]. Two popular options are (i) TauFactor, an open-source MATLAB (MathWorks) application for efficiently calculating the tortuosity factor, volume fraction, surface area and triple phase boundary density from image-based microstructural data [86] and (ii) OpenPNM, an open-source software package developed in Python for performing physics-based pore network simulations [87]. Other image-based mesostructure simulations also account for the CBD in electrode models, to elucidate effective transport properties and extensions into mechanical stresses [88, 89].

EIS is an electrical technique useful to probe various phenomena within LiBs [90]. This section describes the technique but focuses on geometric parameters; sections 4.2.6 and 5.2.1 will address material and interfacial properties.

The EIS method consists of applying a small sinusoidal voltage or current perturbation and recording the respective current or voltage response, whose resulting amplitude and phase shift quantify the impedance; an impedance spectrum is then recorded as a function of the frequency of the excitation signal. Although frequencies ranging from µHz to MHz are in principle accessible, most information is gathered in the range between 0.1 Hz and 100 kHz, where signal-to-noise ratios are good and reproducibility is easier to achieve. EIS provides multiscale resolution, because different physicochemical processes have distinct characteristic timescales, excited at different frequencies within the sampled spectrum.

Diffusive transport—both in electrode active particles and the electrolyte—has characteristic frequencies in the range below 1 Hz. Interfacial effects, such as charge-transfer resistance (kinetic rate constants), double-layer capacitance, and possibly effects arising from surface passivation are resolved at intermediate frequencies, between 1 Hz and 1 kHz. At higher frequencies, above 1 kHz), ohmic and ionic resistance and inductance can be observed. Notably, the linkage between impedance data and material or geometric parameters requires an appropriate model be fitted to the observed spectrum. Often the choice is an ECM whose parameters relate heuristically to physical parameters [91, 92], but formally linearised electrochemical models may also be used [16, 93].

Generally, division of the effective electrolyte conductivity extracted from EIS by the ionic conductivity independently measured for unobstructed liquid directly yields the transport efficiency $\mathcal{B}_k$ discussed in section 3.1.5. Landesfeind et al [65] surveyed many specific details about methods by which Macmullin numbers or transport efficiencies can be determined. Using fits of EIS data with transmission-line models, they determined transport efficiencies in a range of porous electrodes and compared these with values calculated from the Bruggemann correlation typically used in DFN models.

Laser diffraction particle sizers that can measure the distribution of particle radii within powders are commercially available. Typically they use photodetectors to measure the angle and intensity of light scattered through a sample. The collected pattern is then transformed into a particle-size distribution via the Mie theory of scattering. It is common to describe the spread in size distribution with D₁₀, D₅₀ and D₉₀ diameters, representing the size at the 10th, 50th and 90th percentiles, respectively. This technique has been widely applied to active-material samples [60, 94, 95]. A flaw of particle-size measurement via laser diffraction is that it assumes scattering from perfect spheres and generally does not provide accurate distributions from the irregularly shaped particles that make up some battery materials [96].

Mercury intrusion porosimetry is a destructive technique used to measure porosity, volume of pores, pore size distribution and surface area. In a measurement, increasing pressure is applied to force mercury into smaller and smaller void spaces in the porous medium being studied, and recording the volume of intruded mercury [97].

Smaller pores may not always be fully accessible to mercury intrusion porosimetry, owing to the higher pressures necessary [98]. As such, the method is only suitable for battery electrode coatings with larger pore sizes [17, 19, 99]. It can also be used to measure particle radii within electrodes by exploring the contact angle, hydrostatic pressure, and surface tension of the mercury [99]. Similar to laser diffraction, the principles of mercury intrusion porosimetry derive from the assumption of perfectly packed spheres [100]. Geometrical parameters for inhomogenous active material particles may be more precisely extracted using imaging methods.

The Brunauer–Emmett–Teller (BET) adsorption isotherm can be used to measure surface areas within porous media by flushing the pore structure with nonadsorptive gas such as helium, then inducing the adsorption of a surface monolayer of an adsorptive gas, such as nitrogen [101]. Based on formulas derived from the BET theory, the quantity of gas adsorbed correlates with the total surface area of the particles. By tracking the hysteresis of nitrogen adsoprtion and desoprtion, experimental equipment that performs measurements with the BET method can often also be used to infer pore-size distributions.

BET adsorption is a prevalent tool, widely used to extract solid electrode-material surface areas under varying synthesis conditions and degradation/cycling protocols [102–105]. Shortcomings of the technique include the high sensitivity to sample preparation, degassing and lengthy measurement time. BET adsorption cannot deliver very accurate parameters for DFN models, because the technique can only be used on dry powders, and therefore does not characterise the actual areas in chemical contact or available for electron exchange within porous electrodes.

The OCV of a full cell is defined as the potential difference between the positive and negative electrodes when there is no current flow and the cell has relaxed to equilibrium. The OCV is not a constant, but is instead a function of the thermodynamic states of the materials in the electrodes. Generally, the OCV varies significantly with respect to its degree of lithiation (i.e. its SoC), and also slightly with temperature.

To characterise the full-cell OCV, the variation of potential of each of the active materials at equilibrium is measured, both relative to the same reference-electrode reaction [107]—most commonly lithium redox on lithium metal—taking the difference between these two voltages (often termed half-cell OCVs) gives the full-cell OCV. This data can then be supplied to the model either in the form of a look-up table or a suitable mathematical function fitted to the raw data. OCV data is commonly fitted to either (i) functions that are empirical but convenient to work with (e.g. piecewise linear functions [108]), see e.g. [109], or (ii) physically motivated functions [106, 110]. Some of the challenges associated with data-fitting in practice have been discussed by Plett and colleagues [111, 112]. The original modelling work by Newman and co-workers used either a form of the Nernst equation [5] or an empirical expression including hyperbolic tangent terms, exponentials and a constant [6]. Similar approaches have been widely adopted in many other more recent works [44, 113, 114]. Some authors have also fitted a polynomial to the measured data, which may be sufficient for many applications and is simpler to implement [115].

Zhang et al [110], Birkl et al [106] and Verbrugge et al [116] all demonstrated that a realistic function for each half-cell OCV can be derived from first principles using a substitutional lattice gas model. The derivation of these functions is similar to that of a Nernst equation, but is extended to include the energetics of short-range interactions and local equilibria among sites with different energies. The model was demonstrated to accurately represent a wide range of common materials, including graphite, layered nickel-manganese-cobalt oxide, manganese oxide, and iron phosphate. The function is a summation of sigmoidal terms for each plateau in the OCV. Zhang et al [110] and Birkl et al [106] independently derived an equation for the fraction of sites within the active material occupied by lithium, θ_k , in terms of the OCV, U_k :

$$\begin{align} \theta_k(U_k) = \sum_{i = 1}^N \frac{\Delta\theta_i}{1 + \exp{\left[ \left(U_k - V_{0,i}\right) a_i e / kT\right]}}. \end{align} \tag{ 31 }$$

The voltage of each plateau is $V_{0,i}$; the fraction of the SoC range in which the material resides on the ith plateau is $\Delta \theta_i$, and a_i is a fitting parameter that approximates the interaction energy between neighbouring ions in the lattice. The expressions used by Zhang et al [110] and Verbrugge et al [116] are similar to (31), but include a lattice-energy parameter equivalent to the reciprocal of Birkl et al's a_i . Being analytic, (31) is smooth and is therefore apt for use in numerical computations. Unlike the familiar Nernst equation, however, the composition is given explicitly, rather than voltage being an explicit function of the composition. Although expression (31) is generally invertible because θ_k always decreases monotonically with U_k , the inverse function $U_k(\theta_k)$ cannot generally be found analytically, posing a minor implementational inconvenience.

Hysteresis or path dependence effects are usually observed in OCVs, manifesting as a different OCV function being measured during charge, compared to discharge. The reasons for this are widely debated, but generally accepted to be complex and linked to the formation of metastable states within active materials. Hysteresis may be examined using isothermal calorimetry [117], and has been modelled empirically by constructing a differential equation that describes how the hysteresis voltage changes with respect to SoC and history [118].

The maximum lithium concentration, $c^{\max}_{k}$, dictates the extent to which a solid electrode material can incorporate lithium. Materials that can store large amounts of Li per unit weight or volume are generally sought after, since they increase the energy density of the battery. The parameter $c^{\max}_{k}$ is material-dependent and typically lies between 20 000–55 000 mol m⁻³ for metal-oxide cathodes, and 15 000–35 000 mol m⁻³ for graphitic anodes. This parameter appears in model equation (21), defining the Butler-Volmer reaction rates. Notably, for intercalation compounds that form solid solutions, the exchange current density, $j_{k0}$, vanishes either when the lithium concentration at the electrode surface approaches its maximum or when it reaches zero.

The theoretical capacity and $c_k^{\textrm{max}}$ for material k may be calculated, through Faraday's law, by using the active material crystal density, $\rho_k^\textrm{crystal}$ and the molecular mass of the lithiated active material, m_k . Chen et al [19] estimated similar values of $c_k^{\textrm{max}}$ using both the material coating mass per unit area $M_{\textrm{coat},k}$ and active material crystal density, through

$$\begin{align} c_k^{\textrm{max}} = \frac{\rho^\textrm{crystal}}{m_k} \approx \frac{M_{\textrm{coat},k}}{m_kL_k\varepsilon_{\textrm{act},k}}. \end{align} \tag{ 32 }$$

Individual electrodes can typically access maximum and minimum lithium concentrations that lie beyond the capacity ranges accessed in a full cell, as depicted in figure 7. Most positive electrode materials, such as LiCoO₂ and nickel-manganese-cobalt oxide, LiNi_x Mn_y Co_z O₂ (NMC), undergo irreversible changes of their crystal structure at low stoichiometry, which should be avoided to minimise degradation [119]. Full lithiation of the positive electrode is also not usually possible because cells are fabricated with fully lithiated cathode materials; some capacity is inevitably lost due to the loss of lithium inventory during initial SEI formation. For the negative electrode, it is important to leave some margin to avoid overcharging or overdischarging a cell.

Zoom In Zoom Out Reset image size

The literature is often unclear whether stoichiometry is defined in terms of (i) a fractional of the cell's SoC (ranging between 0 and 1) that is accessible in practice, but that does not necessarily correspond to concentrations ranging between 0 and $c_k^{max}$, or (ii) the true theoretical maximum. Adopting a perspective in which lithiation limits are imposed, the fractional SoC, $\Theta(t)$, can be defined as

$$\begin{align} \Theta(t) = \Theta(0)-\frac{A\int_0^t i_\mathrm{circuit}(t^{\prime})\,\mathrm{d}t^{\prime}}{Q_\mathrm{cell}}, \end{align} \tag{ 33 }$$

where $Q_\mathrm{cell}$ is the capacity of the cell when fully charged (but the anode material may not necessarily be fully lithiated and, likewise, the cathode material may necessarily not be fully delithiated), so that $\Theta\gt1$ corresponds to overcharge and $\Theta\lt0$ corresponds to overdischarge. Each electrode's state of lithiation can be expressed using a variable ($\theta_k(U_k)$ in equation (31)) which determines the fraction of crystal-lattice sites occupied by Li⁺. Since $\theta_k(U_k)$ varies in space as well as over time, it is useful to define a volume-averaged counterpart $\langle\theta_k\rangle$:

$$\begin{align} \langle\theta_k\rangle(t) = \frac{\int \theta_k(x,t)\,\mathrm{d}V}{V} = \frac{\int c_{k\mathrm{s}}(x,t)\,\mathrm{d}x}{L_kc_{k,\textrm{max}}}. \end{align} \tag{ 34 }$$

Note that $c_{k\mathrm{s}}$, not c_k , is included in (34), because θ_k is a function of the half-cell OCV U_k , which is a function of $c_{k\mathrm{s}}$ and not c_k . The full-cell SoC $\Theta(t)$ and the half-cell lithiation fractions, $\langle\theta_\mathrm{n}\rangle(t)$ and $\langle\theta_\mathrm{p}\rangle(t)$, are linked by

$$\begin{align} \Theta(t) = \Theta(0) + \frac{FA}{Q_\textrm{cell}} \left\{L_\mathrm{p} c_{\mathrm{p},\textrm{max}}[\langle\theta_\mathrm{p}\rangle(t)-\langle\theta_\mathrm{p}\rangle(0)] - L_\mathrm{n} c_{\mathrm{n},\textrm{max}}[\langle\theta_\mathrm{n}\rangle(t)-\langle\theta_\mathrm{n}\rangle(0)] \right\}. \end{align} \tag{ 35 }$$

Observe that the electrode-domain thicknesses $L_\mathrm{p}$ and $L_\mathrm{n}$ appear here. One important use of DFN models is to optimize electrode balancing by illustrating how choices of $L_\mathrm{p}$ and $L_\mathrm{n}$ affect performance.

Because the full range of individual electrode stoichiometries, $0\lt\theta_\mathrm{p} \lt1$ and $0\lt\theta_\mathrm{n} \lt1$, extends beyond the accessible full-cell SoC range, $0\lt\Theta\lt1$, full-cell OCV measurements are not sufficient to determine $\langle\theta_k\rangle$; in principle half-cell thermodynamic data from the literature or half-cell measurements with harvested electrodes are also necessary. To match half-cell OCV curves with full-cell OCV data, one can examine turning points in differential voltage measurements [17], perform least-squares regression [99], or implement three-electrode-cell measurements [19, 120]. The initial values of $\langle\theta_k\rangle$ relate to the mass balance between the anode and cathode, the first cycle loss observed during the formation process (usually taking place during manufacturing) and the state of health.

Figure 7 illustrates the range of OCV curves reported in the literature for graphitic anodes and NMC cathodes. There is a lack of consensus between reporting of cell-level SoC, $\Theta(t)$, and materials-level lithiation scales, $\langle\theta_k\rangle(t)$, shown particularly by the variation in lithium stoichiometry and measured half-cell OCV for NMC cathodes. Many papers are careful to distinguish between the full-cell state of charge Θ and electrode stoichiometry θ_k , but some [30, 121, 122] use the term SoC for both, which is confusing for readers trying to parameterise a battery model. figure 7 also highlights an example of electrode stoichiometric balancing, from the work of Birkl et al [106], showing how the relative positions of half-cell OCV curves and their minimum/maximum lithiation limits combine to form the full-cell voltage curve.

Reaction kinetics in porous electrodes are typically taken to follow the Butler–Volmer kinetic equations (20) and (21), which relate the interfacial current density j_k to the surface overpotential η_k . The exchange current density $j_{k0}$, and hence the reaction rate, is proportional to the Butler–Volmer kinetic rate constant, k_k . A sample distribution of kinetic rate constants for electrode materials reported from the LiionDB parameter database are shown in figure 8. Experimentally obtained exchange current densities vary by several orders of magnitude, with some variation owing to how they are defined [123, 124]. As highlighted by Dickinson and Wain, this variability may stem from differences in how variables in the Butler–Volmer equation are interpreted by electroanalytical chemists and electrochemical engineers [125]. Differences also arise due to parameter values assumed when processing experimental data, such as effective electrode surface area, discussed in section 3. Caution is advised when comparing published values of $j_{k0}$ and k_k , since the Faraday constant, F, may be omitted from the definition of the interfacial or exchange current densities in some papers, and the effective surface area may not be taken into consideration [5, 44, 126].

Zoom In Zoom Out Reset image size

Ecker et al [17] found that equation (20) was a good fit to exchange current densities measured by EIS for LiNi_0.4Co_0.6O₂ (NMC406), but not for graphite. Schmalstieg et al [99] also found that (20) was a good fit to exchange current densities measured by EIS measurements for NMC811, but the exchange current density had virtually no SoC dependence for graphite. These findings are consistent with SEI-layer properties interfering with the measurements on negative electrodes. Conversely, the same paper by Ecker et al [17] found that the Butler–Volmer equation (8) provided a good fit for EIS measurements of charge-transfer resistance on graphite, but not for NMC406. Similarly, Ko et al [61] found that $k_\mathrm{p}$ extracted from galvanostatic intermittent titration technique (GITT) measurements on NMC111 is not constant, but instead varies with SoC. Recent literature has also suggested that alternative kinetic mechanisms better describe electrochemical reactions at high rates and overpotentials, or with advanced battery formats, such as lithium metal and anode-free cells [123, 124, 127, 128]. A rigorous justification for the common assumption that $\alpha_{\textrm{a} k} = \alpha_{\textrm{c} k} = 0.5$ is also lacking [4, 129], but the factors $\alpha_{\textrm{a} k}$ and $\alpha_{\textrm{c} k}$ are notoriously difficult to measure, though recent developments in nonlinear EIS has shown sensitivity to these reaction symmetry parameters [130]. Deviations in these factors from 0.5 has been used to explain the fact that LiB discharging is faster than charging at a given magnitude of cell voltage [131].

Thermal activation of intercalation kinetics in materials such as graphite and NMC is widely accepted, with Arrhenius relationships assumed between exchange current density and temperature having been reported throughout the literature [17, 19, 99, 132]. Despite being called into question from both theoretical and experimental perspectives, the theory of thermally activated Butler–Volmer kinetics is almost always used in LiB models, for two reasons. First, it requires only two parameters: the rate constant in a reference state, k_k , and associated activation energy. Second, sensitivity analyses [133–135] have shown that DFN models are more sensitive to particle radii, r_k , and diffusivities, D_k , than to reaction kinetics.

In the active materials, the Li⁺ diffusivity, D_k , controls the transport of Li⁺ within the particles. It becomes very important at high currents, where concentration gradients within the active material may become the main source of overpotential if diffusivities are not large enough. When diffusivities are too small, particle surfaces easily become depleted of or saturated with lithium, which shuts down reaction kinetics and thereby strongly impacts battery operation [133–135].

The first DFN models assumed that solid phase diffusion follows Fick's second law [8], i.e. they used (17) and took D_k to be constant. Studies on parameter identifiability have shown that diffusivity is intrinsically linked to the particle radii, $r_k^2 / D_k$ [15, 16]. It should be emphasised that many of the experimental techniques, such as GITT, infer diffusivities under the assumption that a spherical form of Fick's second law applies and so this model assumption would be required for consistency.

Despite the constancy of D_k assumed in the original DFN models, measurements indicate that $D_\textrm{k}(c_k)$ can vary by orders of magnitude, as the Li⁺ concentration changes [17, 19, 99, 121, 136–139]. The data plotted for positive and negative electrodes in figure 9 show that Li⁺ typically diffuses through the active materials in the range of 10⁻⁸ to 10⁻¹⁴ cm² s⁻¹. For graphite negative electrodes, $D_\mathrm{n}(c_\mathrm{n})$ is large for concentrations corresponding to material phase transitions between plateaus on the half cell OCV curve and is small for other stoichiometries. Active materials are typically polycrystalline and both modelling [140] and experimental [141] studies on NMC-type materials agree that diffusion is faster along grain boundaries than within crystals. Cracks form much more quickly during cycling in polycrystalline materials than in single crystal materials, however, causing a further decrease in diffusivity [141, 142]. Nanoscale LFP shows very fast Li-ion transport, but larger LFP crystals have much slower transport, owing to the transport along the rapid b channels becoming blocked by dislocations. Thus, active material quality has a significant impact on battery behaviour.

Zoom In Zoom Out Reset image size

Concentration-dependent diffusivities have been used in models, which have been shown to match experimental data better at various C-rates and temperatures [17]. They have also been used to investigate how phase changes affect differential voltage signals with and without lithium plating [129]. The temperature dependence of diffusion is well reproduced by an Arrhenius response [17, 139]. While Fickian diffusion is by far the most widely used model of mass transport within particles, there is increasing evidence that other models that account for phase change, such as the Cahn-Hilliard equations [143, 144], may be more appropriate for graphite [145] and LFP [45], although parameters such as the characteristic transition-region length do not have clear meanings in the context of classical thermodynamics.

Four-point surface probes are commonly used to measure the electronic conductivity of porous electrodes. A current is applied across the sample via two of the probes, while the other two measure the potential drop. The probe spacing is fixed, but varies from device to device: it can range from a few micrometers to a few millimeters [146]. By relating the applied current and the measured voltage, the electronic conductivity can be estimated [147, 148], using suitable corrections based on sample size and shape [149]. Current collectors can also be separated from the electrode coating to reduce their influence on measurements by strategies such as the use of strong adhesives to delaminate the electrode material or solvents to dissolve the current collector foil [19, 71, 150]. Due to the anisotropy of the electrode coatings, conductivity depends on the orientation of the sample [151]. In an operating battery, the through-plane conductivity (in the direction perpendicular to the current collectors) is most relevant for use in the DFN model, as this relates to the movement of electrons from the surface of the electrode particles to the current collector [152]. Bespoke methods have also been developed to measure the combined interfacial resistance [153, 154].

Through-plane conductivity can be measured using a two-electrode method [155, 156]. This involves sandwiching the sample between two probes, then evaluating the direct-current resistance. Deconvolution of contact resistances can be a challenge in the two electrode set-up; bespoke four-point probe designs offer greater accuracy [157].

During cycling of an active material, the overall voltage response is comprised of the equilibrium OCV of each electrode plus the overpotentials between each electrode and the electrolyte, and the resistive potential difference across the electrolyte and solid components of the electrodes. Each of these contributions to the potentials, except the OCVs, generally increase with current throughput. The galvanostatic pseudo-OCV method minimises the non-OCV contributions by sweeping through a cell's stoichiometric range at a slow rate (e.g. less than C/20—i.e. a current at which the rated capacity of the cell would charge or discharge in 20 h) to determine the cell OCV. It should be noted that this method never truly reaches an open-circuit condition, but can be made to come arbitrarily close to it by decreasing the C-rate. Chen et al [19] compared the pseudo-OCV approach to the GITT technique described in section 4.2.3, concluding that polarisation, kinetically stable phase transitions and hysteresis are still observed in most pseudo-OCV data.

Weppner and Huggins [158] introduced the GITT in 1977 to characterise diffusivities in dense planar electrodes. The method has since been adapted for porous electrodes. GITT involves applying a short and weak current pulse, followed by an open-circuit (no applied current) rest period which is sufficiently long that a cell voltage close to the equilibrium potential is obtained. The charge transferred should be small enough such that the voltage response remains proportional to the current applied, and data should only be processed outside the timescale dominated by capacitive relaxation of the electrolytic double layer. An exemplary set of GITT data is provided in figure 10(a).

Zoom In Zoom Out Reset image size

For intercalation electrodes, GITT is used to establish electrode OCVs at different states of lithiation [19, 60, 61, 159]. After each current pulse and subsequent relaxation step, the equilibrium OCV, U_k , is given (to a good approximation) by the cell voltage, and the lithiation concentration $c_{k0}$ is known from the pulse current and duration through Faraday's law of electrolysis. It should be emphasized that thermodynamic equilibrium is only well approximated by GITT-derived OCVs that have been recorded with sufficient rest periods after each pulse.

The quantification of the solid-state diffusivity from GITT data is based on an analysis with Fick's second law. Assuming linear diffusion limited by the solid electrode material, and that the pulse duration is short enough that diffusion within the electrode film can be assumed semi-infinite, the diffusivity D_k of lithium within the solid is given by

$$\begin{align} D_k = \frac{4}{\pi} \left( \frac{i_\mathrm{circuit} A}{S_k F} \right)^2~\left( \frac{U^{^{\prime}}_k(c_{k0})}{\frac{\partial V}{\partial\sqrt{t}}} \right)^2, \end{align} \tag{ 36 }$$

where A is the geometrical area of the porous electrode, S_k is the area of the electrode-electrolyte interface within it, $c_{k0}$ is the Li⁺ concentration in the solid at the beginning of the pulse, $U^{^{\prime}}_k(c_{k0})$ is the derivative of OCP with respect to concentration at $c_{k0}$, and $\tfrac{\partial V}{\partial\sqrt{t}}$ is the derivative of the terminal voltage with respect to the square root of time. Note that this analysis assumes a semi-infinite film penetration model of diffusion [160]; it is important to be aware of whether assumptions concerning the nature of diffusion are applicable to the system being studied [132].

The GITT method was originally derived for compact thin-film electrodes [158], but remains valid for porous electrodes, as long as S_k is taken to be the effective surface area available for charge transfer within the porous material and the configuration is such that spherical diffusion within the electrode particles does not become rate limiting. Generally, S_k relates to the DFN surface-to-volume parameter b_k , through

$$\begin{align} S_k = b_k A L_k , \end{align} \tag{ 37 }$$

where L_k is the thickness of the porous electrode. Assuming that the electrode comprises spherical particles, all of radius r_k , that occupy a volume fraction $\varepsilon_{\mathrm{act},k}$ of the active material, and that the distribution of particles throughout the film depth is uniform then, it follows that $b_k = 3~\varepsilon_{\mathrm{act},k}/ r_k$. Application of GITT to porous-electrode models is well detailed by Nickol et al [161].

A reaction-kinetic overpotential is present during any galvanostatic measurement, hence the GITT set-up should be chosen to minimise any overpotentials observed on the working electrode, maintaining accurate measurement of the thermodynamic states through the OCV. In two-electrode measurements care should be taken to ensure that these overpotentials are small enough that they do not appreciably contaminated the results. Even when using Li metal counter electrode half cells, the Li metal may contribute an additional overpotential, leading to inaccuracies in the observed cell voltages and underestimation of diffusion properties. Three- or four-electrode cells, incorporating one or two reference electrodes, produce more accurate voltage responses because they allow the contributions from the overpotential to be quantified and subsequently subtracted off of the data [161].

For sufficiently small currents and short time intervals, $\mathrm{d} V/\mathrm{d} \sqrt{t}$ is constant, experimental data can be used to estimate the change in OCV and the charge injected with the pulse directly relates to the concentration change in the solid through Faraday's law of electrolysis. Thus, one can use voltages from the data in figure 10(a) to simplify equation (36), yielding

$$\begin{align} D_k = \frac{4}{\pi \tau} \left( \frac{r_k}{3} \right)^2~\left( \frac{V_4 - V_0}{V_2 - V_1} \right), \end{align} \tag{ 38 }$$

where, apart from the measured voltages V_j , the only inputs are the particle radius, r_k , and the pulse duration, τ. This approach removes some error associated with overpotentials, since the equilibrium voltage is recorded without an internal resistance contribution. Poor estimates of diffusivity can still occur, however, particularly when GITT is undertaken in an SoC region with a relatively flat OCV response, or if equilibrium is not (approximately) reached by the end of each rest period. GITT measurements directly parameterise the diffusional time constant, $r_k^2 / D_k$, with the extracted diffusivity having an inherent dependence on the measured radii $r_k^2$. This link also makes it difficult to interpret GITT results for electrodes with a wide range of particle sizes. An alternative approach is to fit the voltage transient versus $\sqrt{t}$ during the current pulse, and calculate the diffusion coefficient directly using equation (36). This approach does not eliminate any overpotential errors. More elaborate methods to calculate S_k can be used, as described in section 3, but the value should always be consistent with that used in the DFN model.

Recently, there have been several developments in the application and analysis of this technique. A modified GITT technique, using a 'staircase' current profile, has been reported to determine other kinetic parameters, including the exchange current density and charge transfer coefficients, by fitting the observed overpotentials to a modified Butler-Volmer equation [162]. Lacey and colleagues have also modified the GITT technique to reduce experimental time and provided validation based on DFN model simulations [163, 164]. Other developments in the estimation of diffusion coefficients have focused on the methods to fit the model to the data [161, 165, 166].

Cyclic voltammetry (CV) is a technique commonly employed to investigate reduction and oxidation processes [160, 167], using a three electrode system—working, counter and reference electrodes—and measuring the current response to a linearly cycled potential sweep between two given voltage limits, as depicted in figure 10(d). In sweep voltammetry, an applied potential is varied linearly with time while recording the resultant current. This potential sweep is characterised by a particular scan rate, ν (with units V s⁻¹). Data is typically reported in the form of a cyclic voltammogram, which plots current versus voltage, both parametric functions of time. Analysis of cyclic voltammograms provides both qualitative and quantitative information about electrode reaction kinetics and mass transport. For electrode parameterisation, CV can be used to obtain kinetic rate constants of (de)lithiation [168, 169] and the diffusion coefficient of lithium within electrodes [170]. These techniques are discussed in detail in Kim et al [171].

The quantitative analysis of cyclic voltammograms varies according to whether the electrochemical reaction is classified as reversible (i.e. kinetically fast), quasi-reversible, or irreversible (i.e. kinetically sluggish). For a reversible reaction, the forward and reverse electron-transfer steps are in equilibrium and voltage follows a Nernst-like equation. Li-ion intercalation is often considered to be reversible [172], though quasi-reversible analyses have also been applied when rates of charge transfer and mass transport are similar [173]. Given a diffusion-controlled, reversible single-electron transfer process, the peak current is proportional to the square root of the scan rate ($I_\textrm{peak} \propto \sqrt{\nu}$).

For porous electrodes, the surface area and morphology of the surface also affects the shape of the voltammogram. As many common models that underpin CV analysis are derived based on assumptions that the working electrode is small compared to the counter and that liquid-phase diffusion is semi-infinite, translation of the theoretical equations to represent porous electrodes is difficult. The effects on peak voltage and reversibility are particularly pronounced on thicker porous electrodes [173]. The impacts of intercalation-particle geometry and dispersion on CV profiles have also been modelled [169]. It was found that accounting for the real particle-size distribution may be crucial for determining the voltammetric response. When mass-transfer contributions derive from both free electrolyte adjacent to the electrode and pore-filling electrolyte within it, reversible reactions may also produce an apparent shift in peak voltage position, proportional to the scan rate [173].

In general, thin electrode coatings minimise the observed kinetic limitations and the effect of the porous electrode geometry. Given a reversible CV response, the relationship between the peak current, $I_\textrm{peak}$, and the scan rate, ν, see figure 10(c), can be fit with the following equation to elucidate the solid-phase diffusivity of lithium, D_k [160, 173, 174]:

$$\begin{align} I_\textrm{peak} = 0.4463~F c_k S_k \sqrt{\frac{F \nu D_k}{RT}}. \end{align} \tag{ 39 }$$

Here c_k is the concentration on the surface of the particle and S_k is the effective surface area per unit mass. The effective surface area can be further refined to account for the electrochemically active surface area. Since practical battery electrode materials more often exhibit quasi-reversible or irreversible characteristics, CV parameterisation techniques may only offer qualitative comparisons between effective kinetic and transport properties.

Potentiostatic methods involve applying a voltage step and observing the current response. Examples include a voltage sweep, such as in CV, or constant voltage measurements, such as chronoamperometry. Chronoamperometry can also be applied over a range of concentrations or states of charge of the electroactive material and, at each concentration (or voltage), the current transient can be used to calculate the solid-state diffusion coefficient.

In the potentiostatic intermittent titration technique (PITT), a stepped constant voltage is applied and the current decay is tracked, at each voltage, to a set current limit. A rest period can also be introduced in between each voltage step, after the current decays such that OCV is (approximately) reached [175, 176]. The exclusion of rest periods accelerates data collection. However, including the rest period in PITT means the titration method can be used to calculate the diffusion coefficients while minimising overpotentials that contribute errors to the concentration estimation. Care must be taken with any overpotential observed in the measurement and, again, 3-electrode measurements with a reference electrode are significantly more accurate. A typical PITT data set is depicted schematically in figure 10(b).

A model based on Fickian diffusion is used to estimate diffusivity from current data [176]. In the case where $t \ll {r_k^2}/{D_k}$ (short times compared to the time scale of solid-state diffusion), and where interfacial reaction kinetics is reversible, it is appropriate to apply the Cottrell equation [160]:

$$\begin{align} D_k = \pi t \left(\frac{i(t)}{FS_k(c_{k\mathrm{s}}-c_{k0})}\right)^2, \end{align} \tag{ 40 }$$

whereas for $t \gg {r_k^2}/{D_k}$ (long times) the current decay becomes exponential and the solid-state Li diffusion coefficient can be determined from [176, 177]:

$$\begin{align} D_k = -\frac{4~r_k^2}{\pi^2} \frac{\partial \ln i}{\partial t} . \end{align} \tag{ 41 }$$

This has been expressed in terms of a particle radius, r_k , and, within the caveats described in section 4.2.3, it is true to say that $b_k = 3~\varepsilon_\mathrm{act,k}/ r_k$.

Malifarge et al [178] have provided guidelines for the implementation of PITT and its underlying assumptions. Less conventionally, PITT can also evaluate the exchange current density [179] but, as already mentioned (section 3.2.2), EIS is more often used. Diffusion coefficients via PITT and GITT have been shown to corroborate one another [180, 181]. The potentiostatic methods offer a speed advantage over galvanostatic methods, but additional measures may be needed to minimise the potential difference between the working electrode and the reference or counter electrode, and the voltage drop across the electrolyte. Not properly accounting for these overpotentials can lead to overestimates of the solid-phase diffusion coefficients.

The EIS technique was described in section 3.2.2 of this review and is revisited again in section 5.2.1. In the context of electrode parameterisation, EIS may be used to estimate charge-transfer resistance, double-layer capacitance and diffusion timescales, usually by fitting parameters of a simple circuit model to EIS data, by means described below and further detailed in [17, 33, 168]. A typical model is the Randles circuit which, per electrode, consists of a double-layer capacitance in parallel with the series combination of a Warburg element, representing diffusion, and a charge-transfer resistance, as shown in figure 10(c). Further series resistances may be added to represent additional ohmic and ionic voltage drops, as well as additional elements to represent electrolyte dynamics, if required.

4.2.6.1. Charge transfer

An expression for charge-transfer resistance may be derived by linearising the Butler-Volmer equation (20) for small values of the overpotential, i.e. about $\eta_k = 0$. On doing so, one finds that the charge transfer resistance (per electrode), $R_\textrm{ct}$, relates to exchange current density through:

$$\begin{align} j_{k0} = \frac{RT}{S_k FR_\textrm{ct}} = \frac{RT}{A L_k b_k FR_\textrm{ct}}. \end{align} \tag{ 42 }$$

Measurements of $R_\textrm{ct}$ by EIS have thus been used to determine kinetic parameters [93, 182].

Large variances among calculated exchange current densities $j_{k0}$ stem, in part, from the choice of surface area, with different groups using different definitions, as described by Chang et al [183] and Koch [184]. The exchange current density relates to the Butler–Volmer rate constant, k_k , through equation (21).

The identification of the $R_\textrm{ct}$ value associated with the main intercalation reaction for individual electrodes is a challenging task. Without a reference electrode, it is not possible to distinguish between the separate half-cell reactions (recall that in the model, there are in fact two Butler-Volmer equations, one for each electrode). Another issue is that the impedance response of the active materials and the response associated with the constantly evolving passivation surface films tend to occur in similar frequency ranges [185]. Changing porosity and surface roughness also leads to frequency dispersion [186]. The majority of strategies to fit $R_\textrm{ct}$ from overall spectra involve performing EIS measurements under conditions that cause the kinetic and surface-layer behaviour to diverge. This includes observing dependence on temperature [19], changes with conditioning cycles to develop SEI [187], or varying SoC [188]. Insight can also be gained by comparing EIS spectra obtained from full-cell and half-cell configurations [99]. Non-linear EIS has been implemented to analyse higher-order harmonics [130]. The temperature dependence of $R_\textrm{ct}$ has been modelled in terms of the activation energy for lithium transport through the SEI layer [189].

4.2.6.2. Diffusion

While charge transfer occurs at medium to high frequencies, the low-frequency EIS response of a porous electrode provides information about the solid-state diffusion coefficient in each electrode [190]. The considerations mentioned above also apply here, namely that, without a reference electrode, it is not possible to distinguish between the separate electrodes (note that this issue is not unique to EIS, it also applies to GITT and other techniques mentioned in this section). In the frequency domain, diffusion transport is characterised by a Warburg element, which may be derived for various different geometric configurations, including thin films, as well as planar, cylindrical, or spherical particles (see Song and Bazant [191]). For spherical particles, the diffusion-related impedance for each electrode k may be expressed as [16]:

$$\begin{align} Z_{\textrm{diff,sp,k}} = U^{^{\prime}}_k(c_{k0}) c_k^{\max} \frac{\tau_k^d}{3 Q_k} \frac{\tanh{\left(\sqrt{s \tau_k^d}\right)}}{\tanh{\left(\sqrt{s \tau_k^d}\right)} - \sqrt{s \tau_k^d}} , \end{align} \tag{ 43 }$$

where s is the Laplace variable, $U^{^{\prime}}_k(c_{k0})$ is the gradient of half cell OCV, with respect to concentration, $\tau_k^d$ is a diffusion timescale ${r_k^2}/{D_k}$ for each electrode, and $Q_k = F A L_k \varepsilon_{\textrm{act},k} c^{\max}_{k}$ is the maximum capacity of each electrode.

Levi and Aurbach [136] demonstrated good agreement between the diffusivity values obtained by PITT and EIS. Ecker et al [17] demonstrated good agreement between the diffusion values obtained by GITT and EIS, except for the highly lithiated phases of graphite in which too few EIS measurements were taken to identify the spikes in $D_\mathrm{n}$ corresponding to the phase transitions. Further comparison of the techniques for characterising diffusion coefficients can be found in Deng and Lu [192] and Ivanishchev et al [193]. Applying titration techniques for SoC-dependent diffusivities is also particularly time-consuming (ca. 17 days [194]) and EIS may be a more time-efficient alternative. Similar to most parameters with Arrhenius-like responses, the activation energy can be determined by fitting repeat measurements at different temperatures [17, 99].

The ionic conductivity, $\sigma_\mathrm{e}$, quantifies the mobility of ions under an electric field, in the absence of concentration gradients. It appears in the generalised Ohm's law (equation (4)) and determines the ionic current in response to the gradient of the electrochemical potential. Ionic conductivity is a strong function of temperature and concentration. In dilute electrolytes, conductivity rises with salt concentration, as the number of mobile ions in solution increases. As salt content increases past a certain threshold, however, species–species interactions, such as ion pairing, cause $\sigma_\mathrm{e}$ to fall [200]. The dependence of conductivity on salt concentration at a given temperature can be fit well with simple polynomial forms [201] or by phenomenologically inspired forms, such as the square-root dependence of equivalent conductance ($\sigma_\mathrm{e}/(Fc_\mathrm{e})$), predicted by Debye–Hückel limiting behaviour [202–204]. Other empirical functional forms have also been employed, to account for temperature dependence [205, 206].

Ionic conductivities, $\sigma_\mathrm{e}$, for a variety of LiPF₆-containing electrolytes from the property database are presented in figure 11. Generally, ion transport improves with temperature. The temperature dependence of electrolyte conductivity in practical conditions is well-described by a simple Arrhenius relation [17, 200]. Consistent results may be attained by interpolating conductivities using activation energies fit from a range of temperatures. Extrapolating to wider temperature ranges may require compensating for thermal effects on the dielectric constant and viscosity [207, 208].

Zoom In Zoom Out Reset image size

In LiBs, the electrolytic solution permeates porous phases in the separator and electrodes. Therefore, before they can be utilised in the DFN model, effective transport parameters must be obtained by scaling with a transport efficiency factor, $\mathcal{B}$. These extrinsic details were discussed in section 3.1.5.

The salt diffusivity, $D_\mathrm{e}$, determines the contribution to the ionic flux by salt concentration gradients, as given by equation (2). The thermodynamic diffusivity, $\mathscr{D}$, also describes diffusion under the true driving force of chemical potential gradients [196]. A direct conversion between $\mathscr{D}$ and $D_\mathrm{e}$—the value usually experimentally measured—exists via the thermodynamic factor, $1+\partial\ln f_{\pm}/\partial\ln c_{\textrm{e}}$ (discussed below in section 5.1.4):

$$\begin{align} D_\mathrm{e} = c_\textrm{T}\overline{V}_0\left(1+\frac{\partial\ln f_{\pm}}{\partial\ln c_\mathrm{e}}\right)\mathscr{D}. \end{align} \tag{ 44 }$$

The additional term, $c_\textrm{T}\overline{V}_0$, which is the total molarity, $c_\textrm{T}$, multiplied by the solvent partial molar volume, $\overline{V}_0$, arises as a drift factor that puts $D_\mathrm{e}$ in reference to the solvent velocity. It is this diffusion coefficient, $D_\mathrm{e}$, which most commonly appears in DFN model equation (2). Diffusion coefficients have been demonstrated to decrease with increasing salt concentration and solvent viscosity, as measured by a range of techniques. The rate of diffusion also increases with temperature. Figure 12 shows the salt diffusivities, $D_\mathrm{e}$, reported in the literature across varying solvent compositions and salt concentrations. Diffusion coefficients are typically well fitted by simple polynomial expressions or exponentials.

Zoom In Zoom Out Reset image size

Salt diffusion coefficients can be measured experimentally (see restricted diffusion below, in subsection 5.2.2, [204, 209]). Models that use a constant, average diffusivity will inherently force concentration gradients in solution to develop symmetrically [210]. This may cause an over- or underestimation of the concentration polarisation. Many modeling efforts use approximate diffusion coefficients. Schmalstieg et al [99] use the Nernst–Einstein relationship, derived from ideal, infinitely dilute assumptions, to estimate $D_\mathrm{e}$ from $\sigma_\mathrm{e}$. The Nernst-Einstein model is not well justified for concentrated solutions, however, and deviations from this ideal behaviour occur well below the salt concentrations of around 1 M that are typical in LiB electrolytes [197]. Self-diffusivities attained by nuclear magnetic resonance (NMR) or molecular dynamics, while similar in magnitude to $D_\mathrm{e}$, require additional interpretation to be used consistently in the DFN model [211, 212].

The choice of model used to extract the diffusion coefficient from experimental data strongly influences the inferred diffusivity. For example, following Fuller et al [6], the ion flux (equation (2) in our macroscopic relationships) is composed of contributions from migration and diffusion, with convection assumed to be negligible. The majority of DFN battery models have adopted this convention. Indeed, if convection, driven by the density changes which accompany concentration changes, is accounted for, as discussed by the works of Liu and Monroe [213, 214], this may skew the apparent diffusion coefficient and transference number.

The transference number, $t_+^0$, is the fraction of the conductivity contributed by cation motion. This parameter is typically reported relative to a reference frame set by the solvent velocity, which is taken to be zero in conventional DFN models. Cation and anion transference number must sum to unity. From figure 13 it is clear that the lithium ion carries a minority of the charge during typical battery operation. Many battery models choose to use a simplified average $t_+^0$ value in the 0.3–0.4 range.

Zoom In Zoom Out Reset image size

A low transference number implies poor power efficiency. Diederichsen et al [215] demonstrated, via a sensitivity analysis on the DFN model, that raising the transference number by 0.2 can have significant impact on accessible capacity during charge and discharge. The limiting current—the maximum current at which charge carriers are completely depleted at an electrode—is also a strong function of transference and diffusivity.

While it is clear that accuracy of the transference number has a strong impact on the accuracy of battery models at high C-rate, consensus is low on the expected magnitudes and composition dependencies of $t_+^0$. The varied experimental approaches and processing techniques, using dilute or concentrated solution theory, as discussed in section 5.2, also contribute to the spread in the data across all electrolyte formulations and salt compositions [216]. For example, experimentally obtained negative cation transference numbers have been rationalised by the presence of charged or neutral aggregates [217], or by moving boundaries due to lithium plating [218]. The choice of reference velocity (e.g. solvent or volume-averaged) used can also impact apparent transference [196, 213].

The ultimate purpose of these electrolyte parameters is to predict the lithium ion concentration, as a function of position and time. It is important to reiterate that the specific formulation of a DFN model will impact the interpretation of experiments designed to measure $t_+^0$. For example, models that incorporate dendrite growth [219] will predict vastly different transference numbers to models with ion pairing [218]. Similarly, models that do not include solute-volume effects [213] will overpredict transference when convective-polarisation suppression is wrongly accounted for as a migration effect [204, 209]. Finally, it should be borne in mind that, at low C-rates, simplifying assumptions about $t_+^0$ have minimal impact on simulation validity, but that these assumptions may be violated at higher C-rates.

The thermodynamic factor, $(1+\partial\ln f_{\pm}/\partial\ln c_{\textrm{e}})$, is key to translating a concentration difference into an OCV drop, as can be seen from the generalised form of Ohm's law in equation (4). Thermodynamic factors measure solution non-ideality, accounting for deviations from Nernstian behaviour. The thermodynamic factor is a grouping of terms which arises when the derivative of the chemical potential, $\mu_\mathrm{e}$, with respect to concentration, is computed [195]. This leads to the appearance of the derivative of the mean molar salt activity coefficient, $f_{\pm}$, with respect to salt molarity, $c_\mathrm{e}$. In the literature, the gradient of the activity coefficient is also often defined relative to composition in a molal ($m_\mathrm{e}$) basis, which can be mapped from the molar basis via the mean molal activity coefficient, $\gamma_\pm$, solvent molarity, c₀, and partial molar volume of solvent, $\overline{V}_0$:

$$\begin{align} 1+\frac{\partial\ln f_{\pm}}{\partial\ln c_\mathrm{e}} = \frac{1}{c_0\overline{V}_0}\left(1+\frac{\partial \ln{\gamma_\pm}}{\partial \ln{m}}\right) = \frac{1}{RT}\frac{\partial \mu_\mathrm{e}}{\partial \ln{c_\mathrm{e}}}. \end{align} \tag{ 45 }$$

By definition, the thermodynamic factor approaches unity at infinite dilution. Debye-Hückel theory predicts a decrease in thermodynamic factor with respect to concentration near infinite dilution [202, 203, 220], but this concentration range is far from the practical concentrations used in lithium-ion electrolytes. The LiB literature, which focuses on moderately dilute to concentrated electrolytes, typically reports $1+\partial\ln f_{\pm}/\partial\ln c_{\textrm{e}}$ values above unity that increase across the observed concentration range, in agreement with a condition where ion–solvent or ion–ion interactions are significant [209]. There also appears to be a weaker dependence on temperature than for other parameters [206]. Polynomials or power series in $\sqrt{c}$, following Debye-Hückel theory, have been used to fit correlations for $1+\partial\ln f_{\pm}/\partial\ln c_{\textrm{e}}$ [204, 209].

The measurements summarised in figure 14 indicate that conventional lithium-ion electrolytes with compositions above 0.5 M should not be treated as ideal [206, 209, 221]. The original DFN publications included $1+\partial\ln f_{\pm}/\partial\ln c_{\textrm{e}}$ in their model formulation, but assumed ideality because thermodynamic data was unavailable [5, 6]. Many subsequent DFN models for conventional LiBs have dropped $1+\partial\ln f_{\pm}/\partial\ln c_{\textrm{e}}$ [17, 19, 99]. Inspection of the DFN model formulation reveals that neglect of the thermodynamic factor has no effect on compositions predicted by equations (1) and (2), but exclusion of $1+\partial\ln f_{\pm}/\partial\ln c_{\textrm{e}}$ will introduce errors in the concentration overpotential, as seen in equation (4). In modelling, this will impact predicted battery capacity, due to the propensity to underpredict cell-voltage cutoffs during charge/discharge. Ideal assumptions only have a strong influence at medium-high C-rates, however, when ohmic losses are not dominant. In situations where battery operation induces a high degree of concentration polarisation, or for cell formulations with particularly high electrolyte molarities, inclusion of an accurate $1+\partial\ln f_{\pm}/\partial\ln c_{\textrm{e}}$ is warranted.

Zoom In Zoom Out Reset image size

Partial molar volumes are thermodynamic properties that describe how volume is distributed across the components in solution ( $\overline{V}_\mathrm{e}$ for salt and $\overline{V}_0$ for solvent). Doyle et al [5] simplified their treatment of concentrated solution theory by ascribing zero molar volume to the dissolved salt in solution. Partial molar volumes have since been included in DFN-type models [201, 209, 213, 222].

In the one-dimensional ion conservation equation (1), the divergence of flux $\frac{\partial N_\textrm{e}}{\partial x}$ yields an excluded-volume factor $\left(1-\frac{d\ln c_0}{d\ln c_\mathrm{e}}\right)$. As solvent molarity c₀ is typically a weak function of salt content, $c_\mathrm{e}$, in typical lithium-ion electrolytes, Newman [195] use this assumption to neglect the excluded-volume factor, hence setting $\overline{V}_\mathrm{e}$ to zero. Partial molar volumes relate to the excluded-volume factor through the solvent concentration, c₀, and solution density, ρ [195]:

$$\begin{align} 1-\frac{d\ln c_0}{d\ln c_\mathrm{e}} = \frac{\rho}{c_0~M_0}\left(1-\frac{d\ln\rho}{d\ln c_\mathrm{e}}\right) = \frac{1}{c_0\overline{V}_0}, \end{align} \tag{ 46 }$$

with M₀ being the solvent molar mass. One can measure $\overline{V}_\mathrm{e}$ and $\overline{V}_0$ readily with densitometry, as discussed in section 5.2.5. Recent spectroscopic techniques have also allowed for the probing of solvent concentrations [223, 224].

We have previously discussed the use of EIS in inferring geometrical and electrode intrinsic parameters in sections 3.2.2 and 4.2.6, respectively. Here, we cover its use for elucidating electrolyte intrinsic properties. Commercially-available conductivity probes which operate on frequency-based principles are commonly used to measure ionic conductivities for many common lithium-ion electrolyte formulations [200, 201, 205, 209, 225–227].

EIS is also performed with symmetric planar cells with blocking (stainless steel) or non-blocking (lithium) electrodes. The ionic conductivity, $\sigma_\mathrm{e}$, may be extracted by measuring the high-frequency resistance value, $R_\textrm{b}$, in the bulk electrolyte, relative to the EIS cell geometry, and then leveraging the equation:

$$\begin{align} \sigma_\mathrm{e} = \frac{L_\textrm{s}}{R_\textrm{b}A}, \end{align} \tag{ 47 }$$

where $L_\textrm{s}$ is the planar inter-electrode spacing and A is the geometrical electrode area. This is often done with an electrolyte soaked porous separator between electrodes and can be used to calculate the effective conductivity, based on the high frequency bulk resistance intercept on a Nyquist plot [204, 228, 229]. The low-frequency impedance response can also determine contributions from the ionic diffusion coefficient, $D_\mathrm{e}(c_\mathrm{e})$, and transference number, $t_+^0(c_\mathrm{e})$ [228, 230]. This approach is akin to a frequency-domain analysis of repeated polarisation-relaxation experiments, discussed in the following section.

Polarisation-relaxation cells refer to symmetric planar lithium metal electrode cells, used to characterise electrolyte properties by performing a sequence of steps where a current or voltage (or lack thereof) is applied. The concentration gradients resulting from polarisation-relaxation steps are then tracked, often being inferred through concentration overpotentials (also called diffusion potentials), which scale with the concentration gradient [197]. The next few sections outline common polarisation-relaxation experimental techniques used for electrolyte parameterisation.

5.2.2.1. Restricted diffusion

The electrolyte salt diffusion coefficient, $D_\mathrm{e}(c_\mathrm{e})$, is best measured with a restricted-diffusion experiment in a sealed vertical cell, where the ion flux directly opposes the direction of gravity. Orientation is important, especially when no porous separator is employed to prevent natural convection [195] or buoyancy effects [231] from marring the method. During restricted diffusion, a nonuniform concentration profile is induced across an electrolyte, often galvanostatically or potentiostatically, before monitoring the relaxation of the concentration difference between two fixed points. As concentration gradients relate to the OCV via equation (4), the diffusion coefficient can be inferred from the exponential relaxation of the cell voltage [196, 204–206, 209, 232]:

$$\begin{align} \ln V_\textrm{OCV} \propto \Delta c_\mathrm{e} = -\frac{\pi^2 D_\mathrm{e}^\textrm{effective}}{L_\textrm{s}^2} \cdot t + \textrm{const}, \end{align} \tag{ 48 }$$

where $V_\textrm{OCV}$ is the restricted diffusion cell voltage measured during the open-circuit relaxation, $\Delta c_\mathrm{e}$ is the difference in concentration between the points where $V_\textrm{OCV}$ is measured, $L_\textrm{s}$ is the separator or bulk electrolyte thickness and t is time. A formal analysis by Newman and Chapman [196] showed that this linear time dependence of the logarithm of voltage can be used to process data accurately, even when diffusivity varies with respect to concentration. Thus, when concentration differences are sufficiently small, Fickian diffusion with a constant diffusion coefficient can be assumed and the diffusivity resulting from a restricted-diffusion measurement should be interpreted as the value at the electrolyte's equilibrium concentration.

Ehrl et al [233] found that the restricted-diffusion method through steady-state polarisation and long-term relaxation was less susceptible to the influence of double-layer relaxation effects than pulse-polarisation methods. Error due to surface effects may also be avoided by using a four-electrode cell, containing two reference electrodes positioned at some distance from the reactive metal interfaces, to track diffusivities using cell voltages, as presented by Farkhondeh et al [234].

When performing restricted diffusion measurements, one should ensure that the inter-electrode length, $L_\textrm{s}$, between electrodes and duration of voltage relaxation, t, are such that the voltage relaxes linearly. Other indicators of composition such as conductance and spectroscopic response can also been used to track concentration gradients between two points [196, 221, 235, 236].

5.2.2.2. Galvanostatic polarisation

The galvanostatic polarisation, or current interrupt method, originally implemented by Ma et al [7] for polymer electrolytes, can be used to measure the transference number. This analysis is based on semi-infinite linear diffusion [160], where the concentration gradient or OCV measured immediately after a galvanostatic pulse is directly related to the transport and thermodynamic properties of the electrolyte. By varying the polarisation time, $t_\textrm{i}$, and current density, $i_\mathrm{circuit}$, the slope, m, of cell potential, V, after the moment of current cutoff with respect to a transformed time variable $i_\mathrm{circuit}t_\textrm{i}^{1/2}$, can be related to the transference number via:

$$\begin{align} t_+^0 = 1-\frac{mFc_{\mathrm{e}0}\sqrt{\pi D_\mathrm{e}}}{4\frac{dV}{d\ln c_\mathrm{e}}}. \end{align} \tag{ 49 }$$

Here F is the Faraday constant, $c_{\mathrm{e}0}$ is the bulk salt concentration, $D_\mathrm{e}$ is the diffusivity and $\frac{dV}{d\ln c_\mathrm{e}}$ is the concentration dependence of the derivative of the measured voltage (the diffusion potential) with respect to the logarithm of the concentration. Note that this method requires independent concentration cell and diffusion coefficient measurements to determine $D_\mathrm{e}$ and $\frac{dV}{d\ln c_\mathrm{e}}$, and so care must be taken to preclude the propagation of experimental errors. While valid for non-ideal electrolytes, this relationship is derived using concentrated solution theory, assuming a binary electrolyte with electrode reactions under anion-blocking conditions. Several groups have used the galvanostatic polarisation method to measure transport properties in non-aqueous LiB electrolytes [201, 206, 216, 233, 237]. Farkhondeh et al [234] have also fitted their four-electrode galvanostatic-pulse relaxation experiments with a transport model, to extract electrolyte transport properties.

5.2.2.3. Potentiostatic polarisation

The potentiostatic polarisation, or steady-state direct current, method was popularised by Bruce and Vincent [238] for transference measurements in polymer electrolytes. The property measured by this method is better described as a current fraction [239] or transport number, $t_+$ [204, 209, 225], since infinite dilution is assumed in the original derivation of equation (50)—implying that ionic species are non-interacting. It should be emphasised that the interpretation of parameters $t_+$ and $t_+^0$ is distinct. Analysis by Bergstrom et al [232] has shown that potentiostatic polarisation measurements with liquid electrolytes are highly susceptible to the instability of lithium metal electrode surfaces. As such, $t_+$ cannot be used to parameterise models that accurately predict ion flux through the electrolyte phase.

The original Bruce–Vincent method applied a small constant potential bias (10 mV) and compared the initial and steady-state currents, after subtracting interfacial resistance effects, to determine the transference number via:

$$\begin{align} t_+ = \frac{I_\textrm{ss}(V - I_\textrm{i} R_\textrm{i})}{I_\textrm{i}(V - I_\textrm{ss} R_\textrm{ss})}. \end{align} \tag{ 50 }$$

Here, $I_\textrm{i}$ is the initial current measured when a bias voltage, V, is applied, and $i_\textrm{ss}$ is the steady-state current density. $R_\textrm{i}$ and $R_\textrm{ss}$ are the interfacial resistances (measured with EIS) from the electrodes, before and after polarisation. Hou and Monroe [204] have demonstrated that the propensity to underpredict concentration polarisation due to $t_+$ may be corrected for by accounting for the partial molar volume of salt. Recently, Balsara and Newman [240] have rederived the relationship between the steady-state current and transference number from concentrated solution theory, injecting the dependence of the diffusion coefficient and thermodynamic factor for higher salt compositions. Comparisons between galvanostatic and potentiostatic polarisation methods for transference numbers have been provided by Zugmann et al [216] and Ehrl et al [233].

5.2.2.4. Hittorf method

The Hittorf method measures the change in moles of cations, compared to anions, across the cell electrolyte cavity, with the passage of a set amount of charge [220, 241]. This ratio is then directly related to the transference number. Hou and Monroe [204] designed a densitometric Hittorf cell with stopcocks, allowing the isolation of anodic and cathodic chambers after polarisation. Once gradients in the isolated chambers relax, the changes in density between chambers could be used to determine the concentration difference. The transference number is related to the concentration difference via:

$$\begin{align} t_+^0 = 1 + \frac{FV_\textrm{chamber}\Delta c_\mathrm{e}}{Q(1-\overline{V}_\mathrm{e} c_\mathrm{e})}, \end{align} \tag{ 51 }$$

where $V_\textrm{chamber}$ is the volume of the cell chamber where ion accumulation or depletion occurs, $\Delta c_\mathrm{e}$ is the final change in concentration between the neutral chamber and cathodic or anodic chambers, Q is the charge passed and $\overline{V}_\mathrm{e}$ is the partial molar volume of salt attained from the density-molarity correlation. The overall change in composition from the Hittorf method may also by tracked conductometrically and potentiometrically [205].