School of Mathematics and Physics, University of Portsmouth, Lion Terrace, PO1
3HF, United Kingdom School of Mathematics and Physics, University of
Portsmouth, Lion Terrace, PO1 3HF, United Kingdom

# Learning Optimal Forms of Constitutive Relations Characterizing Ion
Intercalation from Data in Mathematical Models of Lithium-ion Batteries

Lindsey Daniels  Department of Mathematics and Statistics, McMaster
University, Hamilton, Ontario, Canada L8S 4K1  Smita Sahu  The Faraday
Institution, Quad One, Becquerel Avenue, Harwell Campus, Didcot, OX11 0RA,
United Kingdom  Kevin J. Sanders  Department of Chemistry, McMaster
University, Hamilton, Ontario, L8S 4L8, Canada  Gillian R. Goward  Department
of Chemistry, McMaster University, Hamilton, Ontario, L8S 4L8, Canada  Jamie
M. Foster  The Faraday Institution, Quad One, Becquerel Avenue, Harwell
Campus, Didcot, OX11 0RA, United Kingdom  Bartosz Protas  Department of
Mathematics and Statistics, McMaster University, Hamilton, Ontario, Canada L8S
4K1<EMAIL_ADDRESS>*

###### Abstract

Most mathematical models of the transport of charged species in battery
electrodes require a constitutive relation describing intercalation of
Lithium, which is a reversible process taking place on the interface between
the electrolyte and active particle. The most commonly used model is the
Butler-Volmer relation, which gives the current density as a product of two
expressions: one, the exchange current, depends on Lithium concentration only
whereas the other expression depends on both Lithium concentration and on the
overpotential. We consider an inverse problem where an optimal form of the
exchange current density is inferred, subject to minimum assumptions, from
experimental voltage curves. This inverse problem is recast as an optimization
problem in which the least-squares error functional is minimized with a
suitable Sobolev gradient approach. The proposed method is thoroughly
validated and we also quantify the reconstruction uncertainty. Finally, we
identify the universal features in the constitutive relations inferred from
data obtained during charging and discharging at different C-rates and discuss
how these features differ from the behaviour predicted by the standard Butler-
Volmer relation. We also identify possible limitations of the proposed
approach, mostly related to uncertainties inherent in the material properties
assumed known in the inverse problem. Our approach can be used to
systematically improve the accuracy of mathematical models employed to
describe Li-ion batteries as well as other systems relying on the Butler-
Volmer relation.

##  1 Introduction

The modelling of many physical processes relies on fundamental understanding
of the interface between two interacting materials (or phases). It may also
require the knowledge of specific material properties, which can be
challenging to obtain or measure 1. Interfacial models can be complex due to
the coupling of different physical processes at the interface which makes
probing their parametric dependence more difficult. In an effort to aid in the
parametrization challenge, we propose an inverse modelling technique to
optimally reconstruct constitutive relationships describing an important class
of interfacial processes from data. Inverse modelling is an approach used to
learn optimal forms of material properties from data by combining a
mathematical model of the system with experimental measurements.

The present work proposes a computational approach alternative to machine
learning techniques which have recently become popular2, 3, 4. Instead, we use
tools from calculus of variations to design a numerical procedure to recover
an optimal form of the constitutive relation, where “optimal” means that the
reconstructed constitutive relation achieves the best possible fit to the
data. Such methods have already proved useful for solving a range of inverse
problems in electrochemistry, which is an area abounding with systems whose
material properties are notoriously difficult to measure directly5, 6, 7, 8,
9, 10, 1. Many of these applications are driven by recent advances in the
mathematical modelling of transport processes in Lithium-ion (Li-ion)
batteries, which are crucial to the development of electric vehicles and
customer electronics11, 12. While most of the work on inverse modelling in
electrochemistry has focused on inferring material properties defined in the
bulk of the electrolyte or active material, we extend this approach to
constitutive relations characterizing interfacial phenomena. More
specifically, we focus on constitutive relations describing the intercalation
of ions from the electrolyte into the active material, which play a central
role in modelling electrochemical systems13, but also find applications in
other area of chemistry and physics as well as in biology. For example, such
interfacial models have been used to relate the current density to the
overpotential in bio-anode polarization curve14, to calculate the reaction
rate from voltammograms15, to describe the electrode dynamics in
supercapacitors 16, and to describe the electrochemistry of a neuron 17.

An inverse problem is usually formulated to minimize the least-squares error
between the values predicted by the model and the experimentally measured
quantities. When studying a distributed system (with dependence on space
and/or time) and when the reconstructed property is a function, rather than a
set of scalars, there are two main types of inverse problems. In the more
traditional problems, the reconstructed material property is assumed to be a
function of independent variables, usually space coordinate and/or time18, 19,
20. On the other hand, in less standard modern inverse problems21, 6, 22, 23,
9, 10, the inferred properties are functions of the dependent variables in the
system. The interfacial inverse problem we consider here belongs to this
second category.

Inverse problems tend to be ill-posed, as they usually do not admit a unique
exact solution, but rather a large number (often infinite) of approximate
solutions19. Variation in experimental measurements together with weak
dependence of the model on some parameters can result in significant changes
in the reconstructed solution which can be amplified by experimental and
numerical errors. Yet, inverse modelling can help to elucidate optimal forms
of constitutive relations and assess the validity of different models, as
solutions violating thermodynamic principles can point to potential
inconsistencies in the models used7, 9.

The structure of the paper is as follows: in the next section we discuss the
Butler-Volmer model; in Section 3, the experimental setup and the measurement
data are introduced; Section 4 details the mathematical framework and
optimization approach, and outlines the sensitivity analysis; the
reconstruction results obtained with the inverse modelling approach are
presented in Section 5; and finally, the conclusions of this work are
presented in Section 6, whereas some additional technical material is
collected in three appendices.

##  2 The Butler-Volmer Relation and its Generalization

Below we state the Butler-Volmer relation in a form that is commonly used,
whereas the multiplicity of forms in which this relation is encountered is
discussed in a recent review24. A complete derivation of the Butler-Volmer
relation can be found in the literature25, 26, 13. Hereafter, the
subscripts/superscripts “s” and “e” will denote quantities defined in the
solid phase and in the electrolyte, respectively.

The net current density $j$ is the difference between the backward and forward
current densities at the interface, and the dependence of this quantity on
electrochemical and thermodynamic parameters is given by the Butler-Volmer
relation

$j=Fk_{0}c_{\rm{e}}^{\alpha_{\rm{a},m}}\left(c_{{\rm{s}},i,{\rm{max}}}-c_{{\rm{s}},i}\right)^{\alpha_{\rm{a},m}}c_{{\rm{s}},i}^{\alpha_{\rm{c},m}}\,\left[e^{\frac{\alpha_{\rm{a},m}F}{RT}\eta}-e^{\frac{-\alpha_{\rm{c},m}F}{RT}\eta}\right],\quad\textrm{where}\quad\eta=\phi_{\rm{s}}-\phi_{\rm{e}}-{U_{\rm{eq}}^{i}}(c_{{\rm{s}},i}),$
(1)

in which $\eta$ is the overpotential, $c_{{\rm{s}},i}$ and $c_{\rm{e}}$ are
the Lithium concentrations in the solid phase and in the electrolyte with
$i={a,c}$ indicating the anode or cathode respectively, $\phi_{\rm{s}}$ and
$\phi_{\rm{e}}$ are the potentials in the solid phase and in the electrolyte
(both measured with respect to a Li reference electrode), whereas
${U_{\rm{eq}}^{i}}(c_{{\rm{s}},i})$ is the concentration-dependent equilibrium
potential of the electrode, $\alpha_{{a},m}$ and
$\alpha_{{c},m}=1-\alpha_{{a},m}$ are the anodic and cathodic charge transfer
coefficients with $m=-,+$ indicating the anode or cathode respectively (we do
not use $a$ and $c$ here to avoid repeated identical indices), $k_{0}$ is the
reaction rate constant, and $c_{{\rm{s}},i,{\rm{max}}}$ is the maximum Lithium
concentration in the solid. As is evident from (1), the Butler-Volmer relation
can be naturally split into an algebraic factor depending on the Lithium
concentration only, typically referred to as the exchange current, and the
remaining exponential factor depending on both the potential difference
$\Delta\phi=\phi_{\rm{s}}-\phi_{\rm{e}}$ and the Lithium concentration in the
electrolyte, the latter dependence through the equilibrium potential
${U_{\rm{eq}}^{i}}(c_{{\rm{s}},i})$. This allows us to factor the Butler-
Volmer relation (1) as

$\displaystyle j$
$\displaystyle=i_{0}(c_{{\rm{s}},i})\,\psi(\Delta\phi,c_{{\rm{s}},i}),\qquad\textrm{where}$
(2a) $\displaystyle i_{0}(c_{{\rm{s}},i})$
$\displaystyle=Fk_{0}c_{\rm{e}}^{\alpha_{{a},m}}\left(c_{{\rm{s}},i,{\rm{max}}}-c_{{\rm{s}},i}\right)^{\alpha_{{a},m}}c_{{\rm{s}},i}^{\alpha_{{c},m}},$
(2b) $\displaystyle\psi(\Delta\phi,c_{{\rm{s}},i})$
$\displaystyle=\exp\left[\frac{\alpha_{{a},m}F}{RT}\left(\Delta\phi-{U_{\rm{eq}}^{i}}(c_{{\rm{s}},i})\right)\right]-\exp\left[-\frac{\alpha_{{c},m}F}{RT}\left(\Delta\phi-{U_{\rm{eq}}^{i}}(c_{{\rm{s}},i})\right)\right].$
(2c)

While the exponential term (2c) has had a long and well-documented history in
electrochemistry 27, 25, 26, additional assumptions are required to justify
the form of the concentration-dependent algebraic factor (2b)24. As will be
discussed in more detail below, we focus on the dependence of the exchange
current on $c_{{\rm{s}},i}$ and its dependence on $c_{\rm{e}}$ in (2b) will be
omitted so that our results are most applicable to situations in which the
electrolyte concentration is approximately uniform (e.g., at relatively low
currents or for thin electrodes).

In the present study, we thus consider a generalization of the Butler-Volmer
relation (2) given

$j=\iota(c_{{\rm{s}},i})\,\psi(\Delta\phi,c_{{\rm{s}},i}),$ (3)

where $\psi(\Delta\phi,c_{{\rm{s}},i})$ is assumed known and given by (2c),
whereas the factor $i_{0}(c_{{\rm{s}},i})$, cf. (2b), is replaced by a general
function $\iota(c_{{\rm{s}},i})$ whose optimal (in a mathematically precise
sense to be specified below) form will be inferred from standard
electrochemical measurements. This inverse problem will be solved using a
computationally robust and thermodynamically consistent numerical optimization
procedure where optimal functional forms of the factor $\iota(c_{{\rm{s}},i})$
will be sought in a very general form subject to minimum only assumptions.
These assumptions are28:

  * A1.

the function $\iota=\iota(c_{{\rm{s}},i})$ must be sufficiently smooth (this
will be made precise in Section 4.2),

  * A2.

when the Lithium concentration in the electrode goes to zero
($c_{{\rm{s}},i}\rightarrow 0$), the Lithium flux out of the active material
must vanish to prevent nonphysical negative Lithium concentrations in the
electrode, and

  * A3.

when the electrode is saturated ($c_{{\rm{s}},i}\rightarrow
c_{{\rm{s}},i,{\rm{max}}}$), the Lithium flux into the active material must
also vanish to prevent nonphysical overfilling of Lithium.

These represent constraints on the behaviour of the factor
$\iota(c_{{\rm{s}},i})$ for small ($c_{{\rm{s}},i}\rightarrow 0$) and large
($c_{{\rm{s}},i}\rightarrow c_{{\rm{s}},i,{\rm{max}}}$) concentrations. We
note that this limiting behaviour is identical as in the standard Butler-
Volmer model, cf. (2b). However, since the equilibrium potential
$U_{\rm{eq}}^{i}(c_{{\rm{s}},i})$ is unbounded in these limits, cf. Supporting
Information (Appendix A), the minimum and maximum concentrations of 0 and
$c_{{\rm{s}},i,{\rm{max}}}$ are in practice never attained28. Application of
the proposed approach to voltage measurements performed during charge and
discharge experiments produces optimal forms of the exchange current
$\iota(c_{{\rm{s}},i})$, which are consistent between charge and discharge and
reveal features making them quite different from expression (2b). We also
characterize the uncertainty of the reconstructions and compare our results to
a simpler approach where the inverse problem is solved by calibrating a small
number of scalar parameters in expression (2b).

##  3 Experimental Datasets

In the formulation of our inverse problem, we will use several datasets from
experiments performed on Li-ion batteries under galvanostatic conditions
during both charge and discharge at different rates29, 30. Below, we first
describe the two experimental setups we considered and present the data
obtained. Then, we demonstrate that, for all datasets, the fraction of the
signal (the time-dependent voltage) containing information about the
interfacial reactions is sufficiently large both to be used in inverse
modelling and for the main assumptions underlying our mathematical model to be
satisfied. Within the electrolyte, there is a potential drop between the anode
and the cathode, which is accounted for when solving the complete Doyle-
Fuller-Newman (DFN) model13. However, this potential drop is often either
neglected in the standard Single-Particle model (SPM), cf. Section 4.1, based
on the assumption that it is small when compared to the contributions from the
electrode kinetics or is accounted for through a linearized model based on the
assumption that the electrolyte can be considered ohmic if the current is
small31, 10. The SPM is typically an accurate reduction of the DFN model for
C-rates up to approximately $1$C 10 and given this well-known limitation of
the SPM, we focus here on small and moderate charge/discharge rates only.

In order to assess whether the experimental datasets are consistent with the
SPM, we determine the corresponding time-dependent concentration and potential
profiles in the cells based on solutions of the complete DFN model13 which are
obtained using the suite DandeLiion32. This allows us to determine the
relevant voltage drops for each component of the cell. Since the potential
drop across the electrolyte is a function of both space (i.e., the distance
$x$ from the anode) and the (dis)charge time, we first compute the total drop
as $\phi_{\rm{e}}(L,t)-\phi_{\rm{e}}(0,t)$, where $L$ is the total length of
the cell thereby taking a maximal estimate of the potential drop across the
electrolyte, and then take the time average.

###  3.1 Slow Charge Rates

For slow charge rates, we consider two datasets corresponding to C-rates of
C/10 and C/3. These datasets were obtained using an electrochemical cell
constructed from an unused commercial Lithium-ion battery containing a double-
sided coated graphite anode and LiNi0.6Mn0.2Co0.2O2 (NMC$622$) cathode 29. The
commercial battery was disassembled, the electrodes were washed, cut, and
dried, and reassembled in a replaceable cartridge-like battery container that
mimics the single-layer pouch cell configuration. The assembled cell consists
of a copper (Cu) current collector plate, graphite anode, Celgard $2325$
separator, NMC$622$ cathode, and an aluminum (Al) current collector plate. The
assembled cell was rested for $24$ hours before being subject to cycling
protocols and experimental measurements. Figure 1 shows each individual
component of the cell before assembly. Full experimental details and set up
can be found by Sanders et al 29.

Figure 1: Image of the cell assembly. From left to right: bottom part of the
cartridge with copper (Cu) metal current collector plate, graphite anode,
Celgard $2325$ separator, NMC$622$ cathode, and top part of the cartridge with
aluminum (Al) current collector plate. The EPDM o-ring used to provide a gas-
tight seal and $8$ M$2$ PEEK screws used to close the cartridge are also shown
around the top part on the right-hand side.

In Figure 2, the experimental voltage curves for the anode and cathode
equilibrium voltages are plotted against the normalized concentration in the
two electrodes. The graphite voltage curve in Figure 2a was obtained by
charging a half-cell with a graphite electrode to the maximum voltage at a
rate of C$/140$. The voltage curve for the NMC$622$ electrode in Figure 2b was
obtained by charging a half-cell with an NMC$622$ electrode to $4.2$V at a
rate of C$/140$. These measurements will be used as the equilibrium potentials
in our model. The voltage curves obtained with the C-rates of C$/10$ and C$/3$
will be shown in Section 5 where we discuss our results.

Figure 2: The equilibrium voltage curve for (a) the graphite anode and (b) for
the NMC$622$ cathode in the cell used in the experiments described in Section
3.129.

Table 1 details the percentage breakdown of the time-averaged voltage drops
obtained by solving numerically the complete DFN model using the material
parameters listed in Tables 5 and 6 in Supporting Information (Appendix C.1).
In all cases, the overpotentials $\eta_{a}$ and $\eta_{c}$ have approximately
the same share of the total potential as the time-averaged drop across the
electrolyte. Table 1 also shows the ratios between the different potential
drops corresponding to the two C-rates considered. We see that for all
potential drops (except for ${V_{c}-V_{a}}$, as expected) these ratios are
close to the ratio of the two C-rates, which is $10/3$, confirming that the
response is close to linear at small C-rates.

Charge rate | $V_{c}-V_{a}$ | $\eta_{a}$ | $\eta_{c}$ | Electrolyte
---|---|---|---|---
C/$10$ | $99.1755$ | $0.2368$ | $0.3344$ | $0.2533$
C/$3$ | $97.1927$ | $1.0208$ | $1.0195$ | $0.7670$
Ratio |  |  |  |
C/$3$:C/$10$ | $0.9800$ | $4.3103$ | $3.0492$ | $3.0281$
Table 1: Percentage breakdowns of the time-averaged potential drops
corresponding to the main components of the cell used in the experiments
described in Section 3.129. The results were obtained by solving the complete
DFN model using the material parameters listed in Tables 5 and 6 in Supporting
Information (Appendix C.1).

###  3.2 Moderate Discharge Rates

Experimental data for discharge rates of C$/2$ and $1$C were obtained from
Chen et al.30. Their experiments utilized a cylindrical LG M$50$$21700$
(LGM$50$) cell containing an LiNi0.8Mn0.1O1 (NMC$811$) cathode and a graphite
anode. A commercial cell was taken apart and assembled into two half cells to
characterize the respective electrode and cell parameters. The obtained
electrodes were then assembled into a three-electrode full cell configuration
with a Celgard $2325$ separator and subjected to cycling protocols under a
variety of environments. Details can be found in the work by Chen et al.30

Table 2 details the percentage breakdown of the time-averaged voltage drops
associated with the main components of the cell which were obtained by solving
numerically the complete DFN model using the material parameters listed in
Tables 7 and 8 in Supporting Information (Appendix C.2). Table 2 also shows
the ratios between the different potential drops corresponding to the two
C-rates considered. We see that when the C-rate increases from C$/2$ to $1$C,
these ratios are approximately $2$ for all potential drops (except for the
equilibrium potential ${V_{c}-V_{a}}$, as expected), which is consistent with
the linear behaviour expected at low C-rates. Chen et al.30 also obtained
results at a higher discharge rate of 2C, but since this regime is outside the
range of validity of the SPM, we will not consider this data in our study.

Charge rate | $V_{c}-V_{a}$ | $\eta_{a}$ | $\eta_{c}$ | Electrolyte
---|---|---|---|---
C/$2$ | $97.2504$ | $1.5916$ | $0.1907$ | $0.9674$
$1$C | $94.8758$ | $2.6239$ | $0.4245$ | $2.0758$
Ratio |  |  |  |
$1$C:C/$2$ | $0.9756$ | $1.6486$ | $2.2262$ | $2.1459$
Table 2: Percentage breakdowns of the time-averaged potential drops
corresponding to the main components of the cell used in the experiments
described in Section 3.230. The results were obtained by solving the complete
DFN model using the material parameters listed in Tables 7 and 8 in Supporting
Information (Appendix C.2).

We remark that some experiments have made use of thin electrodes to ensure
that the potential drop in the electrolyte is negligible33, 34. In our case,
we emphasize that in both the slow and moderate charge/discharge rate data we
see a non-negligible potential drop in the electrolyte, which is on the order
of the contributions from overpotentials $\eta_{a}$ and $\eta_{c}$, cf. Tables
1 and 2. Thus, an SPM with a simple ohmic model accounting for the electrolyte
as discussed in the next section will be suitable for our purposes.

##  4 Mathematical Model

In this section, we first introduce the SPM framework, then formulate our
inverse problem, describe its numerical solution, and finally discuss the
sensitivity analysis. The SPM framework relies on several simplifying
assumptions used to describe the construction of the Li-ion battery. Namely,
the microstructure of the cell in each electrode is represented in terms of a
representative spherical particle, the electrodes are spatially homogeneous
and the intercalation reaction is uniform across the electrode35, 36, 13. In
practice, these assumptions may not reflect the true behaviour of the cell.
For example, during the manufacturing process, the active particles are
usually covered with binder containing carbon and therefore active material
may not be in a direct contact with the electrolyte37, 38. Electrodes are
typically produced with spatial inhomogeneities such that not all particles
undergo intercalation at the same rate39, 40. In addition, the electrode
particles are not strictly spherical, which can lead to preferential
intercalation in some parts of the electrode41, 42. We note here that some of
these limitations can be circumvented by employing idealized nonporous
geometries43, 44, 45, 46, 47. However, the aforementioned assumptions are
generally invoked when modelling cells under low-power operation (generally up
to a C-rate of $1$C)48, 31, 28, 49, 10. We add that in our study the standard
SPM is augmented to include an Ohmic potential drop in the electrolyte and
features concentration-dependent diffusivity describing transport in the
active particles. As indicated in Section 3, the experimental data utilized in
this work is within the realm of validity for this model.

###  4.1 Single Particle Model (SPM)

A standard description of the transport of Lithium in battery electrodes
usually relies on the DFN model 35, 36, 13. However, this model is
computationally expensive and involves a large number of material parameters,
which makes it difficult to calibrate. Both of these issues are partially
alleviated by various simplified models such as the SPM, which assumes that
the electrode particles are uniformly distributed spheres and that
electrochemical reactions are uniform throughout the electrode31. These
assumptions reduce the model to a transport problem in a single representative
particle. The SPM can be derived as a certain asymptotic limit of the DFN
model and the price to be paid for its simplicity is that it is valid for
small C-rates only — it is usually assumed that it does not provide an
accurate description beyond 1C48, 31, 28, 49, 10. In our study, a suitably
modified SPM is applied to describe transport both in the anode and the
cathode.

In the SPM, the Lithium transport in the solid state is described by a
nonlinear diffusion equation written in spherical coordinates35, 36, 13

$\frac{\partial c_{{\rm{s}},i}}{\partial
t}=\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left[D_{\rm{s}}(c_{{\rm{s}},i})r^{2}\frac{\partial c_{{\rm{s}},i}}{\partial
r}\right]\qquad\textrm{in}\ (0,{R_{i}})\times(0,t_{f}],$ (4)

where $r$ measures the radial distance from the center of the particle outward
to the particle’s surface at ${R_{i}}$, $D_{\rm{s}}{(c_{{\rm{s}},i})}$ is the
(concentration-dependent) diffusivity in the electrodes, $t$ is time ranging
from $0$ to the final time $t_{f}$, and $i={a,c}$ corresponds to the anode and
cathode, respectively. The accompanying boundary conditions are

$\displaystyle\frac{\partial c_{{\rm{s}},i}}{\partial r}\bigg{|}_{r=0}$
$\displaystyle=0$ $\displaystyle\textrm{in}\ (0,t_{f}],$ (5) $\displaystyle-
D_{\rm{s}}(c_{{\rm{s}},i})\frac{\partial c_{{\rm{s}},i}}{\partial
r}\bigg{|}_{r={R_{i}}}$ $\displaystyle=\frac{j_{i}}{F}$
$\displaystyle\textrm{in}\ (0,t_{f}],$ (6)

where $F$ is the Faraday constant, and the initial condition is

$c_{{\rm{s}},i}|_{t=0}=c_{s,i;0},$ (7)

where $c_{s,i;0}$ is the initial concentration of Lithium in the solid. Given
that transport in the electrolyte is not explicitly modelled in the SPM, the
current density $j_{i}$ in the boundary condition (6) can be deduced directly
from the total current $I(t)$ applied to the cell. The conservation of charge
at the anode and cathode implies

$4\pi R_{i}^{2}j_{i}=\pm\frac{I(t)}{n_{i}},$ (8)

in which the total number of particles in the electrode is given by
$n_{i}=A_{i}\,L_{i}\,\epsilon_{i}/(4/3\pi R^{3}_{i})$, where $A_{i}$ is the
cross-sectional area of the electrode, $L_{i}$ is the thickness of the
electrode, and $\epsilon_{i}$ is the volume fraction in the electrode. The
right-hand side is positive for $i={a}$ and negative for $i={c}$. Finally, the
total current to the cell is $I(t)$, with $I(t)<0$ for the anode during
charging, $I(t)>0$ for the cathode during charging, and vice versa during
discharging. Relation (8) ensures that there is enough (de)intercalation
reaction to meet the current demand. We emphasize that, while the Butler-
Volmer relation (2) is not used in the solution of the solution of the SPM, it
is required in order to determine the corresponding overpotential $\eta$ which
is done by post-processing the solution of the SPM. This property will be the
basis for the inverse problem introduced in the next section.

Since in our experimental data the potential drop in the the electrolyte is of
approximately the same magnitude as the overpotentials (cf. Tables 1 and 2),
we augment the expression for the total potential drop predicted by SPM with a
term representing the potential drop in the electrolyte which is assumed to be
ohmic 50 and proportional to the applied current $I$ via a constant
resistivity $R_{el}$. This assumption is supported by the data in Tables 1 and
2 which indicate that the potential drop in the electrolyte does increase
approximately linearly with the applied current. The total potential drop
predicted by the model thus takes the form

$V_{cell}={U_{\rm{eq}}^{c}}(c_{{\rm{s}},c}\big{|}_{r=R_{c}})-{U_{\rm{eq}}^{a}}(c_{{\rm{s}},a}\big{|}_{r=R_{a}})+\eta_{c}-\eta_{a}+R_{el}I,$
(9)

where the additional subscripts $a$ and $c$ denote anode and cathode,
respectively. Relation (9) will serve as the basis for the inverse modelling
introduced in Section 4.2. We note that the SPM is usually defined through
(4)–(8) and will refer to the model resulting from the inclusion of relation
(9) as an “augmented” SPM.

###  4.2 Inverse Problems

Solution of the SPM given above yields concentration profiles
$c_{{\rm{s}},i}(t)$ as functions of time $t$. Within the SPM (4)–(8), the
concentration $c_{{\rm{s}},i}(t)$ does not depend on the constitutive relation
(3), however, the corresponding overpotential does. The idea behind the
inverse problems is thus to infer the optimal form of the factor
$\iota(c_{{\rm{s}},i})$ by minimizing the discrepancy between the
overpotential predicted by (3) and the overpotential deduced from the
experimental data. The total potential measured experimentally as discussed in
Section 3 will be denoted $\overline{V}_{cell}(t)$.

In our first, simpler, inverse problem we assume the constitutive relation is
given by the Butler-Volmer equation (2) and calibrate the scalar parameters
$\alpha_{{a},m}$ and $k_{0}$ parameterizing this relation for both the anode
and the cathode as well as the electrolyte resistivity $R_{el}$ in (9). Since
the total potential is dominated by contributions from the equilibrium
potentials, cf. relation (9) and Tables 1–2, we consider the following error
functional

${\mathcal{J}}_{1}=\frac{1}{2}\int\limits_{0}^{t_{f}}\left[\Delta\eta(t)-\Delta\overline{\eta}(t)\right]^{2}\,dt,$
(10)

where $\Delta\eta(t)=\eta_{c}(t)-\eta_{a}(t){+R_{el}I(t)}$ and
$\Delta\overline{\eta}(t)=\overline{V}_{cell}(t)-\left[{U_{\rm{eq}}^{c}}(c_{{\rm{s}},c}\big{|}_{r=R_{c}})-{U_{\rm{eq}}^{a}}(c_{{\rm{s}},a}\big{|}_{r=R_{a}})\right]$
in which $\eta_{a}(t)$ and $\eta_{c}(t)$ are obtained by inverting relation
(1) defined on the anode and the cathode (each of these relations depends on
the parameters $\alpha_{\rm{a},m}$ and $k_{0}$), whereas $c_{{\rm{s}},a}$ and
$c_{{\rm{s}},c}$ are the solutions of system (4)–(8) obtained at the two
electrodes (additional information how ${\mathcal{J}}_{1}$ evaluated is
provided in the next section). The error functional (10) is a least-square
measure of the discrepancy between the portion of the total potential directly
depending on the constitutive relation describing ion intercalation, and is
therefore numerically better conditioned than a least-square measure of the
difference $V_{cell}(t)-\overline{V}_{cell}(t)$. This leads to the finite-
dimensional optimization problem

$\textrm{P1:}\qquad[\widehat{k}_{0,-},{\widehat{\alpha}_{a,-}},\widehat{k}_{0,+},{\widehat{\alpha}_{a,+}},\widehat{R}_{el}]=\underset{[k_{0,-},{\alpha_{a,-}},k_{0,+},{\alpha_{a,+}},R_{el}]\in\mathbb{R}^{5}}{\textrm{argmin}}~{}\mathcal{J}_{1}(k_{0,-},{\alpha_{a,-}},k_{0,+},{\alpha_{a,+}},R_{el}),$

where hats ($\widehat{\cdot}$) denote optimal values of the material
parameters and subscripts $-$ and $+$ on $k_{0}$ and $\alpha_{{a},m}$ denote
quantities defined on the anode and the cathode, respectively.

In contrast to Problem P1 where we seek optimal parameters in relations (2)
defined both on the anode and the cathode, in our main inverse problem,
Problem P2, we will focus on inferring the optimal form of the factor
$\iota(c_{{\rm{s}},a})$ in the constitutive relation (3) on the anode only
while assuming that on the cathode the constitutive relation has the standard
form (2) with parameters $\widehat{k}_{0,-}$ and ${\widehat{\alpha}_{a,-}}$
obtained by solving Problem P1. Thus, the second error functional is defined
to measure the error in the overpotential on the anode as

$\displaystyle{\mathcal{J}}_{2}$
$\displaystyle=\frac{1}{2}\int\limits_{0}^{t_{f}}\left[\eta_{a}(t)-\overline{\eta}^{P2}(t)\right]^{2}\,dt,\qquad\textrm{where}$
(11a) $\displaystyle\eta_{a}$
$\displaystyle=\psi^{-1}\left(\frac{I}{F\iota(c_{{\rm{s}},a})}\right),$ (11b)
$\displaystyle\overline{\eta}^{P2}(t)$
$\displaystyle=-\left[\overline{V}_{cell}(t)-\left({U_{\rm{eq}}^{c}}(c_{{\rm{s}},c}\big{|}_{r=R_{c}})-{U_{\rm{eq}}^{a}}(c_{{\rm{s}},a}\big{|}_{r=R_{a}})+\eta^{P1}_{c}+\widehat{R}_{el}I\right)\right],$
(11c)

in which the superscript “P1” denotes quantities computed using information
from the solution of Problem P1 and, with a slight abuse of notation,
$\psi^{-1}(\cdot)$ denotes the function
$\Delta\phi=\Delta\phi(c_{{\rm{s}},i})$ which is the solution of the equation
$I/(F\iota(c_{{\rm{s}},a}))-\psi(\Delta\phi(c_{{\rm{s}},a}),c_{{\rm{s}},a})=0$
for $c_{{\rm{s}},a}\in[0,c_{{\rm{s}},a,{\rm{max}}}]$, cf. (3) (the existence
of such a function is guaranteed by the implicit function theorem). This
notation emphasizes the dependence of the expression in (11b) on the factor
$\iota(c_{{\rm{s}},a})$ in the constitutive relation (3) describing ion
intercalation in the anode. As stipulated by assumptions A1–A3 in Section 2,
see also Supporting Information (Appendix A), we must require the function
$\iota(c_{{\rm{s}},a})$ to vanish at both vanishing and maximum
concentrations, i.e., for $c_{{\rm{s}},a}=0$ and for
$c_{{\rm{s}},a}=c_{{\rm{s}},a,{\rm{max}}}$. We also need the function
$\iota(c_{{\rm{s}},a})$ to be continuous with square-integrable
(distributional) derivatives, such that it will belong to the Sobolev space51
$H_{0}^{1}([0,c_{{\rm{s}},a,{\rm{max}}}])$ (the subscript “0” indicates that
functions belonging to this space vanish at the boundary). The optimal
functional form of the exchange current $\iota(c_{{\rm{s}},a})$ on the anode
can then be obtained as a solution of the optimization problem

$\textrm{P2:}\qquad[\widehat{\iota}(c_{{\rm{s}},a})]=\underset{\iota(c_{{\rm{s}},a})\in{H}_{0}^{1}([0,c_{{\rm{s}},a,{\rm{max}}}])}{\textrm{argmin}}~{}\mathcal{J}_{2}(\iota(c_{{\rm{s}},a})),$

where the corresponding constitutive relation (3) takes the optimal form
$\widehat{j}(c_{{\rm{s}},a})=\widehat{\iota}(c_{{\rm{s}},a})\psi(\Delta\phi,c_{{\rm{s}},a})$.
We note that the framework outlined above could also be utilized to
reconstruct the optimal form of the exchange current on the cathode. In this
case, the cost functional (11) would include $\eta_{c}(t)$, equation (11b)
would involve $\eta_{c}$ as a function of $\iota(c_{{\rm{s}},c})$ whereas
equation (11c) would use $\eta_{a}^{P1}$ (with a correction to the sign on the
right-hand side, cf. equation (9)). In such case the optimal form of the
constitutive relation would be
$\widehat{j}(c_{{\rm{s}},c})=\widehat{\iota}(c_{{\rm{s}},c})\psi(\Delta\phi,c_{{\rm{s}},c})$,
cf. (3). Moreover, one could solve these two problems iteratively making it
possible to refine the reconstructions of the exchange current on both
electrodes. Our approach to solving Problems P1 and P2 is described next.

###  4.3 Computational Approach

Evaluation of the error functionals (10) and (11) requires solution of the
SPM, cf. (4)–(9), which is done using a standard finite element method 10, 32.
Problem P1 represents unconstrained minimization over $\mathbb{R}^{5}$ and is
solved numerically in a straightforward manner using the function fminunc in
MATLAB. The key step is evaluation of the error functional (10) and the
procedure is outlined in Algorithm 1.

Algorithm 1 Evaluation of the error functional ${\mathcal{J}}_{1}$ in the
solution of Problem P1
Input:
$\overline{V}_{\textrm{cell}}(t)$ — experimental voltage as a function of time
$k^{o}_{0,-}$, $k^{o}_{0,+}$ — initial guesses for reaction rate for anode and
cathode
$\alpha^{o}_{a,-}$, $\alpha^{o}_{a,+}$ — initial guesses for transfer
coefficients for anode and cathode
$R^{o}_{el}$ — initial guess for electrolyte resistance
Output:
${\mathcal{J}}_{1}({k}_{0,-},{k}_{0,+},{\alpha}_{a,-},{\alpha}_{a,+},{R}_{el})$

Solve the SPM equations (4)-(8) for $c_{{\rm{s}},a}$ and $c_{{\rm{s}},c}$

Compute
$\Delta\bar{\eta}(t)=-\left[\overline{V}_{cell}(t)-\left({U_{\rm{eq}}^{c}}(c_{{\rm{s}},c}\big{|}_{r=R_{c}})-{U_{\rm{eq}}^{a}}(c_{{\rm{s}},a}\big{|}_{r=R_{a}})\right)\right]$

Compute ${\eta}_{a}$ and ${\eta}_{c}$ by inverting relations 1 and 8 with
unknowns ${k}_{0,-}$, ${k}_{0,+}$,${\alpha}_{a,-}$, and ${\alpha}_{a,+}$

Compute $\Delta{\eta}(t)={\eta}_{c}-{\eta}_{a}+{R}_{el}I(t)$

Determine
${\mathcal{J}}_{1}({k}_{0,-},{k}_{0,+},{\alpha}_{a,-},{\alpha}_{a,+},{R}_{el})$
by evaluating the integral in (10)

Problem P2 represents minimization over an infinite-dimensional function space
$H^{1}_{0}$ and therefore requires a specialized approach. We compute the
optimal exchange current using the discrete gradient flow as
$\widehat{\iota}(c_{{\rm{s}},i})=\lim_{n\rightarrow\infty}\iota^{(n)}(c_{{\rm{s}},i})$
where the approximations $\iota^{(n)}(c_{{\rm{s}},i})$ are computed
iteratively by

$\displaystyle\iota^{(n+1)}(c_{{\rm{s}},i})$
$\displaystyle=\iota^{(n)}(c_{{\rm{s}},i})-\tau^{(n)}\nabla\mathcal{J}_{2}(\iota^{(n)}(c_{{\rm{s}},i})),$
(12) $\displaystyle\iota^{(0)}(c_{{\rm{s}},i})$
$\displaystyle=\iota_{0}(c_{{\rm{s}},i}),$

in which $\iota_{0}(c_{{\rm{s}},i})$ is a suitable initial guess,
$\nabla\mathcal{J}_{2}(\iota^{(n)}(c_{{\rm{s}},i}))$ is the Sobolev gradient
of the error functional (11) at the $n$th iteration, and $\tau^{(n)}$ is the
step size along the gradient direction. The initial guess
$\iota_{0}(c_{{\rm{s}},i})$ is chosen to vanish at $c_{{\rm{s}},i}=0$ and at
$c_{{\rm{s}},i}=c_{{\rm{s}},i,{\rm{max}}}$, so that the same limiting
behaviour will be inherited by the optimal solution
$\widehat{\iota}(c_{{\rm{s}},i})$, provided the gradients
$\nabla{\mathcal{J}}_{2}(c_{{\rm{s}},i})$ also vanish at $c_{{\rm{s}},i}=0$
and $c_{{\rm{s}},i}=c_{{\rm{s}},i,{\rm{max}}}$, cf. Supporting Information
(Appendix A). In practice, the initial guess $\iota_{0}(c_{{\rm{s}},i})$ is
taken as the factor $i_{0}$ in the Butler-Volmer relation (2b) with optimal
parameters found by solving Problem P1 which, by construction, ensures the
correct behaviour at $c_{{\rm{s}},i}=c_{0},c_{{\rm{s}},i,{\rm{max}}}$.
Computation of the Sobolev gradient
$\nabla\mathcal{J}_{2}(\iota(c_{{\rm{s}},i}))$, which is a key element of
approach (12), is described in Supporting Information (Appendix B) and
together with the properties of the initial guess $\iota_{0}(c_{{\rm{s}},i})$
ensures the required smoothness of the optimal exchange current
$\widehat{\iota}(c_{{\rm{s}},i})$ and its correct behaviour at
$c_{{\rm{s}},i}=c_{0},c_{{\rm{s}},i,{\rm{max}}}$. The optimal step size in
(12) is found by solving the line-minimization problem

$\tau^{(n)}=\underset{\tau>0}{\operatorname{argmin}}~{}\mathcal{J}_{2}(\iota(c_{{\rm{s}},i})-\tau\nabla^{{H}^{1}}\mathcal{J}_{2}(\iota(c_{{\rm{s}},i}))),$
(13)

which is done conveniently using Brent’s algorithm52. Expression (12) is
evaluated using the collocation approach on the domain
$[0,c_{{\rm{s}},i,{\rm{max}}}]$ discretized with a uniform grid based on
$N=500$ grid points. We also considered finer meshes, but they did not produce
appreciable differences in the obtained results, hence the aforementioned
resolution was used as a trade-off between accuracy and the computational
cost. When determining the Sobolev gradients, one must specify the value of
the Sobolev parameter $0<\ell<\infty$ which appears in the definition of the
inner product (24) and determines the relative smoothness of the gradients.
More specifically, $\ell$ may be regarded as a characteristic, or cut-off,
value of the concentration below which all information in the gradient
$\nabla\mathcal{J}_{2}(\iota(c_{{\rm{s}},i}))$ is smoothed out. At the
validation stage, values of $\ell$ spanning several orders of magnitude were
considered to determine the value resulting in the fastest rate of convergence
of iterations (12). Based on these tests, we chose $\ell=3000000$ and
$\ell=100000$ for the case of slow charge and moderate discharge rates,
respectively. The steps leading to the solution of Problem P2 are summarized
as Algorithm 2.

The computational approach was systematically validated by solving Problems P1
and P2 with “manufactured” measurement data. More specifically, we assumed a
certain form of the constitutive relation ((2) with some values of
$\alpha_{{a},m}$ and $k_{0}$ for Problem P1 and (3) with some function
$\iota(c_{{\rm{s}},i})$ for Problem P2) which was used to generate
measurements. The true constitutive relations were then reconstructed using
Algorithm 2. Excellent reconstruction accuracy was achieved with errors
essentially at the level of errors due to numerical discretization. We add
that Problems P1 and P2 are non-convex and we have found some evidence for the
presence of nonunique (local) minimizers when different initial guesses are
used in Algorithms 1 and 2. However, in all instances these local minimizers
involved physically inconsistent material properties (e.g., with incorrect
signs), which made it easy to eliminate them. Analysis of the sensitivity of
solutions to Problems P1 and P2 to perturbations of measurement data is
discussed in the next section.

Algorithm 2 Solution of Problem P2 using discrete gradient flow (12)
Input:
$\epsilon$ — tolerance
$\ell$ — Sobolev parameter
$\overline{V}_{\textrm{cell}}(t)$ — experimental voltage as a function of time
$k^{o}_{0,-}$, $k^{o}_{0,+}$ — initial guesses for reaction rate for anode and
cathode
$\alpha^{o}_{a,-}$, $\alpha^{o}_{a,+}$ — initial guesses for transfer
coefficients for anode and cathode
$R^{o}_{el}$ — initial guess for electrolyte resistance
Output: $\widehat{k}_{0,-}$, $\widehat{k}_{0,+}$,$\widehat{\alpha}_{a,-}$,
$\widehat{\alpha}_{a,+}$, $\widehat{R}_{el}$, and
$\widehat{\iota}(c_{{\rm{s}},i})$

Solve Problem P1 using fminunc

Compute
$\overline{\eta}(t)\leftarrow-\left[\overline{V}_{cell}(t)-\left({U_{\rm{eq}}^{c}}(c_{{\rm{s}},c}\big{|}_{r=R_{c}})-{U_{\rm{eq}}^{a}}(c_{{\rm{s}},a}\big{|}_{r=R_{a}})+\eta^{P1}_{c}+\widehat{R}_{el}I\right)\right]$

Set
$\iota_{0}(c_{{\rm{s}},i})=F\widehat{k}_{0,-}c_{\rm{e}}^{\widehat{\alpha}_{a,-}}\left(c_{{\rm{s}},i,{\rm{max}}}-c_{{\rm{s}},i}\right)^{\widehat{\alpha}_{a,-}}c_{{\rm{s}},i}^{(1-\widehat{\alpha}_{a,-})}$

repeat

Evaluate $\nabla^{L^{2}}\mathcal{J}(\iota(c_{{\rm{s}},i})))$

Compute $\nabla^{H^{1}}\mathcal{J}(\iota(c_{{\rm{s}},i})))$ given
$\nabla^{L^{2}}\mathcal{J}(\iota(c_{{\rm{s}},i})))$

Perform a Polak-Ribiere conjugate gradient update
$\nabla^{H^{1}}\mathcal{J}(\iota(c_{{\rm{s}},i})))$

Find the step size $\tau^{(n)}$ by solving the line-search problem (13) with
Brent’s method

$\iota^{(n+1)}(c_{{\rm{s}},i})\leftarrow\iota^{(n)}(c_{{\rm{s}},i})-\tau^{(n)}\nabla\mathcal{J}(\iota^{(n)}(c_{{\rm{s}},i}))$

until
$\frac{|\mathcal{J}(\iota^{(n)}(c_{{\rm{s}},i})))|-|\mathcal{J}(\iota^{(n-1)}(c_{{\rm{s}},i})))|}{|\mathcal{J}(\iota^{(n)}(c_{{\rm{s}},i})))|}<\epsilon$

###  4.4 Sensitivity Analysis

Inverse problems are generally ill-posed and variations in the experimental
data can result in significant shifts in the reconstructed solutions 19. To
explore the effects of small perturbations in the experimental data on the
solution of Problems P1 and P2, a Monte-Carlo analysis is developed where the
measurement data is perturbed as

$\widehat{V}_{cell}(t)=V_{cell}(t)+\sum\limits_{k=1}^{3}\frac{a_{k}}{k^{2}}\sin\left(\frac{2\pi
kt}{t_{f}}\right)+\sum\limits_{k=1}^{3}\frac{b_{k}}{k^{2}}\cos\left(\frac{2\pi
kt}{t_{f}}\right),\qquad t\in[0,t_{f}],$ (14)

in which $a_{k}$ and $b_{k}$, $k=1,2,3$, are normally distributed random
variables with zero mean and a prescribed standard deviation. Problems P1 and
P2 with measurement data perturbed as in (14) are then solved $M$ times ($M\gg
1$), each time using a different sample of the random variables $a_{k}$ and
$b_{k}$. This allows us to compute the statistics of the reconstruction errors
and of the optimal constitutive relations $\widehat{\iota}(c_{{\rm{s}},a})$
corresponding to perturbed measurements which provide information about their
sensitivity to noise.

Figure 3: PDFs of the error functionals (a) ${\mathcal{J}}_{1}$ and (b)
${\mathcal{J}}_{2}$ obtained in Monte-Carlo sensitivity analysis. The vertical
lines represent the values of ${\mathcal{J}}_{1}$ and ${\mathcal{J}}_{2}$
obtained in solutions of Problems P1 and P2 with no noise.

Our sensitivity analysis is then performed based on the manufactured
measurement data at 1C (which was also used to validate Algorithm 2) and using
$M=3000$ normally distributed noise samples with zero mean and $0.0025$
standard deviation. Using the empirical $3$-sigma rule53, the relative
variation of the measurement data due to noise can be estimated as $6\times
0.0025/1.3\approx 0.012$ or $1.2\%$, where $1.3$ in the denominator is the
approximate total voltage drop during the (dis)charge cycle. The effect of
this noise is thus comparable to the potential drop in the electrolyte, cf.
Tables 1 and 2, which is accounted for in the voltage (9) predicted by the
augmented SPM. The probability density functions (PDFs) of the error
functionals (10) and (11) obtained by solving Problems P1 and P2 with
measurement data perturbed as described above are shown in Figure 3. The
corresponding mean inferred exchange currents are shown in Figure 4 together
with their 95% confidence intervals and the standard Butler-Volmer relation
(1) with some reference parameter values that were used to “manufacture” the
measurements. We observe that optimal reconstructions obtained by solving both
Problem P1 and P2 capture the general functional form of the true constitutive
relation quite well, but the magnitude is off by a bit more than 10% even
though the noise amplitude is rather small. It ought to be emphasized however
that such level of uncertainty is not uncommon in parametrization of
electrochemical systems9, 10. The corresponding distributions of optimal anode
overpotentials are shown in Figure 5 where deviations from the original data
are rather small. We also mention that additional Monte-Carlo simulations were
performed using the random coefficients $a_{k}$ and $b_{k}$ in (14) with
several decreasing standard deviations and in this limit of vanishing noise
amplitude, the deterministic solution is recovered, confirming the consistency
of the calculations. These observations allow us to conclude that a modest
amount of noise in the measurement data has some quantitative and little
qualitative effect on the obtained reconstructions.

Figure 4: Optimal reconstructions of the exchange current (black dashed lines)
$i_{0}(c_{\rm{s}})$ and (blue circles) $\widehat{\iota}(c_{{\rm{s}},a})$
obtained by solving Problems P1 and P2 with measurement data perturbed as in
(14) for (a) the entire concentration range and (b) intermediate
concentrations indicated by the dotted box in (a). The shaded bands indicate
the $95\%$ confidence intervals whereas the vertical lines represent the
bounds $c_{<}$ and $c_{>}$ of the identifiability interval ${\mathcal{I}}$.
The true exchange current used to manufacture the measurements is shown in
green.

Figure 5: Anode overpotentials corresponding to the optimal exchange currents
shown in Figure 4 for (a) the entire normalized time window and (b)
magnification around intermediate times indicated by the dotted box in (a).
The line types and colour-coding are the same as in Figure 4 with an added red
solid line to represent the average value of the overpotential when the total
potential is perturbed with noise.

##  5 Results

We compare the optimal constitutive relations reconstructed by solving
Problems P1 and P2 using data for slow charge and moderate discharge rates
against the standard form of the Butler-Volmer relation (1) for a graphite
electrode. In the latter, we set $\alpha_{{a},m}=\alpha_{{c},m}=1/2$, which is
a standard choice for the transfer coefficients13, 26, and in the problem with
slow charge rates we set $k_{0}=1.8742\times 10^{-12}$ which is taken from a
similar set-up of a graphite anode with an NMC$622$ cathode54, whereas in the
problem with moderate charge rates we set $k_{0}=6.72\times 10^{-12}$ which
was measured experimentally when collecting the data for a graphite anode with
an NMC$811$ cathode 30. Hereafter we will refer to these values of
$\alpha_{{a},m}$, $\alpha_{{c},m}$ and $k_{0}$ as “nominal”.

###  5.1 Slow Charge Rates

We first consider solutions of Problem P1 where scalar parameters are
calibrated in expression (2b) before analyzing the results obtained by solving
Problem P2 in which the exchange current $\iota(c_{{\rm{s}},a})$ is sought in
a fairly general form. In Figure 6, we show the values of the error functional
(10) corresponding to solutions of Problem P1 and see that compared to the
values of ${\mathcal{J}}_{1}$ obtained using the standard Butler-Volmer
relation (1), they are reduced only mildly, by less than one order of
magnitude for both C-rates. The solution to Problem P1 suggests that a small
deviation from the standard functional form of the exchange current density
already provides a relatively good fit.

Figure 6: Cost functional ${\mathcal{J}}_{2}$ as function of iterations in the solution of Problem P2 using Algorithm 2 at the C-rate of C$/10$ (blue lines with circles) and C$/3$ (black lines with circles). The horizontal lines represent the values of the error functional ${\mathcal{J}}_{1}$ corresponding to the standard Butler-Volmer equation (1) (dotted lines) and the terminal values obtained in the solutions of Problem P1 at C-rates of C$/10$ (blue solid lines) and C$/3$ (black solid lines).  Charge Rate | Parameter | Value
---|---|---
C$/10$ | anode reaction rate constant, $\widehat{k}_{0,-}$ | $1.076\times 10^{-12}$
| cathode reaction rate constant, $\widehat{k}_{0,+}$ | $4.416\times 10^{-12}$
| anode anodic transfer coefficient, ${\widehat{\alpha}_{a,-}}$ | $0.5667$
| cathodic anodic transfer coefficient, ${\widehat{\alpha}_{a,+}}$ | $0.4948$
| electrolyte resistance, $\widehat{R}_{el}$ | $0.5813$
C$/3$ | anode reaction rate constant, $\widehat{k}_{0,-}$ | $2.428\times 10^{-12}$
| cathode reaction rate constant, $\widehat{k}_{0,+}$ | $4.93\times 10^{-12}$
| anode anodic transfer coefficient, ${\widehat{\alpha}_{a,-}}$ | $0.4999$
| cathodic anodic transfer coefficient, ${\widehat{\alpha}_{a,+}}$ | $0.4959$
| electrolyte resistance, $\widehat{R}_{el}$ | $0.5625$
Table 3: Material parameters obtained by solving Problem P1.

The material parameters inferred by solving Problem P1 are summarized in Table
3. We see that the optimal parameter values obtained at the two C-rates are
quite similar and from Figure 7, we see that the dependence on
$c_{{\rm{s}},a}$ is similar to the standard Butler-Volmer relation, but with a
reduction in magnitude due to values of $\widehat{k}_{0}$ being smaller by two
orders of magnitude.

Figure 7: Optimal reconstructions of the exchange current obtained by solving
Problem P1 (dashed lines) and Problem P2 (solid lines) as function of the
normalized concentration $c_{\rm{s}}/c_{\rm{s,max}}$ for (a) the entire
concentration range and (b) small values of the exchange current. Panel (a)
also shows the standard Butler-Volmer relation (1) with the nominal parameter
values (dotted lines), whereas the vertical lines represent the bounds of the
identifiability intervals at the two C-rates.

The reconstructed overpotentials and the total cell potentials are compared as
functions of time to the experimental data in Figure 8. While the anode
overpotential corresponding to the Butler-Volmer relation with parameters
calibrated by solving Problem P1, cf. Table 3, has values coinciding with a
subset of the values observed in the actual measurements, its dependence on
time $t$ is incorrect. The same behaviour is also evident, though less
pronounced due to the presence of large equilibrium potentials, in Figure 8b
where we show the total potentials.

We now move on to discuss the results obtained by solving Problem P2. In
Figure 6, we see that after a few hundred iterations in Algorithm 2, the value
of the error functional ${\mathcal{J}}_{2}$ is reduced by several orders of
magnitude. The obtained optimal forms $\widehat{\iota}(c_{{\rm{s}},a})$ of the
exchange current shown in Figure 7 are strikingly similar for both C-rates.
The optimal solutions are skewed towards smaller concentrations where larger
exchange currents are predicted, unlike solutions of Problem P1. In addition,
they feature localized peaks at concentration values near the smallest
resolved values. As was the case with solutions of Problem P1, the optimal
functions $\widehat{\iota}(c_{{\rm{s}},a})$ inferred at the two C-rates are
rather similar, except that the peak at low concentrations is steeper at the
higher C-rate. Figure 8 shows that, as expected from the small errors evident
in Figure 6, both the overpotentials and the total potentials predicted by the
constitutive relation (3) with the optimal exchange currents
$\widehat{\iota}(c_{{\rm{s}},a})$ shown in Figure 8 provide a better fit to
the measured data.

Figure 8: (a) Anode overpotentials and (b) total potentials predicted by the
constitutive relation (3) using optimal reconstructions of the exchange
current obtained by solving Problem P1 (dot-dashed lines) and Problem P2
(dashed lines) as function of the normalized time $t/t_{tot}$. Solid lines and
dotted line represent the experimental measurements and predictions obtained
using the standard Butler-Volmer relation (1) with the nominal parameter
values, respectively.

###  5.2 Moderate Discharge Rates

Overall, the results obtained by solving Problems P1 and P2 using the data for
moderate discharge rates are qualitatively similar to the results presented in
the previous section. In our analysis below, we will therefore focus on
highlighting the differences.

Charge Rate | Parameter | Value
---|---|---
C/2, discharge | anode reaction rate constant, $k_{0,-}$ | $3.165\times 10^{-12}$
| cathode reaction rate constant, $k_{0,+}$ | $1.138\times 10^{-12}$
| anode anodic transfer coefficient, ${\alpha_{a,-}}$ | $0.5001$
| cathodic anodic Transfer Coefficient, $\alpha_{c,a}$ | $0.4999$
| electrolyte resistance, $R_{el}$ | $0.0010$
1C, discharge | anode reaction rate constant, $k_{0,-}$ | $2.414\times 10^{-12}$
| cathode reaction rate constant, $k_{0,+}$ | $0.867\times 10^{-12}$
| anode anodic transfer coefficient, ${\alpha_{a,-}}$ | $0.5001$
| cathodic anodic transfer coefficient, $\alpha_{c,a}$ | $0.4999$
| electrolyte resistance, $R_{el}$ | $0.0010$
Table 4: Material parameters obtained by solving Problem P1.

The optimal parameters inferred by solving Problem P1 and collected in Table 4
are in the same range as those obtained at low charging rates, cf. Table 3,
and also show little variability between the two C-rates considered. The
values of the error functional (10) shown in Figure 9 reveal a similar
improvement resulting from solving Problem P1 as for the slow discharge rates,
cf. Figure 6. However, while in solutions of Problem P2 the reduction of the
least-square error ${\mathcal{J}}_{2}$ is still substantial, the errors
saturate at the level of $\mathcal{O}(10^{-2}-10^{-3})$.

The optimal form of expression (2b) inferred by solving Problem P1 and shown
in Figure 10 is again similar to the nominal form, but has a smaller
amplitude, roughly by a factor of 2\. On the other hand, the optimal forms of
the exchange current $\widehat{\iota}(c_{{\rm{s}},a})$ inferred by solving
Problem P2 reveal significant differences with respect to the corresponding
functions obtained at low charging rates, cf. Figure 7. Most importantly, in
the present case large values of the exchange current
$\widehat{\iota}(c_{{\rm{s}},a})$ occur at large concentrations close to
$c_{{\rm{s}},a,{\rm{max}}}$. The data in Figure 11 confirm the small errors
evident in Figure 9 are achieved with total potentials closely matching the
measurement data. We note that by reducing the value of parameter $\ell$ in
(24) it is in fact possible to obtain overpotentials and total potentials that
are indistinguishable (within graphical resolution) from the measured data.
However, the price to be paid for that would be increasingly irregular forms
of the optimal exchange current $\widehat{\iota}(c_{{\rm{s}},a})$.

Figure 9: Cost functional ${\mathcal{J}}_{2}$ as function of iterations in the
solution of Problem P2 using Algorithm 2 at the C-rate of C/2 (blue lines with
circles) and 1C (black lines with circles). The horizontal lines represent the
values of the error functional ${\mathcal{J}}_{1}$ corresponding to the
standard Butler-Volmer equation (1) (dotted lines) and the terminal values
obtained in the solutions of Problem P1 at C-rates of C/2 (blue solid lines)
and 1C (black solid lines).  Figure 10: Optimal reconstructions of the
exchange current obtained by solving Problem P1 (dashed lines) and Problem P2
(solid lines) as function of the normalized concentration
$c_{\rm{s}}/c_{\rm{s,max}}$ for the entire concentration range. The standard
Butler-Volmer relation (1) with the nominal parameter values is shown using a
dotted line whereas the vertical lines represent the bounds of the
identifiability intervals at the two C-rates.

Figure 11: (a) Anode overpotentials and (b) total potentials predicted by the
constitutive relation (3) using optimal reconstructions of the exchange
current obtained by solving Problem P1 (dot-dashed lines) and Problem P2
(dashed lines) as function of the normalized time $t/t_{tot}$. Solid lines and
dotted line represent the experimental measurements and predictions obtained
using the standard Butler-Volmer relation (1) with the nominal parameter
values, respectively.  Figure 12: Comparison of optimal reconstructions
$\widehat{\iota}(c_{\rm{s}})$ of the exchange current obtained by solving
Problem P2 using all four datasets considered.

Finally, we attempt to extract some universal features of the optimal exchange
currents determined using measurements acquired at the different
charge/discharge rates and, to this end, in Figure 12 we compare the optimal
reconstructions $\widehat{\iota}(c_{{\rm{s}},a})$ obtained by solving Problem
P2 for all considered datasets. We see that there is a range of concentrations
$0.25\times 10^{4}\lessapprox c_{{\rm{s}},a}\lessapprox 2\times 10^{4}$ where
the inferred exchange currents are qualitatively similar revealing some
possibly universal features which include a steep increase of
$\widehat{\iota}(c_{{\rm{s}},a})$ at lower concentrations $c_{{\rm{s}},a}$
followed by a gentle and nonmonotonic decrease at higher concentrations.
Outside that range the optimal reconstructions
$\widehat{\iota}(c_{{\rm{s}},a})$ exhibit distinct behaviour specific to each
dataset.

##  6 Discussion and Concluding Remarks

In this work, we have developed an inverse modelling approach to infer an
optimal form of the concentration-dependent exchange current
$\widehat{\iota}(c_{{\rm{s}},a})$ appearing in a generalization (3) of the
Butler-Volmer relation describing intercalation of Li ions into an active
material. While many inverse problems have already been solved in the area of
electrochemistry21, 6, 22, 23, 9, 10, this is to the best of our knowledge the
first attempt to infer constitutive relations describing an interfacial
phenomenon. Although the proposed approach relies on many simplifications,
most importantly, on the use of the SPM to describe the transport of Li, it is
shown to be computationally robust and thermodynamically consistent. In
particular, the model equipped with the inferred optimal exchange currents
$\widehat{\iota}(c_{{\rm{s}},a})$ predicts very accurately the time evolution
of the cell voltage which strongly suggests that ansatz (3) and assumptions
A$1$–A$3$ are physically consistent (we note that such a good match with
experimental data need not occur if the model used in the inverse problem does
not accurately describe the physics of the problem9). Interestingly, by
exploiting the coupling between different scales in the model, our approach
allows us to reconstruct a constitutive relation defined on the microscale
based on macroscale measurements. The variational formulation of the inverse
problem involves averaging of the reconstruction errors over time and
concentrations such that the vanishing of the exchange current for small and
large concentrations, cf. Assumptions A$2$ and A$3$, does not cause problems.

We emphasize that the observations made below are predicated on the accuracy
of the material properties which were assumed known in the solution of
Problems P1 and P2, cf. Supporting Information (Appendix C). While a
significant effort was made to determine these parameters accurately, their
inherent uncertainties may have some effect on the obtained results. In
particular, one material property known to be difficult to determine with high
accuracy is the diffusivity $D_{\rm{s}}(c_{\rm{s}})$ appearing in (4) 55. To
probe the effect of its uncertainty, we have solved Problems P1 and P2 in
which this diffusivity was allowed to vary over several orders of magnitude.
While quantitative details of the reconstructed constitutive relations were
different, the key qualitative features of the optimal exchange currents as
discussed below were essentially unchanged. That being said, the present study
represents only a proof of the concept and more work needs to be done to
understand the extent to which our observations about optimal forms of the
exchange current could generalize to other battery materials.

The optimal exchange currents $\widehat{\iota}(c_{{\rm{s}},a})$ reconstructed
by solving Problem P2 based on datasets corresponding to different
charging/discharging rates show a qualitatively similar dependence on
$c_{{\rm{s}},a}$ for intermediate concentrations $0.25\times 10^{4}\lessapprox
c_{{\rm{s}},a}\lessapprox 2\times 10^{4}$, cf. Figure 12. This behaviour is
different from the form of the exchange current predicted by the standard
Butler-Volmer relation (2b), even after its parameters are calibrated by
solving Problem P1. Outside that range, the optimal exchange currents exhibit
non-universal features specific to individual cases. In particular, the
optimal exchange currents inferred from data representing charging and
discharging feature sharp peaks at, respectively, low and high concentrations,
cf. Figures 7 and 10, indicating possible hysteresis effects in the
constitutive relations. In addition, the results shown in Figure 12 indicate
that the exchange current density tends to be higher during discharging which
is consistent with the observations made in earlier studies 56 where this
effect was partly attributed to a higher de-solvation energy barrier.

The sharp peaks characterizing the optimal exchange current at low
concentrations when charging and at high concentrations when discharging
suggest that the electrode kinetics are faster at the beginning of the cycle.
This might be due to some interfacial phenomena not accounted for in our model
such as, e.g., capacitance effects and variation in the activity
coefficient57. It is known that the structure of the double layer and ion
adsorption can affect electrode kinetics26. These phenomena have been
implicated in causing anomalies in rate constants and the exchange current, as
well as variations in current-potential plots26. The behaviour observed at
early times in the charge and discharge cycles could be compensating for
residual effects produced by the double layer at the electrode interface. In
this context, modelling efforts have produced corrections to the Butler-Volmer
equation, such as the Frumkin correction which allows the rate constant to be
dependent on potential, in order to take into account the effects of the
electric double layer and ion adsorption in electrode kinetics26. On the other
hand, a thermodynamically consistent derivation of the exchange current
density shows that its functional form is dependent on the rate constant,
concentration, and the activity coefficient 57. At the same time, the activity
of the relevant species is often assumed to be $1$, which may not be
reflective of the dynamics in the system and could lead to its non-ideal
behaviour26.

Another possible reason for the presence of large exchange currents at low and
high concentrations could be due to the system moving away from equilibrium at
the beginning of the cycle protocol when there is a uniform distribution of Li
ions in the electrolyte throughout the cell26, 13. When charging begins, there
is a richness of Li near the particle and as the protocol continues, Li moves
in the electrolyte towards the anode. The spike observed at early times of the
charging protocol could possibly be attributed to the richness of Li near the
electrode, which could contribute to a larger exchange current and,
consequently, a larger current density at the interface. In the case of a
discharge protocol, there is a richness of Li ions in the anode near the
surface and as the cycle continues, Li ions must diffuse from deeper in the
electrode to the surface. This is not accounted for in the SPM because the
concentration $c_{\rm{e}}$ in the electrolyte is assumed constant, cf. Section
4.1.

The standard formulation of the Butler-Volmer model is based on a simple,
single step reaction and as such does not allow for formation of short-lived
intermediate products. This is a key assumption which may not be valid in
experimental settings; for example, Li-ion batteries often suffer from Li-ion
plating, which is a result of fast charging and Li metal plating out onto the
surface of the electrode58. While this plated material may absorb into the
electrode after a certain relaxation time, the side reaction is not captured
by the Butler-Volmer kinetic relation 59, 60, 61. Some modelling efforts have
attempted to extend the Butler-Volmer to account for such effects62, 63, 64,
65.

It is also possible that with peaks at low and high Li concentrations the
optimal exchange currents $\widehat{\iota}(c_{{\rm{s}},a})$ may be
compensating for certain effects other than interfacial phenomena which are
not accounted for in the model. One such possibility has to do with the nature
of Li transport in the active material, which under certain circumstances is
known to lead to phase transformations and in turn can significantly affect
the diffusivity66. The spike in the reconstructed exchange current is located
at approximately $12$–$15\%$ capacity, cf. Figure 7, which could be a result
of a rapid increase in the diffusion directly after the phase transition in
graphite. When comparing the plots of ${U_{\rm{eq}}^{a}}$ and
$\frac{d{U_{\rm{eq}}^{a}}}{d(c_{{\rm{s}},a}/c_{{\rm{s}},a,{\rm{max}}})}$ vs.
$c_{{\rm{s}},a}/c_{{\rm{s}},a,{\rm{max}}}$ (the latter not shown for brevity),
there does appear to be agreement between the minima/maxima observed in the
optimal exchange current and the remnants of these phase transitions captured
by the equilibrium potential, ${U_{\rm{eq}}^{a}}$. The reconstructed
constitutive relations could be therefore compensating for an inadequate
transport model used for graphite. It is known that diffusion models are
difficult to quantify experimentally and that the physics of transport in the
solid is complicated by interfacial phenomena9. The literature often reports
use of different models of transport in the solid phase based on Fickian
diffusion, non-linear diffusion and the Cahn-Hilliard phase separation 9.
However, the anomalous behaviour of optimal exchange currents at low and high
concentrations persisted even when a strongly non-monotonic diffusivities
$D(c_{{\rm{s}},i})$, similar to the properties measured or inferred for
graphite66, 9 were used in the transport equation (4).

More work based on additional datasets is thus needed in order to understand
the precise origins of the universal features in the optimal exchange currents
$\widehat{\iota}(c_{{\rm{s}},i})$ inferred based on different datasets. In
particular, it would be beneficial to understand what these universal features
can tell us about the validity of the different assumptions used to derive the
standard Butler-Volmer relation (1), so that these assumptions could be
suitably refined to produce more accurate constitutive models. Even though we
have given proof of concept for the approach and its application to
interfacial phenomena in electrochemistry, the question of whether it can be
used to learn a general improvement to the widely-used “standard” exchange
current density in the Butler-Volmer relation (1) remains open. Next steps
towards answering this question might include the application of the method
using the full DFN framework, thereby generalizing the physics and allowing us
to examine higher C-rates. Or, arguably better, one would like to apply this
approach to experimental data specifically designed for our purposes, e.g.,
this might involve a non-porous electrode such that there are fewer device and
model parameters, thereby spotlighting the dependence on the interfacial
constitutive relation. In the best case, it may then be possible to learn
something about the fundamental electrochemistry occurring at material
interfaces.

As is evident from the discussion in Section 4.4, a key limiting factor of our
approach is its sensitivity to noise. The main issue is that the useful
information about the interfacial reaction kinetics is encoded in the
overpotential $\eta$ which at low C-rates represents only a small fraction of
the measured cell voltage, cf. Tables 1–2. On the other hand, at higher
C-rates where this proportion increases, the SPM loses its validity. A
preferred solution would thus be to formulate the inverse problem based on the
full DFN model (or one of its reductions, but more complete than the SPM), as
this would make it possible to use data obtained at higher C-rates where the
overpotential $\eta$ represents a larger proportion of the measured voltage.
Such an inverse problem will be mathematically more complicated since, unlike
the case of the SPM, the solutions of the DFN model depend on the form of the
constitutive relation describing ion intercalation. Another limitation of the
proposed approach is that the optimal forms of the constitutive relation
depend on the C-rate, cf. Figures 10 and 12. While this is not uncommon when
an inverse problem is solved using multiple experimental datasets, this
complicates application of the method in practice when one needs to model
battery cells at different C-rates. A possible solution is to infer an
“average” optimal exchange current by simultaneously fitting voltage curves
corresponding to several C-rates. When extending the proposed analysis to
alternative materials, considerations such as multiphase intercalation67, 68
and mixed conduction effects44 would need to be taken into account. The former
would necessitate redefining diffusivity in the bulk whereas the latter would
require simultaneous solution of the diffusion and conduction problems at the
particle level. We plan to pursue such extensions of the present approach in
the near future.

In principle, the framework developed in the present investigation could be
used to study other interfacial phenomena, where experimental data can be
leveraged to learn the optimal forms of the relevant constitutive relations.
Some interesting problems of this type include the reconstruction of
constitutive relations governing the Hodgkin-Huxley dynamics of ion transfer
through a cell membrane 17, the Shockley-Read-Hall recombination of electrons
in a semiconductor 69, 70, or the charge density of an adsorbed layer formed
between a charged interface and an electrolyte 71.

LD, KJS, GRG and BP were supported by a Collaborative Research & Development
grant # CRD494074-16 from Natural Sciences & Engineering Research Council of
Canada. JF and SS were supported by the Faraday Institution, UK (grant number
FIRG003).

##  Appendix A Limiting Behaviour of Constitutive Relationships at Small and
Large Concentrations

We derive conditions that must be satisfied by the function $\iota(c_{s})$ for
$c_{s}\rightarrow 0,c_{s,max}$ in order for the constitutive relation (3) to
be physically consistent. We note that $c_{s}=\iota^{-1}(\Delta\phi)$ relates
the potential $\Delta\phi$ to the concentration $c_{s}$ and the range of
$\iota^{-1}(\Delta\phi)$ is $0\leq\iota^{-1}(\Delta\phi)\leq c_{s,max}$. The
equilibrium potential can be separated into a bounded portion $u_{eq}$ and an
unbounded part given by 28

${U_{eq}^{i}}(\Delta\phi)=\frac{RT}{F(\alpha_{a}+\alpha_{c})}\log\left(\frac{k_{c}c_{e}(c_{s,max}-\iota^{-1}(\Delta\phi))}{k_{a}\iota^{-1}(\Delta\phi)}\right)+u_{eq}.$
(15)

Plugging this expression into (2), we obtain

$j=F\iota(c_{s})\left[e^{\frac{\alpha_{a}F}{RT}\left(\Delta\phi-
u_{eq}\right)}\left(\frac{k_{c}c_{e}(c_{s,max}-c_{s})}{k_{a}c_{s}}\right)^{-\frac{\alpha_{a}}{\alpha_{a}+\alpha_{c}}}-e^{\frac{\alpha_{c}F}{RT}\left(\Delta\phi-
u_{eq}\right)}\left(\frac{k_{c}c_{e}(c_{s,max}-c_{s})}{k_{a}c_{s}}\right)^{\frac{\alpha_{c}}{\alpha_{a}+\alpha_{c}}}\right],$
(16)

which should vanish when $c_{s}\rightarrow 0$ as the electrode cannot be
depleted at vanishing concentrations, and when $c_{s}\rightarrow c_{s,max}$ as
the electrode cannot continue to be filled once the maximum concentration is
reached. In order to overcome the divergence of the factors in the parentheses
in these limits, we must require the function $\iota(c_{s})$ to vanish
sufficiently rapidly as $c_{s}\rightarrow 0,c_{s,max}$. More specifically,
expression (16) will vanish in these limits provided

$\iota(c_{s})=\begin{cases}\mathcal{O}\left[\left(\frac{k_{a}c_{s}}{k_{c}c_{e}(c_{s,max}-c_{s})}\right)^{\omega}\right],&\omega>\frac{\alpha_{c}}{\alpha_{a}+\alpha_{c}},\quad\textrm{as}\
c_{s}\to 0,\\\
\mathcal{O}\left[\left(\frac{k_{c}c_{e}(c_{s,max}-c_{s})}{k_{a}c_{s}}\right)^{\gamma}\right],&\gamma>\frac{\alpha_{a}}{\alpha_{a}+\alpha_{c}},\quad\textrm{as}\
c_{s}\to c_{s,max}\\\ \end{cases}.$ (17)

Condition (17) must be incorporated into the formulation of Problem P2.

##  Appendix B Computation of Gradients

We begin by computing the Gâteaux (directional) derivative of the error
functional (11) defined as
$\mathcal{J}_{2}^{\prime}(\iota(c_{s});\iota^{\prime}(c_{s}))=\lim\limits_{\epsilon\to
0}\epsilon^{-1}\left[\mathcal{J}_{2}(\iota(c_{s})+\epsilon\iota^{\prime}(c_{s}))-\mathcal{J}_{2}(\iota(c_{s}))\right]$,
where $\iota^{\prime}(c_{s})$ is an arbitrary perturbation of $\iota(c_{s})$,
and obtain

$\mathcal{J}_{2}^{\prime}(\iota(c_{s});\iota^{\prime}(c_{s}))=\int\limits_{0}^{t_{f}}\eta^{\prime}(\iota(c_{s}(t));\iota^{\prime}(c_{s}))\left[\eta(\iota(c_{s}(t)))-\overline{\eta}(\iota(c_{s}(t)))\right]\,dt.$
(18)

The Riesz representation theorem then asserts that the Gâteaux differential
(18) can be expressed as

$\mathcal{J}_{2}^{\prime}(\iota(c_{s});\iota^{\prime}(c_{s}))=\left\langle\nabla^{L^{2}}\mathcal{J},\iota^{\prime}\right\rangle_{L^{2}}=\int_{0}^{c_{s,max}}\nabla^{L^{2}}{\mathcal{J}}_{2}(c_{s})\cdot\iota^{\prime}(c_{s})\,dc,$
(19)

where $\nabla^{L^{2}}{\mathcal{J}}_{2}(c_{s})$ is the gradient of the error
functional computed with respect to the $L^{2}$ inner product
$\langle\cdot,\cdot\rangle_{L^{2}}$. To transform (18) into a form consistent
with (19), where $\iota^{\prime}(c_{s})$ appears as a factor and integration
is performed with respect to $c_{s}$, we define
$\iota_{\epsilon}(c_{s})=\iota(c_{s})+\epsilon\iota^{\prime}(c_{s})+\mathcal{O}(\epsilon^{2})$
so that the perturbed Butler-Volmer relation is
$j_{\epsilon}=\iota(c_{s})\psi(\Delta\phi,c_{s})+\epsilon\iota^{\prime}(c_{s})\psi(\Delta\phi,c_{s})+\mathcal{O}(\epsilon^{2})$.
Then, the corresponding cell overpotential, $\eta$, is

$\displaystyle\eta$
$\displaystyle=\psi^{-1}\left(\frac{I}{AF\left[\iota(c_{s})+\epsilon\iota^{\prime}(c_{s})+\mathcal{O}(\epsilon^{2})\right]}\right)$
$\displaystyle=\psi^{-1}\left(\frac{I}{AF\iota(c_{s})}\right)+\epsilon\underbrace{\left(\left[\psi^{-1}\left(\frac{I}{AF\iota(c_{s})}\right)\right]^{\prime}\frac{-I}{(AF)(\iota(c_{s}))^{2}}\right)\iota^{\prime}(c_{s})}_{\eta^{\prime}(\iota(c_{s}(t));\iota^{\prime}(c_{s}))}+\mathcal{O}(\epsilon^{2}).$
(20)

Using (20) in (18), we obtain

$\displaystyle\mathcal{J}_{2}^{\prime}(\iota(c_{s});\iota^{\prime}(c_{s}))$
$\displaystyle=\int\limits_{0}^{t_{f}}\frac{\frac{-I\iota^{\prime}(c_{s}(t))}{(AF)(\iota(c_{s}(t)))^{2}}}{\psi^{\prime}\left(\frac{I}{AF\iota(c_{s}(t))}\right)}\left[\eta(\iota(c_{s}(t)))-\overline{\eta}(\iota(c_{s}(t)))\right]\,dt$
$\displaystyle=\int\limits_{0}^{c_{s,max}}\frac{1}{\psi^{\prime}\left(\frac{I}{AF\iota(c_{s})}\right)}\frac{-I\iota^{\prime}(c_{s})}{(AF)(\iota(c_{s}))^{2}}\left(\eta(\iota(c_{s}))-\overline{\eta}(\iota(c_{s}))\right)\,\left|\frac{dt}{dc_{s}}\right|\,dc_{s},$
(21)

where the integration variable has been changed from $t$ to $c_{s}$ using the
relation $c_{s}=c_{s}(t)$ defined by the numerical solution of the SPM. We
note that (21) is now in a form consistent with the Riesz representation (19),
which allows us to define the $L^{2}$ gradient as

$\nabla^{L^{2}}\mathcal{J}_{2}=\left(\frac{1}{\psi^{\prime}\left(\frac{I}{AF\iota(c_{s})}\right)}\frac{-I}{(AF)(\iota(c_{s}))^{2}}\right)\left(\eta(\iota(c_{s}))-\overline{\eta}(\iota(c_{s}))\right)\,\left|\frac{dt}{dc_{s}}\right|,\qquad
c_{s}\in[0,c_{s,max}].$ (22)

We observe that while expression (22) is defined for $c_{s}\in[0,c_{s,max}]$,
it will be nonzero only for $c_{s}\in[c_{<},c_{>}]={\mathcal{I}}$, where
$c_{<}=\inf_{t\in[0,t_{f}]}c_{s}(t)$, $c_{>}=\sup_{t\in[0,t_{f}]}c_{s}(t)$ and
${\mathcal{I}}$ is referred to as the identifiability interval. In other
words, the sensitivity information can be obtained only for concentration
values in the identifiability interval, which is the range of concentrations
attained in a given experiment. In addition, expression (22) can produce
discontinuous functions of $c_{s}$ which, when used in (12), would not produce
constitutive relations with the required regularity, cf. Assumption A1 in
Section 2. In the light of expression (12), the regularity of the optimal
constitutive relation $\widehat{\iota}(c_{s})$ is determined by the smoothness
of the gradients $\nabla{\mathcal{J}}_{2}(c_{s})$.

Gradients with the required regularity are constructed using the Riesz
identity (19) extended to the Sobolev space $H^{1}([0,c_{s,max}])$ of
functions with square-integrable distributional derivatives, namely,

$\mathcal{J}_{2}^{\prime}(\iota(c_{s});\iota^{\prime}(c_{s}))=\left\langle\nabla^{L^{2}}\mathcal{J}_{2},\iota^{\prime}\right\rangle_{L^{2}}=\Big{\langle}\nabla\mathcal{J}_{2},\iota^{\prime}\Big{\rangle}_{H^{1}},$
(23)

which makes it possible to identify the Gâteaux differential (21) with the
$H^{1}$ inner product defined as

$\Big{\langle}\nabla\mathcal{J}_{2},\iota^{\prime}\Big{\rangle}_{H^{1}}=\int\limits_{0}^{c_{s,max}}\left[\nabla\mathcal{J}_{2}\,\iota^{\prime}+\ell^{2}\frac{d\left(\nabla\mathcal{J}_{2}\right)}{dc_{s}}\frac{d\iota^{\prime}}{dc_{s}}\right]\,dc_{s},$
(24)

where $\ell\in\mathbb{R}$ is a parameter with the meaning of a characteristic
concentration. We note that as long as $0<\ell<\infty$, inner products (24)
with different values of $\ell$ are equivalent in the sense of norm
equivalence. However, the choice of the value of parameter $\ell$ is important
in computations as it controls the degree of regularization in the solution of
Problem P2. Performing integration by parts with respect to $c_{s}$ on the
second term in (24), we obtain

$\Big{\langle}\nabla\mathcal{J}_{2},\iota^{\prime}\Big{\rangle}_{H^{1}}=\int\limits_{0}^{c_{s,max}}\left(Id-\ell^{2}\frac{d^{2}}{dc_{s}^{2}}\right)\nabla\mathcal{J}_{2}\,\iota^{\prime}\,dc_{s}-\ell^{2}\frac{d\left(\nabla\mathcal{J}_{2}\right)}{dc_{s}}\iota^{\prime}\Bigg{|}_{0}^{c_{s,max}}.$
(25)

Noting that the perturbation of the exchange current $\iota^{\prime}(c_{s})$
(i.e., the “test” function) is arbitrary and the Sobolev gradient
$\nabla{\mathcal{J}}_{2}(c_{s})$ needs to vanish at the endpoints, such that
$\iota^{\prime}(0)=\iota^{\prime}(c_{s,max})=0$ (cf. Assumptions A2–A3 in
Section 2 and Appendix A), relation (25) implies

$\displaystyle\left(Id-\ell^{2}\frac{d^{2}}{dc_{s}^{2}}\right)\nabla\mathcal{J}_{2}$
$\displaystyle=\nabla^{L^{2}}\mathcal{{\mathcal{J}}}_{2}$
$\displaystyle\textrm{for}\ c_{s}\in(0,c_{s,max}),$ (26a)
$\displaystyle\nabla{\mathcal{J}}_{2}$ $\displaystyle=0$
$\displaystyle\textrm{at}\ c_{s}=0,c_{s,max},$ (26b)

which shows that the Sobolev gradient can be obtained by solving a boundary-
value problem with the $L^{2}$ gradient (22) appearing as a source term. We
add that this formalism allows us to propagate the sensitivity information
outside the identifiability interval ${\mathcal{I}}$ such that the Sobolev
gradient $\nabla{\mathcal{J}}_{2}(c_{s})$ is a continuous functions not
identically equal to zero for
$c_{s}\in[0,c_{s,max}]\>\backslash\>{\mathcal{I}}$. It is known72 that
extraction of gradients in spaces of smoother functions such as $H^{1}$ can be
interpreted as low-pass filtering of the $L^{2}$ gradients with parameter
$\ell$ acting as the cut-off length-scale of the filter. The value of $\ell$
can significantly affect the rate of convergence of the iterative procedure
(12).

##  Appendix C Experimental Parameters

###  C.1 Slow charge rates

Parameter | Units | Anode | Cathode
---|---|---|---
electrode thickness, $L$, |  $\mu$m | $50$ | $50$
electrode particle radius, $R_{p}$ |  $\mu$m | $9$ | $6.5$
electrode cross-section area, $A$ | m2 | $0.000551$ | $0.000551$
volume fraction of electrolyte, $\epsilon_{l}$ | unitless | $0.1$ | $0.1$
Brunauer-Emmett-Teller surface area, $b_{0}$ | m-1 | $3.0\times 10^{5}$ | $4.1538\times 10^{5}$
conductivity in solid, $\sigma_{0}$ | Sm-1 | $14$ | $68.1$
permeability factor of electrode, $\mathcal{B}_{0}$ | unitless | $0.0493$ | $0.0515$
reaction rate constant, $k_{0}$ | m5/2s-1mol-1/2 | $1.8742\times 10^{-10}$ | $4.9375\times 10^{-11}$
maximum concentration of Lithium in solid, $c_{s,max}$ | mol m-3 | $35920$ | $42580$
initial concentration of Lithium in electrode, $c_{0}$ | mol m-3 | $1000$ | $1000$
activation energy, $E_{i}$ | J mol-1 | $30300$ | $80600$
Table 5: Electrode specific experimental parameters for slow charge rate experiments29. Parameter | Units | Value
---|---|---
temperature, $T$ | K | $295.15$
relative temperature, $\hat{T}$ | K | $296$
universal gas constant, $R$ | J mol-1 K-1 | $8.3144$
C$/10$ current, $I_{0}$ | A | $2.29\times 10^{-3}$
C$/3$ current, $I_{0}$ | A | $0.0076$
transference number, $t^{+}$ | unitless | $0.26$
Table 6: Global experimental parameters for slow charge rate experiments29.

Diffusivities in the electrodes are assumed to depend on Arrhenius temperature
and are given by

$\displaystyle D_{s}(c_{s,i},T)$
$\displaystyle=D_{s}(c_{s,i})e^{\left(\frac{E_{i}}{R\hat{T}}-\frac{E_{i}}{RT}\right)},$
$\displaystyle D_{s}(c_{s,a})$ $\displaystyle=8.4\times
10^{-9}e^{-11.3\left(\frac{c_{s,a}}{c_{s,a,max}}\right)}+8.2\times 10^{-12},$
$\displaystyle D_{s}(c_{s,c})$ $\displaystyle=3.7\times 10^{-13}-3.4\times
10^{-13}e^{-12\left(\frac{c_{s,c}}{c_{s,c,max}}-0.62\right)}$

with units of m2/s.

###  C.2 Moderate discharge rates

Parameter | Units | Anode | Cathode
---|---|---|---
electrode thickness, $L$, |  $\mu$m | $85.2$ | $75.6$
electrode particle radius, $R_{p}$ |  $\mu$m | $5.86$ | $5.22$
electrode cross-section area, $A$ | m2 | $0.10465$ | $0.1027$
volume fraction of electrolyte, $\epsilon_{l}$ | unitless | $0.25$ | $0.335$
Brunauer-Emmett-Teller surface area, $b_{0}$ | m-1 | $3.8396\times 10^{5}$ | $3.8218\times 10^{5}$
conductivity in solid, $\sigma_{0}$ | Sm-1 | $215$ | $0.18$
permeability factor of electrode, $\mathcal{B}_{0}$ | unitless | $0.0177$ | $0.0701$
reaction rate Constant, $k_{0}$ | m5/2s-1mol-1/2 | $0.0672\times 10^{-10}$ | $0.3545\times 10^{-10}$
maximum concentration of Lithium in solid, $c_{s,max}$ | mol m-3 | $33133$ | $63104$
initial concentration of Lithium in electrode, $c_{0}$ | mol m-3 | $1000$ | $1000$
Table 7: Electrode specific experimental parameters for moderate discharge rate experiments30. Parameter | Units | Value
---|---|---
temperature, $T$ | K | $298.15$
C$/2$ current, $I_{0}$ | A | $2.5$
$1$C current, $I_{0}$ | A | $5$
Transference Number, $t^{+}$ | unitless | $0.259$
Table 8: Global experimental parameters for moderate discharge rate
experiments30.

Transport in the electrodes is assumed to be governed by Fickian diffusion
with

$\displaystyle D_{s}(c_{s,a})$ $\displaystyle=3.3\times 10^{-10},$
$\displaystyle D_{s}(c_{s,c})$ $\displaystyle=4\times 10^{-11}$

and with units of m2/s.

##  Appendix D Checklist To Report Theoretical Battery Studies

Here, we complete the checklist for the minimal information set 73.

  1. 1.

Have you provided all assumptions, theory, governing equations, initial and
boundary conditions, material properties (e.g., open-circuit potential) with
appropriate precision and literature sources, constant states (e.g.,
temperature), etc.?

Remarks: Yes, all elements of the model have been stated in the main text and
the supplemental information.

  2. 2.

If the calculations have a probabilistic component (e.g., Monte Carlo, initial
configuration in Molecular Dynamics, etc.), did you provide statistics (mean,
standard deviation, confidence interval, etc.) from multiple ($\geq 3$) runs
of a representative case?

Remarks: Yes, Monte Carlo simulations were run with $3000$ samples. The
methodology is reported in the main text and confidence intervals are included
on the relevant figures.

  3. 3.

If data-driven calculations are performed (e.g., machine learning), did you
specify dataset origin, the rationale behind choosing it, what information it
contains, and the specific portion of it being utilized? Have you described
the thought process for choosing a specific modeling paradigm?

Remarks: While the methodology used is data driven, it relies on calculus of
variations rather than machine learning. The experimental data is described in
Section $3$ whereas the optimization scheme is outlined in Algorithms 1 and 2.

  4. 4.

Have you discussed all sources of potential uncertainty, variability, and
errors in the modeling results and their impact on quantitative results and
qualitative trends? Have you discussed the sensitivity of modeling (and
numerical) inputs such as material properties, time step, domain size, neural
network architecture, etc. where they are variable or uncertain?

Remarks: Yes, this discussion is provided in Sections 4.4 and 6 in the main
text.

  5. 5.

Have you sufficiently discussed new or not widely familiar terminology and
descriptors for clarity? Did you use these terms in their appropriate context
to avoid misinterpretation? Enumerate these terms in the “Remarks”.

Remarks: Yes, the terminology is consistent with the conventions used in the
research fields.

## References

  * Wang et al. 2022 Wang, A.; O’Kane, S. E.; Planella, F. B.; Le Houx, J.; O’Regan, K.; Zyskin, M.; Edge, J. S.; Monroe, C.; Cooper, S.; Howey, D. A. et al. Review of parameterisation and a novel database (LiionDB) for continuum Li-ion battery models. _Progress in Energy_ 2022, _4_ , 032004
  * Lv et al. 2022 Lv, C.; Zhou, X.; Zhong, L.; Yan, C.; Srinivasan, M.; Seh, Z. W.; Liu, C.; Pan, H.; Li, S.; Wen, Y. et al. Machine learning: an advanced platform for materials development and state prediction in lithium-ion batteries. _Advanced Materials_ 2022, _34_ , 2101474
  * Houchins and Viswanathan 2020 Houchins, G.; Viswanathan, V. An accurate machine-learning calculator for optimization of Li-ion battery cathodes. _The Journal of Chemical Physics_ 2020, _153_ , 054124
  * Chemali et al. 2018 Chemali, E.; Kollmeyer, P. J.; Preindl, M.; Emadi, A. State-of-charge estimation of Li-ion batteries using deep neural networks: A machine learning approach. _Journal of Power Sources_ 2018, _400_ , 242–255
  * Klett et al. 2012 Klett, M.; Giesecke, M.; Nyman, A.; Hallberg, F.; Lindström, R. W.; Lindbergh, G.; Furo, I. Quantifying Mass Transport during Polarization in a Li Ion Battery Electrolyte by in Situ Li NMR Imaging. _J. Am. Chem. Soc._ 2012, _134_ , 14654–14657
  * Sethurajan et al. 2015 Sethurajan, A. K.; Krachkovskiy, S. A.; Halalay, I. C.; Goward, G. R.; Protas, B. Accurate characterization of ion transport properties in binary symmetric electrolytes using in situ NMR imaging and inverse modeling. _The Journal of Physical Chemistry B_ 2015, _119_ , 12238–12248
  * Richardson et al. 2018 Richardson, G.; Foster, J. M.; Sethurajan, A. K.; Krachkovskiy, S. A.; Halalay, I. C.; Goward, G. R.; Protas, B. The effect of ionic aggregates on the transport of charged species in lithium electrolyte solutions. _Journal of The Electrochemical Society_ 2018, _165_ , H561–H567
  * Sethurajan et al. 2019 Sethurajan, A. K.; Foster, J. M.; Richardson, G.; Krachkovskiy, S. A.; Bazak, J. D.; Goward, G. R.; Protas, B. Incorporating Dendrite Growth into Continuum Models of Electrolytes: Insights from NMR Measurements and Inverse Modeling. _Journal of The Electrochemical Society_ 2019, _166_ , A1591–A1602
  * Escalante et al. 2020 Escalante, J. M.; Ko, W.; Foster, J. M.; Krachkovskiy, S.; Goward, G.; Protas, B. Discerning models of phase transformations in porous graphite electrodes: Insights from inverse modelling based on MRI measurements. _Electrochimica Acta_ 2020, _349_ , 136290
  * Escalante et al. 2021 Escalante, J. M.; Sahu, S.; Foster, J. M.; Protas, B. On Uncertainty Quantification in the Parametrization of Newman-type Models of Lithium-ion Batteries. _Journal of The Electrochemical Society_ 2021, _168_ , 110519
  * Chen et al. 2019 Chen, W.; Liang, J.; Yang, Z.; Li, G. A review of lithium-ion battery for electric vehicle applications and beyond. _Energy Procedia_ 2019, _158_ , 4363–4368
  * Blomgren 2016 Blomgren, G. E. The development and future of lithium ion batteries. _Journal of The Electrochemical Society_ 2016, _164_ , A5019
  * Newman and Thomas-Alyea 2004 Newman, J.; Thomas-Alyea, K. E. _Electrochemical Systems_ ; John Wiley and Sons, 2004
  * Hamelers et al. 2011 Hamelers, H. V.; Ter Heijne, A.; Stein, N.; Rozendal, R. A.; Buisman, C. J. Butler–Volmer–Monod model for describing bio-anode polarization curves. _Bioresource technology_ 2011, _102_ , 381–387
  * Compton and Banks 2011 Compton, R. G.; Banks, C. E. _Understanding voltammetry, 2nd ed._ ; Imperial College Press London, 2011
  * Oyarzun et al. 2021 Oyarzun, D. I.; Zhan, C.; Hawks, S. A.; Cerón, M. R.; Kuo, H. A.; Loeb, C. K.; Aydin, F.; Pham, T. A.; Stadermann, M.; Campbell, P. G. Unraveling the Ion Adsorption Kinetics in Microporous Carbon Electrodes: A Multiscale Quantum-Continuum Simulation and Experimental Approach. _ACS Applied Materials & Interfaces_ 2021, _13_ , 23567–23574
  * Hodgkin and Huxley 1952 Hodgkin, A. L.; Huxley, A. F. A quantitative description of membrane current and its application to conduction and excitation in nerve. _The Journal of physiology_ 1952, _117_ , 500
  * Banks and Kunisch 2012 Banks, H. T.; Kunisch, K. _Estimation techniques for distributed parameter systems_ ; Springer Science & Business Media, 2012
  * Tarantola 2005 Tarantola, A. _Inverse Problem Theory and Methods for Model Parameter Estimation_ ; Society for Industrial and Applied Mathematics, 2005
  * Isakov 1998 Isakov, V. _Inverse Problems for Partial Differential Equations_ ; Applied Mathematical Sciences 127; Springer, 1998
  * Bukshtynov et al. 2011 Bukshtynov, V.; Volkov, O.; Protas, B. On Optimal Reconstruction of Constitutive Relations. _Physica D: Nonlinear Phenomena_ 2011, _240_ , 1228–1244
  * Sethurajan et al. 2019 Sethurajan, A.; Krachkovskiy, S.; Goward, G.; Protas, B. Bayesian uncertainty quantification in inverse modeling of electrochemical systems. _Journal of Computational Chemistry_ 2019, _40_ , 740–752
  * Bukshtynov and Protas 2013 Bukshtynov, V.; Protas, B. Optimal reconstruction of material properties in complex multiphysics phenomena. _Journal of Computational Physics_ 2013, _242_ , 889 – 914
  * Dickinson and Wain 2020 Dickinson, E. J.; Wain, A. J. The Butler-Volmer equation in electrochemical theory: Origins, value, and practical application. _Journal of Electroanalytical Chemistry_ 2020, _872_ , 114145
  * Latz and Zausch 2013 Latz, A.; Zausch, J. Thermodynamic derivation of a Butler–Volmer model for intercalation in Li-ion batteries. _Electrochimica Acta_ 2013, _110_ , 358–362
  * Bard and Faulkner 2001 Bard, A. J.; Faulkner, L. R. _Electrochemical methods: Fundamentals and applications_ ; Wiley New York, 2001
  * Fletcher 2009 Fletcher, S. Tafel slopes from first principles. _Journal of Solid State Electrochemistry_ 2009, _13_ , 537–549
  * Richardson et al. 2020 Richardson, G. W.; Foster, J. M.; Ranom, R.; Please, C. P.; Ramos, A. M. Charge transport modelling of Lithium-ion batteries. _European Journal of Applied Mathematics_ 2020, 1–49
  * Sanders et al. 2022 Sanders, K. J.; Aguilera, A. R.; Keffer, J. R.; Balcom, B. J.; Halalay, I. C.; Goward, G. R. Transient lithium metal plating on graphite: Operando 7Li nuclear magnetic resonance investigation of a battery cell using a novel RF probe. _Carbon_ 2022, 377–385
  * Chen et al. 2020 Chen, C.-H.; Planella, F. B.; O’regan, K.; Gastol, D.; Widanage, W. D.; Kendrick, E. Development of experimental techniques for parameterization of multi-scale lithium-ion battery models. _Journal of The Electrochemical Society_ 2020, _167_ , 080534
  * Richardson et al. 2020 Richardson, G.; Korotkin, I.; Ranom, R.; Castle, M.; Foster, J. Generalised single particle models for high-rate operation of graded lithium-ion electrodes: systematic derivation and validation. _Electrochimica Acta_ 2020, _339_ , 135862
  * Korotkin et al. 2021 Korotkin, I.; Sahu, S.; O’Kane, S. E.; Richardson, G.; Foster, J. M. DandeLiion v1: An extremely fast solver for the Newman model of lithium-ion battery (dis) charge. _Journal of The Electrochemical Society_ 2021, _168_ , 060544
  * Wu et al. 2012 Wu, S.-L.; Zhang, W.; Song, X.; Shukla, A. K.; Liu, G.; Battaglia, V.; Srinivasan, V. High rate capability of Li (Ni1/3Mn1/3Co1/3) O2 electrode for Li-ion batteries. _Journal of The Electrochemical Society_ 2012, _159_ , A438
  * Nguyen et al. 2022 Nguyen, T.-T.; Delobel, B.; Berthe, M.; Fleutot, B.; Demortière, A.; Delacourt, C. Mathematical Modeling of Energy-Dense NMC Electrodes: I. Determination of Input Parameters. _Journal of The Electrochemical Society_ 2022, _169_ , 040546
  * Doyle et al. 1993 Doyle, M.; Fuller, T. F.; Newman, J. Modeling of Galvanostatic Charge and Discharge of the Lithium/Polymer/Insertion Cell. _Journal of The Electrochemical Society_ 1993, _140_ , 1526–1533
  * Fuller et al. 1994 Fuller, T. F.; Doyle, M.; Newman, J. Simulation and optimization of the dual lithium ion insertion cell. _Journal of the Electrochemical Society_ 1994, _141_ , 1
  * Mistry et al. 2018 Mistry, A. N.; Smith, K.; Mukherjee, P. P. Secondary-phase stochastics in lithium-ion battery electrodes. _ACS applied materials & interfaces_ 2018, _10_ , 6317–6326
  * Mistry et al. 2021 Mistry, A.; Trask, S.; Dunlop, A.; Jeka, G.; Polzin, B.; Mukherjee, P. P.; Srinivasan, V. Quantifying negative effects of carbon-binder networks from electrochemical performance of porous li-ion electrodes. _Journal of The Electrochemical Society_ 2021, _168_ , 070536
  * Mistry et al. 2020 Mistry, A.; Usseglio-Viretta, F. L.; Colclasure, A.; Smith, K.; Mukherjee, P. P. Fingerprinting redox heterogeneity in electrodes during extreme fast charging. _Journal of The Electrochemical Society_ 2020, _167_ , 090542
  * Hein et al. 2020 Hein, S.; Danner, T.; Westhoff, D.; Prifling, B.; Scurtu, R.; Kremer, L.; Hoffmann, A.; Hilger, A.; Osenberg, M.; Manke, I. et al. Influence of conductive additives and binder on the impedance of lithium-ion battery electrodes: effect of morphology. _Journal of The Electrochemical Society_ 2020, _167_ , 013546
  * Ferraro et al. 2020 Ferraro, M. E.; Trembacki, B. L.; Brunini, V. E.; Noble, D. R.; Roberts, S. A. Electrode mesoscale as a collection of particles: coupled electrochemical and mechanical analysis of NMC cathodes. _Journal of The Electrochemical Society_ 2020, _167_ , 013543
  * Mistry et al. 2022 Mistry, A.; Heenan, T.; Smith, K.; Shearing, P.; Mukherjee, P. P. Asphericity can cause nonuniform lithium intercalation in battery active particles. _ACS Energy Letters_ 2022, _7_ , 1871–1879
  * Verbrugge and Koch 1996 Verbrugge, M. W.; Koch, B. J. Modeling lithium intercalation of single-fiber carbon microelectrodes. _Journal of the Electrochemical Society_ 1996, _143_ , 600
  * Johnson et al. 2022 Johnson, I. D.; Mistry, A. N.; Yin, L.; Murphy, M.; Wolfman, M.; Fister, T. T.; Lapidus, S. H.; Cabana, J.; Srinivasan, V.; Ingram, B. J. Unconventional Charge Transport in MgCr2O4 and Implications for Battery Intercalation Hosts. _Journal of the American Chemical Society_ 2022, _144_ , 14121–14131
  * Amin et al. 2015 Amin, R.; Ravnsbæk, D. B.; Chiang, Y.-M. Characterization of electronic and ionic transport in Li1-xNi0. 8Co0. 15Al0. 05O2 (NCA). _Journal of the electrochemical society_ 2015, _162_ , A1163
  * Amin and Chiang 2016 Amin, R.; Chiang, Y.-M. Characterization of electronic and ionic transport in Li1-xNi0. 33Mn0. 33Co0. 33O2 (NMC333) and Li1-xNi0. 50Mn0. 20Co0. 30O2 (NMC523) as a function of Li content. _Journal of The Electrochemical Society_ 2016, _163_ , A1512
  * Amin et al. 2008 Amin, R.; Maier, J.; Balaya, P.; Chen, D.; Lin, C. Ionic and electronic transport in single crystalline LiFePO4 grown by optical floating zone technique. _Solid State Ionics_ 2008, _179_ , 1683–1687
  * Santhanagopalan et al. 2006 Santhanagopalan, S.; Guo, Q.; Ramadass, P.; White, R. E. Review of models for predicting the cycling performance of lithium ion batteries. _Journal of power sources_ 2006, _156_ , 620–628
  * Schmalstieg et al. 2018 Schmalstieg, J.; Rahe, C.; Ecker, M.; Sauer, D. U. Full cell parameterization of a high-power lithium-ion battery for a physico-chemical model: Part I. Physical and electrochemical parameters. _Journal of The Electrochemical Society_ 2018, _165_ , A3799
  * Marquis et al. 2019 Marquis, S. G.; Sulzer, V.; Timms, R.; Please, C. P.; Chapman, S. J. An asymptotic derivation of a single particle model with electrolyte. _Journal of The Electrochemical Society_ 2019, _166_ , A3693
  * Adams and Fournier 2005 Adams, R. A.; Fournier, J. F. _Sobolev Spaces_ ; Elsevier, 2005
  * Press et al. 1986 Press, W. H.; Flanner, B. P.; Teukolsky, S. A.; Vetterling, W. T. _Numerical Recipes: the Art of Scientific Computations_ ; Cambridge University Press, 1986
  * Devore 2011 Devore, J. L. _Probability and Statistics for Engineering and the Sciences_ ; Cengage learning, 2011
  * Ecker et al. 2015 Ecker, M.; Tran, T. K. D.; Dechent, P.; Käbitz, S.; Warnecke, A.; Sauer, D. U. Parameterization of a Physico-Chemical Model of a Lithium-Ion Battery: I. Determination of Parameters. _Journal of The Electrochemical Society_ 2015, _162_ , A1836–A1848
  * Park et al. 2010 Park, M.; Zhang, X.; Chung, M.; Less, G. B.; Sastry, A. M. A review of conduction phenomena in Li-ion batteries. _Journal of Power Sources_ 2010, _195_ , 7904–7929
  * Jow et al. 2018 Jow, T. R.; Delp, S. A.; Allen, J. L.; Jones, J.-P.; Smart, M. C. Factors limiting Li+ charge transfer kinetics in Li-ion batteries. _Journal of the electrochemical society_ 2018, _165_ , A361
  * Christensen and Newman 2006 Christensen, J.; Newman, J. Stress generation and fracture in lithium insertion materials. _Journal of Solid State Electrochemistry_ 2006, _10_ , 293–319
  * O’Kane et al. 2022 O’Kane, S. E.; Ai, W.; Madabattula, G.; Alonso-Alvarez, D.; Timms, R.; Sulzer, V.; Edge, J. S.; Wu, B.; Offer, G. J.; Marinescu, M. Lithium-ion battery degradation: how to model it. _Physical Chemistry Chemical Physics_ 2022, _24_ , 7909–7922
  * Zinth et al. 2014 Zinth, V.; Von Lüders, C.; Hofmann, M.; Hattendorff, J.; Buchberger, I.; Erhard, S.; Rebelo-Kornmeier, J.; Jossen, A.; Gilles, R. Lithium plating in lithium-ion batteries at sub-ambient temperatures investigated by in situ neutron diffraction. _Journal of Power Sources_ 2014, _271_ , 152–159
  * Uhlmann et al. 2015 Uhlmann, C.; Illig, J.; Ender, M.; Schuster, R.; Ivers-Tiffée, E. In situ detection of lithium metal plating on graphite in experimental cells. _Journal of Power Sources_ 2015, _279_ , 428–438
  * Wandt et al. 2018 Wandt, J.; Jakes, P.; Granwehr, J.; Eichel, R.-A.; Gasteiger, H. A. Quantitative and time-resolved detection of lithium plating on graphite anodes in lithium ion batteries. _Materials Today_ 2018, _21_ , 231–240
  * Yang et al. 2017 Yang, X.-G.; Leng, Y.; Zhang, G.; Ge, S.; Wang, C.-Y. Modeling of lithium plating induced aging of lithium-ion batteries: Transition from linear to nonlinear aging. _Journal of Power Sources_ 2017, _360_ , 28–40
  * Yang et al. 2018 Yang, X.-G.; Ge, S.; Liu, T.; Leng, Y.; Wang, C.-Y. A look into the voltage plateau signal for detection and quantification of lithium plating in lithium-ion cells. _Journal of Power Sources_ 2018, _395_ , 251–261
  * O’Kane et al. 2020 O’Kane, S. E.; Campbell, I. D.; Marzook, M. W.; Offer, G. J.; Marinescu, M. Physical origin of the differential voltage minimum associated with lithium plating in Li-ion batteries. _Journal of The Electrochemical Society_ 2020, _167_ , 090540
  * Sahu and Foster 2023 Sahu, S.; Foster, J. M. A continuum model for lithium plating and dendrite formation in lithium-ion batteries: Formulation and validation against experiment. _Journal of Energy Storage_ 2023, _60_ , 106516
  * Levi et al. 2003 Levi, M.; Wang, C.; Markevich, E.; Aurbach, D.; Chvoj, Z. Noteworthy electroanalytical features of the stage 4 to stage 3 phase transition in lithiated graphite. _Journal of Solid State Electrochemistry_ 2003, _8_ , 40–43
  * Verbrugge et al. 2022 Verbrugge, M. W.; Baker, D. R.; Chen, S.; He, M.; Cai, M. Cation Mixing and Capacity Loss in Li—— Ni0. 6Mn0. 2Co0. 2O2 Cells: Experimental Investigation and Application of the Multi-Site, Multi-Reaction Model. _Frontiers in Energy Research_ 2022, _10_ , 671
  * Baker et al. 2019 Baker, D. R.; Verbrugge, M. W.; Gu, W. Multi-species, multi-reaction model for porous intercalation electrodes: Part II. Model-experiment comparisons for linear-sweep voltammetry of spinel lithium manganese oxide electrodes. _Journal of The Electrochemical Society_ 2019, _166_ , A521
  * Shockley and Read Jr 1952 Shockley, W.; Read Jr, W. Statistics of the recombinations of holes and electrons. _Physical review_ 1952, _87_ , 835
  * Hall 1952 Hall, R. N. Electron-hole recombination in germanium. _Physical review_ 1952, _87_ , 387
  * Bousse et al. 1983 Bousse, L.; De Rooij, N. F.; Bergveld, P. The influence of counter-ion adsorption on the $\psi$0/pH characteristics of insulator surfaces. _Surface science_ 1983, _135_ , 479–496
  * Protas et al. 2004 Protas, B.; Bewley, T. R.; Hagen, G. A Computational Framework for the Regularization of Adjoint Analysis in Multiscale PDE Systems. _J. Comput. Phys._ 2004, _195_ , 49–89
  * Mistry et al. 2021 Mistry, A.; Verma, A.; Sripad, S.; Ciez, R.; Sulzer, V.; Brosa Planella, F.; Timms, R.; Zhang, Y.; Kurchin, R.; Dechent, P. et al. A Minimal Information Set To Enable Verifiable Theoretical Battery Research. _ACS Energy Letters_ 2021, _6_ , 3831–3835