\section{Introduction}

With the electrification of vehicles and the spread of portable electronic devices, lithium-ion batteries have become a key technology for energy storage. This has motivated a quest to design batteries that can store more energy, last longer, and operate efficiently and safely. Mathematical models are an invaluable tool for battery design and control as they provide a cheap, safe and fast alternative to experiments. 

There exist multiple approaches to deterministic models for lithium-ion batteries, and they can usually be classified into two broad categories: equivalent-circuit models and physics-based models. Equivalent-circuit models describe the battery by assuming it has the structure of a particular electrical circuit. The morphology of this circuit is a modelling assumption, and the values of the different parameters need to be determined from fitting the model to experimental data (see \cite{Hu2012,Widanage2016} for examples). Equivalent-circuit models are widely used in battery management systems (BMS) because they are computationally cheap to simulate \cite{Farmann2016}. However, they offer little insight on the battery internal states and therefore it is very complicated to extend these models to account for additional physics. A detailed description of equivalent-circuit models can be found in \cite{Plett2015}.

In this work, we focus on physics-based models. These models are derived from fundamental physical laws and have the advantage that they provide significant insight on the internal states of the battery. However, due to their complexity, they are computationally more expensive than equivalent-circuit models. The most widely used physics-based model is the Doyle-Fuller-Newman (DFN) model, originally posed in \cite{Fuller1994,Fuller1994a}. This model accounts for conservation of mass and charge in both electrode and electrolyte, and the four equations are coupled together through the reaction current density term that describes the electrode-electrolyte intercalation reaction. One of the key aspects for the success of the DFN model is that it assumes a simplified geometry that captures the main features of a battery but is much more affordable to solve than trying to resolve all the phenomena at a microscale level \cite{Lu2020}. The geometry of the DFN model assumes that most of the variables (electrolyte concentration and potentials both in electrode and electrolyte) can be resolved at a macroscale level, in which the material is assumed to be homogeneous. In order to account for the effects of the porous material on the transport equations in the homogeneous material, effective transport parameters are used (often using the Bruggeman correlation \cite{Tjaden2016}). The transport of intercalated lithium is described by assuming that at each point of the electrode there is a representative spherical particle in which lithium diffuses. Even though it is not explicitly derived in \cite{Fuller1994}, this is a homogenised model and it can be formalised from the mathematical point of view (see \cite{Hunt2020} for details). A more detailed description of the DFN model can be found in the handbooks \cite{Plett2015,Newman2004}.

The DFN model is still quite complex as it involves a coupled system of differential-algebraic equations (DAEs). Hence, it is not suitable for applications where speed is crucial, such as battery optimization and control, so simpler models have been posed in the literature. One of the best-known examples is the Single Particle Model (SPM). The key idea of this model is that we can use one single particle in each electrode to represent the behaviour of all particles. SPMs may \cite{Moura2017} or may not \cite{Bizeray2018} include electrolyte dynamics. In both cases the model is much simpler than the DFN model as most of the variables of interest can be easily computed from the concentrations in the electrode particles and the electrolyte.

The SPMs are widely used, but in their origin there are still two main challenges. First, most of these models are posed on an \emph{ad hoc} basis and therefore the connection to the DFN (or similar) model is lost. This is far from ideal because, in the process of derivation, inconsistencies might be introduced and certain features lost. In addition, it makes it hard to add new physics in a consistent way. The second issue, closely related to the first, is that we can find multiple Single Particle Models in the literature following different formulations \cite{Moura2017,Bizeray2018,Kemper2013,Perez2016,Marquis2019,Richardson2020}. Over the last few years, there has been a surge in the asymptotic analysis of battery models \cite{Marquis2019,Richardson2020,Marquis2020,Richardson2012,Moyles2019asymptotic,Moyles2018experimental,Hennessy2019} with the goal of obtaining reduced models in a systematic manner. The different authors take different assumptions to reduce the models. One big difference is that \cite{Richardson2012,Moyles2019asymptotic,Moyles2018experimental,Hennessy2019} assume fast diffusion in the electrodes and thus do not reduce to SPM-type models, while \cite{Marquis2019,Richardson2020,Marquis2020} retain this feature. Given that our interest on this work is on SPM-type models, we will focus our attention on the latter.

The key idea behind asymptotic methods is to analyse the dimensionless groupings that appear in the model, determine which can be assumed to be small and perform an expansion in the limit where these parameters tend to zero to obtain a simpler model. More details on these methods can be found in the handbook \cite{Hinch1991}. These methods present multiple advantages. First, they are systematic so they can be applied to all the models, even the very complex ones. Second, they also ensure the consistency between the full and reduced models, and thus the reduced model satisfies the same physical laws as the full one (e.g. conservation of mass). And third, we can validate the assumptions taken and estimate the error of the reduced model \emph{a priori} (i.e. before even solving the model).

So far we referred to electrochemical models only, but a similar description holds for coupled thermal-electrochemical models. These are usually posed as a coupling between the DFN model, to describe the electrochemistry, and some sort of thermal model \cite{Plett2015,Gu2000,Srinivasan2003}. Due to their complexity, simpler models have been posed \emph{ad hoc}, usually based on SPM-type electrochemical models \cite{Perez2016}. The asymptotic analysis of thermal-electrochemical models is even more recent as it is a natural evolution of the asymptotic reduction of pure electrochemical models. To our knowledge, the only references in the literature on asymptotic reductions of thermal-electrochemical models are \cite{Hennessy2019,Timms2020}. In \cite{Hennessy2019} a multilayer cell is considered but with the electrode fast diffusion electrochemical model presented in \cite{Moyles2019asymptotic}; while in \cite{Timms2020} a model for single-layer pouch cells based on DFN and the SPMe of \cite{Marquis2019}.

The goal of this paper is to derive a SPM-type thermal-electrochemical model for a multi-layer cell, using an asymptotic method based on minimal assumptions. The resultant model, which we refer to as Thermal Single Particle Model with electrolyte (TSPMe), is validated against the full thermal DFN model and experimental data to show its validity and proof that it is suitable for practical applications. To the best of our knowledge, this model is the first instance of an asymptotically derived thermal-electrochemical model for a multi-layer battery that also accounts for finite diffusion in the electrode particles.


\section{Discussion of Doyler-Fuller-Newman model and Single Particle Models}\label{sec:SPM_discussion}

Before presenting the formulation for the Thermal Single Particle Model with electrolyte, we discuss the nature of Single Particle Models and compare it with that of the DFN model. Some of the most popular reduced models for lithium-ion batteries are Single Particle Models (SPM), and in the literature we can find several different models referred to as such \cite{Moura2017,Bizeray2018,Kemper2013,Perez2016,Marquis2019,Richardson2020}. These are usually based on slightly different sets of assumptions and therefore their formulations are different as well. Because the different models are referred to by the same name, it makes the term ``Single Particle Model'' very confusing. With the hope of shedding some light into this issue, in this section we discuss what are the fundamental features that constitute a Single Particle Model, understanding it as a type of model rather than a specific formulation. We refer to models of this kind as SPM-type models, which applies to both models with and without electrolyte dynamics. However, given that the model without electrolyte dynamics (SPM) is a particular case to the model with electrolyte dynamics (SPMe), we consider the latter in this analysis. By considering what are the essential features of an SPM-type model we can better understand and classify existing models in the literature.

\subsection{Electrochemical models}\label{sec:EC_models}
We start focusing on purely electrochemical models. By electrochemical models we refer to models that only account for mass and charge transport (in both the electrodes and the electrolyte) but do not include further physical effects such as thermal, mechanical or degradation. The classic example of an electrochemical model is the Doyle-Fuller-Newman (DFN) model \cite{Fuller1994}. To better understand the electrochemical models, we can represent them as a block diagram, as shown in Figure~\ref{fig:block_diagram_EC}. In these diagrams, we have an input (or inputs) that are known and want to calculate an output (or outputs) which are unknown. To convert an input into an output we use a model, represented in the diagram by a block. Note that this model could be anything that converts the input into an output, but here we consider it from the mathematical point of view, in which the model is composed of a set of differential and/or algebraic equations that need to be solved. The most typical setup for electrochemical models is to use the applied current density ($i_\mathrm{app}$) as input and calculate the terminal voltage ($V$) as an output. Other setups, such as inputting the voltage are common as well, but make the reduced models significantly more expensive as we discuss later. In addition, sometimes there are other variables that might be of interest and hence be part of the model output. 

\begin{figure}
    \centering
    \begin{subfigure}{0.49\textwidth}
      \centering
      \includegraphics[scale=1]{img/diagram_DFN.pdf}  
      \caption{DFN model}
    \end{subfigure}\hfill
    \begin{subfigure}{0.49\textwidth}
      \centering
      \includegraphics[scale=1]{img/diagram_SPMe.pdf}  
      \caption{SPM-type}
    \end{subfigure}
    \caption{Block diagram of standard electrochemical models: (a) the Doyler-Fuller-Newman model (DFN) and (b) the Single Particle Model type (SPM-type). Each block is a model (or part of a model) that converts the input variable into the output variable (or variables). In brackets, we define the the type of expressions that compose the model (e.g. PDEs, ODEs, explicit expressions) and, if applicable, the type of expressions after spatial discretisation (e.g. DAEs or ODEs). The colour of the blocks indicates the complexity of the model: red for a system of differential-algebraic equations (DAEs, high complexity), yellow for a system of ordinary differential equations (ODEs, medium complexity) and green for explicit expressions (low complexity).}
    \label{fig:block_diagram_EC}
\end{figure}

We start by considering the DFN model \cite{Fuller1994}. The model is composed of a system of parabolic and elliptic partial differential equations (PDEs), which are represented as a single block in the diagram indicating the need to solve the full model at every step output is required (see Figure~\ref{fig:block_diagram_EC}a). A very common method to solve partial differential equations is to use the method of lines \cite{Schiesser2016}, which consists of discretising the system in space, for example using finite volume methods \cite{LeVeque2002} or finite element methods \cite{Johnson2009}, and leaving the time derivative without discretising. This method transforms parabolic PDEs into ordinary differential equtions (ODEs) and elliptic PDEs into algebraic equations. Therefore, for the DFN model, which is a coupled system of parabolic and elliptic PDEs, we obtain after discretising a system of differential-algebraic equations (DAEs). Solving numerically a system of DAEs is a complex task, hence the model is very expensive to compute.

The key aspect of SPM-type models is that they break this block into simpler blocks (see Figure~\ref{fig:block_diagram_EC}b), and we can group these blocks into two main steps. The first step involves calculating from the input current the concentrations in the representative electrode particles ($c_k$ with $k \in \{\mathrm{n},\mathrm{p}\}$ for negative and positive electrode, respectively), and in the electrolyte ($c_\mathrm{e}$, if not neglected). All the other variables of interest can then be calculated from these concentrations using explicit expressions, or even closed-form expressions in some cases. The first step in the SPM-type models (``$c_k$ model'' and ``$c_\mathrm{e}$ model'' blocks in Figure~\ref{fig:block_diagram_EC}b) involves solving up to three decoupled partial differential equations: two PDEs for the concentration of lithium in the particles (one for each electrode) and one PDE for the ion concentration in the electrolyte. If the electrolyte dynamics are not considered we only need to solve the two PDEs for the particles. Because the PDEs are decoupled, we can solve them independently which results in a simple problem. To numerically solve these PDEs we can use the method of lines again, in order to obtain a system of ODEs which can be solved using standard methods. These systems of ODEs are much simpler to solve than the systems of DAEs that arise from the DFN model. Additionally, further techniques can be used to speed up the calculations of the solution of the ODE system (e.g residue grouping \cite{Smith2008} or balanced truncation \cite{Jun2015}). The second step (``Potentials'' block in Figure~\ref{fig:block_diagram_EC}b) consists of evaluating the explicit expression for the terminal voltage and for any other variable of interest (such as currents or potentials) which takes the electrolyte and electrode concentrations as inputs. Because this step involves evaluating an explicit expression rather than solving an equation, it is computationally cheap. In addition, we can limit this step to evaluating only the variables we are interested in. The main advantage of this approach is that it fully decouples the problem, so each concentration is the solution of an independent PDE and therefore much cheaper to solve than the DFN model. Note, however, that this decoupling only occurs in the case where current is the input of the model. If instead we want to simulate a potentiostatic discharge, we need to impose an algebraic constraint on the voltage and thus we need to solve a system of DAEs. This system, even though it is more complex than the standard decoupled SPM-type models, is still much simpler than the DFN model.

In the literature, we find several instances of SPM-type models, posed either in an \emph{ad hoc} manner \cite{Moura2017,Bizeray2018,Kemper2013,Perez2016} or a systematic manner \cite{Marquis2019,Richardson2020}. In most cases both approaches yield very similar results, but asymptotic methods are systematic and, therefore, more reliable. Moreover, the systematic approach allows us to easily extend the models to account for extra physical phenomena in a rigorous way, ensuring that underlying physical principles of the DFN model still hold for the SPM-type models. These systematically derived models are obtained by applying asymptotic techniques, in which the size of the dimensionless groupings appearing in the model is used to decide which terms of the equations can be neglected. There are different sets of assumptions that lead to the derivation of SPM-type models (see \cite{Marquis2019,Richardson2020} for two different examples) but they yield very similar results \cite{Timms2020Corrigendum}. However, all of these assumptions are based on the same key idea that allows for SPM-type models: the battery operates in a regime where deviations from the equilibrium potential are small, which means that the battery operates at moderate to low C-rates. At high C-rates the discrepancies between DFN and SPM-type models are considerable, mainly due to the effect of the electrolyte, and thus the DFN model should be used \cite{Moura2017,Marquis2019}. Note that the definition of what is considered a ``high C-rate" for the applicability of these models is not an absolute definition, but it depends on the model parameters.

In the scenario of moderate to low C-rates, the analysis is as follows. If the deviations from equilibrium potential are small, then all the particles within an electrode are approximately at the same potential, which is around the equilibrium (or open-circuit) potential, and therefore they behave very similarly. In consequence, we can define a representative particle (or average particle) for each electrode, and then solve the lithium diffusion problem only for those two particles. Because all the particles behave similarly, the reaction current density for the lithium intercalation must be identical for all particles within an electrode and thus, using a conservation of charge argument, we can calculate it directly from the input current. This is a crucial result because the reaction current density was the condition that coupled the four conservation equations in the DFN model (conservation of mass and charge in both electrode and electrolyte). If the reaction current density is known \emph{a priori}, each of the four problems can now be treated separately and, moreover, the potential equations both in the electrodes and the electrolyte can be integrated analytically. Therefore, the only equations to solve numerically are the ones governing the concentrations and everything else can be calculated from explicit expressions, as shown in Figure~\ref{fig:block_diagram_EC}b. This reasoning can be formalised into a rigorous mathematical analysis using asymptotic techniques, as shown in \ref{sec:reduction_EC_model}. If additional assumptions are taken, there might be further simplifications that take place in the model which may yield different expressions for the potentials (``Potentials'' block in Figure~\ref{fig:block_diagram_EC}b). These further simplifications are discussed in Section \ref{sec:further_simplifications}.


\subsection{Thermal-electrochemical models}\label{sec:TEC_models}
After considering electrochemical models, we focus on the coupled thermal-electrochemical models. These models, apart from the physics discussed in the previous section, describe how the temperature evolves in the battery. The two models are strongly coupled because the parameters of the electrochemical model depend on temperature and the heat generated in the battery depends on the electrochemical processes. This coupling adds an extra layer of complexity to thermal-electrochemical models compared to pure electrochemical models. We can visualise this extra layer of complexity in the form of a block diagram, as shown in Figure~\ref{fig:block_diagram_TEC}. Compared to the electrochemical models in Figure~\ref{fig:block_diagram_EC}, the thermal-electrochemical model has additional blocks to process the thermal inputs and outputs. In particular, the inputs for thermal-electrochemical models are applied current density and ambient temperature ($T_\mathrm{amb}$), and the new outputs are voltage and cell temperature ($T$).

\begin{figure}
    \centering
    \begin{subfigure}{0.48\textwidth}
      \centering
      \includegraphics[scale=1]{img/diagram_TDFN.pdf}  
      \caption{TDFN model}
    \end{subfigure}\hfill
    \begin{subfigure}{0.48\textwidth}
      \centering
      \includegraphics[scale=1]{img/diagram_TSPMe.pdf}  
      \caption{TSPM-type}
    \end{subfigure}
    \caption{Block diagram of thermal-electrochemical models: (a) the thermal Doyler-Fuller-Newman model (TDFN) and (b) the Thermal Single Particle Model type (SPM-type). Each block is a model (or part of a model) that converts the input variable into the output variable (or variables). In brackets, we define the the type of expressions that compose the model (e.g. PDEs, ODEs, explicit expressions) and, if applicable, the type of expressions after spatial discretisation (e.g. DAEs or ODEs). The colour of the blocks indicates the complexity of the model: red for a system of differential-algebraic equations (DAEs, high complexity), yellow for a system of ordinary differential equations (ODEs, medium complexity) and green for explicit expressions (low complexity).}
    \label{fig:block_diagram_TEC}
\end{figure}



The classic example of a thermal-electrochemical model would be the DFN coupled with a thermal model as shown in \cite{Plett2015,Gu2000,Srinivasan2003}. This model is even more computationally expensive than the standard DFN model, and another issue arises when considering a battery composed of multiple layers (we refer to each layer as a ``cell''). As discussed in \cite{Hunt2020}, in this situation the electrochemical behaviour needs to be described at a cell level, but the thermal behaviour needs to be described at a battery level (i.e. across multiple cells). This brings up a multi-scale coupling in the model that requires a hierarchy of submodels, similar to the one between electrodes and particles in the DFN model.

\begin{figure}
    \centering
    \includegraphics[scale=1]{img/multiscale_sketch.pdf}
    \caption{Sketch of the multiple scales involved in a thermal-electrochemical model. The thermal model is posed at the battery level, but it is coupled to the cell level through temperature dependence of parameters and heat generation terms. The cell model, in turn, is coupled with the particle model through the interface reaction.}
    \label{fig:multiscale}
\end{figure}

The thermal model is posed at the battery level because heat can flow from one cell to the neighbouring ones. However, ions cannot flow from one cell to the neighbouring ones, so the electrochemical models should be posed independently for each cell. Even in the case for cylindrical and prismatic cells, composed of big electrodes rolled up around a mandrel, the path for ions to go from one layer to the neighbouring ones is so long compared to the electrode thickness that these effects can be neglected. The two problems are strongly coupled, because the electrochemical properties of each cell can depend on the temperature while the heat generated by the electrochemistry impacts the evolution of temperature. The coupling between these two problems at different scales imposes a hierarchy similar to the classic DFN model. There, we solve for the electrode potential and the electrolyte concentration and potential across the cell, and at each point of the electrode domain there is a representative particle that describes the intercalated lithium concentration at the microscale. Here, we solve for temperature across the battery and at each point there is a representative cell that describes the electrochemical behaviour. Moreover, if the chosen electrochemical model is the DFN, then the hierarchy is at three levels: battery, cell and particle. These relations are shown in Figure~\ref{fig:multiscale}.

The thermal-electrochemical coupling at multiple scales results, after spatial discretisation, into a system of DAEs again, even when using an SPM-type electrochemical model, because for each point at the battery level we have a cell to represent the electrochemical behaviour. These cells are connected in parallel and therefore the total applied current is divided over the different cells, but not in an even way as it depends on the cell temperature. Therefore, we need additional algebraic constraints for the terminal voltage and applied current of each cell. Algebraic constraints were one of the main causes of the complexity of the DFN model, therefore, in this case too, we look for methods to simplify the model. The starting point is the DFN model coupled to a thermal model (see \cite{Plett2015,Gu2000,Srinivasan2003} for details). The simplifications on the electrochemical part of the model follow the procedure explained in Section \ref{sec:EC_models}, so here we focus on the thermal part only. In order to get rid of the DAEs again, we need to work in a regime in which the different instances of the electrochemical model (i.e. the different cells) can be considered to be very similar. It is the same idea as taking that all the particles work similarly within each electrode in the SPM model. There are two main assumptions that lead to this decoupling: the assumption that the temperature variation is small enough to not cause significant variations in the electrochemical parameters, and the assumption that the temperature is approximately homogeneous within the battery. The validity of each assumption depends on the battery chemistry, geometry and cooling conditions. In this work, we focus on the latter assumption as it is reasonable for the batteries and the conditions we are interested in: LG M50 cylindrical batteries cycled in a thermal chamber without any active cooling. The physical reasoning behind the reduction of the thermal model is that in this situation the bottleneck in the heat transfer is the heat dissipation to the environment rather than the heat transport inside the battery. Then, the thermal gradients inside the battery are negligible and the temperature is closely captured by the average temperature and we find that each cell of the battery behaves very similarly so we only need to solve one instance of the electrochemical model. The details for the derivation of the reduced thermal model are shown in \ref{sec:derivation_thermal}.

In terms of the block diagram, the TSPMe works in the following way. A system of coupled PDEs and ODEs (the latter for the temperature) is solved from the input variables, which in this case are the input current and ambient temperature. All the other quantities of interest can then be calculated from them. The main difference compared to the purely electrochemical models is that, if the parameters depend on temperature, the system of PDEs and ODEs is now coupled so concentrations and temperature need to be solved simultaneously (see Figure~\ref{fig:block_diagram_TEC}).


\section{Thermal Single Particle Model with electrolyte (TSPMe)}\label{sec:TSPMe}
Having discussed the background and interpretation of both DFN and SPM-type models, we can now present the Thermal Single Particle Model with electrolyte (TSPMe). The details on the full derivation from the DFN model are given in \ref{sec:derivation_TSPMe}. The TSPMe is a reduced thermal-electrochemical model for lithium-ion batteries. It consists of an SPMe for the electrochemistry coupled with a lumped thermal model for the cell, which requires solving three diffusion equations for the concentrations (two for a representative particle in each electrode and one for the electrolyte) and a first-order ODE for the average cell temperature. This model can be further simplified for some particular cases, as presented in Section \ref{sec:further_simplifications}.

\begin{figure}
    \centering
    \includegraphics[scale=1]{img/geometry_sketch.pdf}
    \caption{Geometry of the TSPMe. Note that we do not pose a geometry for the thermal model because it is a lumped model and therefore there is no space distribution for temperature. The subscripts $\mathrm{n}$ and $\mathrm{p}$ denote the parameters for the negative and positive electrodes, respectively.}
    \label{fig:geometry_sketch}
\end{figure}

The geometry of the model is shown in Figure~\ref{fig:geometry_sketch}. The two electrode particles are defined in the domain $0 \leq r \leq R_k$ (for $k \in \{\mathrm{n},\mathrm{p}\}$), where $\mathrm{n}$ denotes the negative electrode/particle and $\mathrm{p}$ denotes the positive electrode/particle. The electrolyte domain spans across the two porous electrodes and the separator ($0 \leq x \leq L$), where the negative electrode is defined in $0 \leq x \leq L_\mathrm{n}$, the separator in $L_\mathrm{n} \leq x \leq L - L_\mathrm{p}$ and the positive electrode in $L - L_\mathrm{p} \leq x \leq L$. For the full cell, given that the heat equation is averaged over the whole cell, we do not need to define a domain.

We need to solve for the concentrations (and we are interested in certain quantities defined) both in the electrodes and the electrolyte, so we identify them with the subscripts n and p to denote the negative and positive electrodes/particles, respectively, and the subscript e to denote the electrolyte.

Then, in each particle the governing equations for intercalated lithium are 
\begin{subequations}\label{eq:SPMe_cs}
\begin{align}
\pdv{c_{k}}{t} &= \frac{1}{r^2} \pdv{}{r} \left(r^2 D_{k} \pdv{c_{k}}{r} \right), & \quad \text{ in } 0 < r < R_k,\\
\pdv{c_{k}}{r} &= 0, & \quad \text{ at } r = 0,\\
- D_{k} \pdv{c_{k}}{r} &= \frac{J_k}{a_k F}, & \quad \text{ at } r = R_k,\\
c_{k} &= c_{k,\mathrm{init}}, & \quad \text{ at } t = 0,
\end{align}
\end{subequations}
with
\begin{align}
J_\mathrm{n} &= \frac{i_\mathrm{app}}{L_\mathrm{n}}, & J_\mathrm{p} &= - \frac{i_\mathrm{app}}{L_\mathrm{p}},
\end{align}
where $c_{k}$ is the intercalated lithium concentration, $D_{k}$ is the diffusion coefficient of intercalated lithium (which may depend on $c_k$), $R_k$ is the particle radius, $J_k$ is the reaction volumetric current density, $a_k$ is the surface area density, $F$ is the Faraday constant, $c_{k,\mathrm{init}}$ is the initial concentration, $i_\mathrm{app}$ is the applied current density, and $L_k$ is the electrode thickness. Recall that the subscript $k \in \{\mathrm{n},\mathrm{p} \}$ denotes the negative and positive electrode, respectively.

The electrolyte problem is given by
\begin{subequations}\label{eq:SPMe_ce}
\begin{equation}
\varepsilon(x) \pdv{c_{\mathrm{e}}}{t} = \pdv{}{x} \left( D_\mathrm{e} \mathcal{B}(x) \pdv{c_{\mathrm{e}}}{x} \right) + 
\begin{cases}
(1 - t^+) \frac{i_\mathrm{app}}{F L_\mathrm{n}}, & \quad \text{ if } 0 \leq x < L_\mathrm{n},\\
0, & \quad \text{ if } L_\mathrm{n} \leq x < L - L_\mathrm{p},\\
-(1 - t^+) \frac{i_\mathrm{app}}{F L_\mathrm{p}}, & \quad \text{ if } L - L_\mathrm{p} \leq x \leq L,
\end{cases}
\end{equation}
\begin{align}
\pdv{c_{\mathrm{e}}}{x} &= 0, & \quad \text{ at } x = 0, L,\\
c_{\mathrm{e}} &= c_{\mathrm{e},\mathrm{init}}, & \quad \text{ at } t = 0,
\end{align}
\end{subequations}
where $c_\mathrm{e}$ is the lithium-ion concentration in the electrolyte, $\varepsilon(x)$ is the porosity (or electrolyte volume fraction), $D_\mathrm{e}$ is the bulk diffusion coefficient of lithium ions (which may depend on $c_\mathrm{e}$), $\mathcal{B}(x)$ is the geometry factor, $t^+$ is the transference number, $L$ is the total cell thickness, and $c_{\mathrm{e},\mathrm{init}}$ is the initial concentration. 

The spatially dependent porosity  $\varepsilon(x)$ and geometric factor $\mathcal{B}(x)$ are defined as
\begin{align}
\varepsilon(x) &= 
\begin{cases}
\varepsilon_\mathrm{n}, & \quad \text{ if } 0 \leq x < L_\mathrm{n},\\
\varepsilon_\mathrm{s}, & \quad \text{ if } L_\mathrm{n} \leq x < L - L_\mathrm{p},\\
\varepsilon_\mathrm{p}, & \quad \text{ if } L - L_\mathrm{p} \leq x \leq L,
\end{cases} &
\mathcal{B}(x) &= 
\begin{cases}
\mathcal{B}_\mathrm{n}, & \quad \text{ if } 0 \leq x < L_\mathrm{n},\\
\mathcal{B}_\mathrm{s}, & \quad \text{ if } L_\mathrm{n} \leq x < L - L_\mathrm{p},\\
\mathcal{B}_\mathrm{p}, & \quad \text{ if } L - L_\mathrm{p} \leq x \leq L,
\end{cases}
\end{align}
where the values in each subdomain $\varepsilon_k$ and $\mathcal{B}_k$ are constant. The subscripts $k \in \{\mathrm{n}, \mathrm{s}, \mathrm{p} \}$ denote the values for the negative electrode, separator and postive electrode, respectively.

The thermal problem is given by
\begin{subequations}\label{eq:TSPMe_T}
\begin{align}
\theta \dv{T}{t} &= - h a_\mathrm{cool} (T - T_\mathrm{amb}) + Q_\mathrm{s} + Q_\mathrm{e} + Q_\mathrm{irr} + Q_\mathrm{rev},\\
T(0) &= T_\mathrm{amb},
\end{align}
with
\begin{align}
Q_\mathrm{s} &= \frac{i_\mathrm{app}^2}{3 L} \left( \frac{L_\mathrm{n}}{\sigma_{\mathrm{n}}} + \frac{L_\mathrm{p}}{\sigma_{\mathrm{p}}} \right),\\
Q_\mathrm{e} &= - (1-t^+) \frac{2 R T}{F} \frac{i_\mathrm{app}}{L} \left( \frac{1}{L_\mathrm{p}} \int_{L-L_\mathrm{p}}^L \log \left( c_\mathrm{e} \right) \dd x - \frac{1}{L_\mathrm{n}} \int_{0}^{L_\mathrm{n}} \log \left( c_\mathrm{e} \right) \dd x \right) + \frac{1}{L} \int_0^L \frac{i_\mathrm{e}(x,t)^2}{\sigma_\mathrm{e} \left( c_\mathrm{e} (x,t) \right) \mathcal{B}(x)} \dd x,\\
Q_\mathrm{irr} &= \frac{2 R T}{F} \frac{i_\mathrm{app}}{L} \left( \frac{1}{L_\mathrm{n}} \int_0^{L_\mathrm{n}} \arcsinh \left( \frac{i_\mathrm{app}}{a_\mathrm{n} L_\mathrm{n} j_\mathrm{n}}\right) \dd x + \frac{1}{L_\mathrm{p}} \int_{L-L_\mathrm{p}}^{L} \arcsinh \left( \frac{i_\mathrm{app}}{a_\mathrm{p} L_\mathrm{p} j_\mathrm{p}}\right) \dd x \right),\\
Q_\mathrm{rev} &= \frac{i_\mathrm{app}}{L} \left( \Pi_\mathrm{n} - \Pi_\mathrm{p} \right),
\end{align}
\end{subequations}
where $T$ is the average temperature of the cell (assumed to be homogeneous in space), $\theta$ is the lumped volumetric heat capacity of the cell, $h$ is the heat exchange coefficient, $a_\mathrm{cool}$ is the cooling surface area of the cell per unit of volume, $T_\mathrm{amb}$ is the ambient temperature, $Q_\mathrm{s}$ is the heat source term due to the electrode (Joule heating), $Q_\mathrm{e}$ is the heat source term due to the electrolyte (both Joule heating and due to concentration gradients), $Q_\mathrm{irr}$ is the irreversible heat source term of the intercalation reaction and $Q_\mathrm{rev}$ is the reversible heat source term of the intercalation reaction. In addition, $\sigma_{k}$ is the electronic conductivity of the electrodes, $\sigma_\mathrm{e}$ is the ionic conductivity of the electrolyte, $R$ is the universal gas constant, $j_k$ is the reaction surface current density, and $\Pi_k$ is the Peltier coefficient (sometimes referred to as entropic term). This coefficient is defined as $\Pi_k = T \pdv{U_k}{T}$ and describes the heat generated or sunk due to the entropy of the intercalation reaction \cite{Plett2015}.

After solving \eqref{eq:SPMe_cs}, \eqref{eq:SPMe_ce} and \eqref{eq:TSPMe_T} numerically, we can use the values of $c_k$, $c_\mathrm{e}$ and $T$ to calculate any other variable of interest from explicit expressions.


For example, the potentials ($\Phi_\mathrm{n}$, $\Phi_\mathrm{p}$ and $\Phi_\mathrm{e}$) and currents ($i_\mathrm{n}$, $i_\mathrm{p}$ and $i_\mathrm{e}$) in the electrodes and the electrolyte can be calculated as
\begin{subequations}\label{eq:SPMe_Phi}
\begin{multline}
\Phi_{\mathrm{n}} = U_\mathrm{n} (\left. c_{\mathrm{n}} \right|_{r = R_\mathrm{n}}) - \frac{i_\mathrm{app}(2 L_\mathrm{n} - x) x}{2 L_\mathrm{n} \sigma_{\mathrm{n}}} + \frac{i_\mathrm{app} L_\mathrm{n}}{3 \sigma_{\mathrm{n}}} - \frac{1}{L_\mathrm{n}} \int_0^{L_\mathrm{n}} \int_0^x \frac{i_\mathrm{e}(s,t) \dd s}{\sigma_\mathrm{e} \left(c_\mathrm{e}(s,t)\right) \mathcal{B}(s)} \dd x \\
+ (1 - t^+) \frac{2 R T}{F} \frac{1}{L_\mathrm{n}} \int_0^{L_\mathrm{n}} \log \left( \frac{c_\mathrm{e}(x,t)}{c_\mathrm{e}(0,t)}\right) \dd x + \frac{2 R T}{F} \frac{1}{L_\mathrm{n}} \int_0^{L_\mathrm{n}} \arcsinh \left( \frac{i_\mathrm{app}}{a_\mathrm{n} L_\mathrm{n} j_\mathrm{n}} \right) \dd x,\\
\end{multline}
\begin{multline}
\Phi_{\mathrm{p}} = U_\mathrm{p} (\left. c_{\mathrm{p}} \right|_{r = R_\mathrm{p}}) + \frac{i_\mathrm{app}(2 (L - L_\mathrm{p}) - x) x}{2 L_\mathrm{p} \sigma_{\mathrm{p}}} - \frac{i_\mathrm{app} \left(2 L_\mathrm{p}^2 - 6 L L_\mathrm{p} + 3 L^2 \right)}{6 L_\mathrm{p} \sigma_{\mathrm{p}}} - \frac{1}{L_\mathrm{p}} \int_{L - L_\mathrm{p}}^{L} \int_0^x \frac{i_\mathrm{e}(s,t) \dd s}{\sigma_\mathrm{e} \left(c_\mathrm{e}(s,t)\right) \mathcal{B}(s)} \dd x \\
+ (1 - t^+) \frac{2 R T}{F} \frac{1}{L_\mathrm{p}} \int_{L-L_\mathrm{p}}^{L} \log \left( \frac{c_\mathrm{e}(x,t)}{c_\mathrm{e}(0,t)}\right) \dd x - \frac{2 R T}{F} \frac{1}{L_\mathrm{p}} \int_{L - L_\mathrm{p}}^{L} \arcsinh \left( \frac{i_\mathrm{app}}{a_\mathrm{p} L_\mathrm{p} j_\mathrm{p}} \right) \dd x,\\
\end{multline}
\begin{align}
\Phi_{\mathrm{e}} &= (1-t^+)\frac{2 R T}{F} \log \left( \frac{c_\mathrm{e}(x,t)}{c_\mathrm{e}(0,t)} \right) - \int_0^x \frac{i_\mathrm{e}(s,t)}{\sigma_\mathrm{e} \left( c_\mathrm{e} (s,t) \right) \mathcal{B}(s)} \dd s,\\
i_{\mathrm{n}} &= \frac{i_\mathrm{app}}{L_\mathrm{n}}(L_\mathrm{n} - x),\\
i_{\mathrm{p}} &= \frac{i_\mathrm{app}}{L_\mathrm{p}}(x - (L - L_\mathrm{p})),\\
i_\mathrm{e} &= \begin{cases}
\frac{i_\mathrm{app}}{L_\mathrm{n}} x, & \quad \text{ if } 0 \leq x < L_\mathrm{n},\\
i_\mathrm{app}, & \quad \text{ if } L_\mathrm{n} \leq x < L - L_\mathrm{p},\\
\frac{i_\mathrm{app}}{L_\mathrm{p}}(L-x), & \quad \text{ if } L - L_\mathrm{p} \leq x \leq L,
\end{cases}
\end{align}
\end{subequations}
with the exchange current densities $j_\mathrm{n}$ and $j_\mathrm{p}$ defined as
\begin{equation}
j_k = m_k \left. \sqrt{c_\mathrm{e} c_{k} \left(c_{k}^{\max} - c_{k} \right)} \right|_{r = R_k}.
\end{equation}
Here, $U_k$ is the open circuit potential, $m_k$ is the intercalation reaction rate, and $c_{k}^{\max}$ is the maximum concentration in the electrode.

The output voltage is given by
\begin{subequations}\label{eq:V_TSPMe}
\begin{equation}
V_\mathrm{out} = \Phi_{\mathrm{p}}(L,t) - \Phi_{\mathrm{n}}(0,t) = U_\mathrm{eq} + \eta_\mathrm{r} + \eta_\mathrm{e} + \Delta \Phi_\mathrm{e} + \Delta \Phi_\mathrm{s},
\end{equation}
where
\begin{align}
U_\mathrm{eq} &= U_\mathrm{p} (\left.c_{\mathrm{p}}\right|_{r = R_\mathrm{p}}) - U_\mathrm{n} (\left. c_{\mathrm{n}}\right|_{r = R_\mathrm{n}}),\\
\eta_\mathrm{r} &= - \frac{2 R T}{F} \left(\frac{1}{L_\mathrm{p}}\int_{L-L_\mathrm{p}}^{L} \arcsinh \left( \frac{i_\mathrm{app}}{a_\mathrm{p} L_\mathrm{p} j_\mathrm{p}} \right) \dd x + \frac{1}{L_\mathrm{n}} \int_{0}^{L_\mathrm{n}} \arcsinh \left( \frac{i_\mathrm{app}}{a_\mathrm{n} L_\mathrm{n} j_\mathrm{n}} \right) \dd x \right),\\
\eta_\mathrm{e} &= (1-t^+) \frac{2 R T}{F} \left( \frac{1}{L_\mathrm{p}} \int_{L-L_\mathrm{p}}^L \log \left( c_\mathrm{e} \right) \dd x - \frac{1}{L_\mathrm{n}} \int_{0}^{L_\mathrm{n}} \log \left( c_\mathrm{e} \right) \dd x \right),\\
\Delta \Phi_\mathrm{e} &= - \left(\frac{1}{L_\mathrm{p}} \int_{L - L_\mathrm{p}}^{L} \int_0^x \frac{i_\mathrm{e}(s,t) \dd s}{\sigma_\mathrm{e} \left(c_\mathrm{e}(s,t)\right) \mathcal{B}(s)} \dd x - \frac{1}{L_\mathrm{n}} \int_0^{L_\mathrm{n}} \int_0^x \frac{i_\mathrm{e}(s,t) \dd s}{\sigma_\mathrm{e} \left(c_\mathrm{e}(s,t)\right) \mathcal{B}(s)} \dd x \right),\\
\Delta \Phi_\mathrm{s} &= - \frac{i_\mathrm{app}}{3} \left( \frac{L_\mathrm{n}}{\sigma_{\mathrm{n}}} + \frac{L_\mathrm{p}}{\sigma_{\mathrm{p}}} \right).
\end{align}
\end{subequations}
We can interpret these terms in the following way: $U_\mathrm{eq}$ is the equilibrium potential to which the terminal voltage converges when no current is applied, while all the other terms are deviations from equilibrium due to different effects. $\eta_\mathrm{r}$ is due to the reaction overpotentials, $\eta_\mathrm{e}$ is due to the concentration gradients in the electrolyte, $\Delta \Phi_\mathrm{e}$ is due to Ohmic losses in the electrolyte, and $\Delta \Phi_\mathrm{s}$ is due to Ohmic losses in the (solid) electrode.

The TSPMe is valid in a particular regime both in terms of electrochemical and thermal submodels, and the physical intuition for them is discussed in Section \ref{sec:SPM_discussion}. The electrochemical submodel is valid when the deviations from the equilibrium potential are slow, which is equivalent to moderate to low C-rates. In terms of the dimensionless parameters, this corresponds to the case when
\begin{align}
    \lambda &= \frac{F \Phi_0}{R T_\mathrm{amb}} \gg 1, & \Sigma_k &= \frac{R T_\mathrm{amb}}{F L i_0} \sigma_k \gtrsim 1, & \Sigma_\mathrm{e} &= \frac{R T_\mathrm{amb}}{F L i_0} \sigma_{\mathrm{e},\mathrm{typ}} \gtrsim 1,
\end{align}
where $i_0$ and $\Phi_0$ are the typical current and electrode potential, respectively, and $\sigma_{\mathrm{e},\mathrm{typ}}$ is the typical value of $\sigma_\mathrm{e}$. The assumption $\lambda \gg 1$ (i.e. $\lambda$ significantly larger than one) means that $\lambda$ must be larger than one for the assumption to hold and, the larger $\lambda$ is, the better the reduced model works. On the other hand, the assumption of $\Sigma_k \gtrsim 1$ (and similarly for $\Sigma_\mathrm{e}$) means that $\Sigma_k$ must be of $\order{1}$ or larger (i.e. not significantly smaller than one). Further details on these dimensionless parameters, as well as the details of the asymptotic reduction, can be found in \ref{sec:derivation_TSPMe}.

The thermal submodel is valid when the heat exchange with environment is much slower than the heat transfer inside of the battery. In terms of dimensionless parameters, this corresponds to the case when
\begin{align}
    \mathrm{Bi} &= \frac{h L_\mathrm{batt}}{\kappa} \ll 1, & \mathcal{K} &= \frac{\kappa t_0}{L^2_\mathrm{batt} \theta} \gg 1,
\end{align}
where $L_\mathrm{batt}$ is the typical length scale of the battery, $\kappa$ is the lumped thermal conductivity of the battery and $t_0$ is the typical discharge time.

To the best of our knowledge, the TSPMe presented here is the first instance of a formal derivation (i.e. using asymptotic techniques) of a TSPM-type model for multi-layer batteries and including finite diffusion effects in the particles. There are other instances of asymptotically derived thermal-electrochemical models in the literature, but they focus on other aspects. The work in \cite{Marquis2020,Timms2020} focuses on single layer pouch cells. On the other hand, the analysis in \cite{Hennessy2019} assumes infinitely fast diffusion, which does capture the relaxation of voltage that we observe in the experimental data when the battery is switched off (see Section \ref{sec:results}). 

In terms of the electrochemical model, our analysis is based on the same assumptions as the work by Richardson et al. \cite{Richardson2020}, but in our case we focus on the homogeneous electrodes case. Therefore, we are able to write the electrolyte potential (\ref{eq:SPMe_Phi}c) as an explicit expression rather than having to solve an elliptic PDE for it. This allows us to write simpler expressions for all the potentials, while retaining all the features in \cite{Richardson2020} for homogeneous electrodes. On the other hand, compared to the model by Marquis et al. \cite{Marquis2019}, our approach directly captures all the nonlinearities in the voltage expression, so there is no need to introduce them \emph{a posteriori} (see \cite{Marquis2020} for details).


\subsection{Further simplifications}\label{sec:further_simplifications}
Notice that the TSPMe is already an SPM-type model, as defined in Section \ref{sec:SPM_discussion}, which can be implemented at a very low computational cost. Therefore, from the point of view of computational efficiency, there is no need to further simplify the model. However, in the literature we find further simplifications to SPM-type models. These further simplified models can also be helpful to understand the physical behaviour of batteries and for applications where the available computational power is really low (such as in battery management systems). In this section we present some further simplifications that can be done by taking additional assumptions to show the links between the TSPMe model and other models found in the literature \cite{Marquis2019,Richardson2020,Marquis2020,Richardson2012,Moyles2019asymptotic,Moyles2018experimental,Hennessy2019}. The full details of the derivation can be found in \ref{sec:further_simplifications_app}.

Here we present the dimensional form of three main different simplifications: quasi-steady-state electrolyte concentration, constant electrolyte conductivity and fast lithium diffusion. For the first two we consider as well a particular case of each that allows even further reduction of the model (small variation of electrolyte concentration and large electrolyte conductivity). Notice that the three simplifications are independent from each other, so we can choose to use a subset of them or all together.

\subsubsection{Quasi-steady state electrolyte concentration}
One simplification for the electrolyte concentration is to take the quasi-steady-state problem, which is a valid assumption when the current varies over a larger time scale than the electrolyte diffusion and electrolyte ion generation time scales. Mathematically, this corresponds to the case when
\begin{align}
    \mathcal{C}_\mathrm{e} &= \frac{L^2}{D_{\mathrm{e},\mathrm{typ}} t_0} \ll 1, & \gamma_\mathrm{e} &= \frac{c_{\mathrm{e},\mathrm{init}}}{c_\mathrm{n}^{\max}} \ll 1,
\end{align}
where $D_{\mathrm{e},\mathrm{typ}}$ is the typical value of $D_\mathrm{e}$ and $c_\mathrm{n}^{\max}$ is the maximum concentration in the negative electrode particles.

The physical meaning of this limit is that electrolyte transient effects are negligible because they happen at a much shorter time scale than variations in the applied current. This is a reasonable assumption for constant-current discharge (or charge) experiments.

With these assumptions, there is no time dependence in the electrolyte concentration equations so their integration is a one-off step. Additionally, if the ion transport properties in the electrolyte ($D_\mathrm{e}$ and $t^+$) are assumed to be constant, we can find the following analytical expression for the concentration
\begin{subequations}\label{eq:ce_ss}
\begin{equation}
    c_\mathrm{e} = c_{\mathrm{e},\mathrm{init}} + \frac{i_\mathrm{app} (1 - t^+)}{6 F D_\mathrm{e} v_\mathrm{pore}} \Delta c (x)
\end{equation}
with
\begin{equation}
    \Delta c (x) = \begin{cases}
        \left(\frac{2 \varepsilon_\mathrm{p} L_\mathrm{p}^2}{\mathcal{B}_\mathrm{p}} + \frac{3 L_\mathrm{s} (2 \varepsilon_\mathrm{p} L_\mathrm{p} + \varepsilon_\mathrm{s} L_\mathrm{s})}{\mathcal{B}_\mathrm{s}} - \frac{3 \frac{v_\mathrm{pore}}{L_\mathrm{n}} (L_\mathrm{n}^2 - x^2) - 2 \varepsilon_\mathrm{n} L_\mathrm{n}^2}{\mathcal{B}_\mathrm{n}} \right), & \text{ if } 0 < x < L_\mathrm{n},\\
        \left( - \frac{2 \varepsilon_\mathrm{n} L_\mathrm{n}^2}{\mathcal{B}_\mathrm{n}} + \frac{2 \varepsilon_\mathrm{p} L_\mathrm{p}^2}{\mathcal{B}_\mathrm{p}} + \frac{6 \frac{v_\mathrm{pore}}{L} (L - L_\mathrm{p} - x) - 6 \varepsilon_n L_\mathrm{n} L_\mathrm{s} - 3 \varepsilon_\mathrm{s} L_\mathrm{s}^2}{\mathcal{B}_\mathrm{s}}\right), & \text{ if } L_\mathrm{n} < x < L - L_\mathrm{p},\\
        \left( - \frac{2 \varepsilon_\mathrm{n} L_\mathrm{n}^2}{\mathcal{B}_\mathrm{n}} - \frac{3 L_\mathrm{s} (2 \varepsilon_\mathrm{n} L_\mathrm{n} + \varepsilon_\mathrm{s} L_\mathrm{s})}{\mathcal{B}_\mathrm{s}} + \frac{3 \frac{v_\mathrm{pore}}{L_\mathrm{p}} \left( (L^2 - L_\mathrm{p}^2) - (2L -x) x \right) + 2 \varepsilon_\mathrm{p} L_\mathrm{p}^2}{\mathcal{B}_\mathrm{p}}\right), & \text{ if } L - L_\mathrm{p} < x < L,
    \end{cases}
\end{equation}
\end{subequations}
where the volume of pores per unit of electrode plate area is given by
\begin{equation}
v_\mathrm{pore} = \int_0^L \varepsilon(x) \dd x  = \varepsilon_\mathrm{n} L_\mathrm{n} + \varepsilon_\mathrm{s} L_\mathrm{s} + \varepsilon_\mathrm{p} L_\mathrm{p},
\end{equation}
and $L_\mathrm{s} = L - L_\mathrm{n} - L_\mathrm{p}$ is the thickness of the separator.

An additional step to the previous simplification occurs when the variation in electrolyte concentration is small. This corresponds to the limit
\begin{align}
    \mathcal{C}_\mathrm{e} &= \frac{L^2}{D_{\mathrm{e},\mathrm{typ}} t_0} \ll 1, & \gamma_\mathrm{e} &= \frac{c_{\mathrm{e},\mathrm{init}}}{c_\mathrm{n}^{\max}} \sim 1,
\end{align}
and it yields the result $c_\mathrm{e} = c_{\mathrm{e},\mathrm{init}}$. Therefore, we obtain a Single Particle Model as we do not need to solve for $c_{\mathrm{e}}$. Notice that this limit is the same to the one studied in \cite{Marquis2019} in terms of electrolyte concentration. The only differences between that model and the one presented here are in the voltage expressions. In their case they assumed large conductivity in the electrodes, while we assume small overpotentials (i.e. $\lambda \gg 1$), but the results are analogous.  

This assumption also helps simplify the expression for the electrolyte potential because we find that $\sigma_\mathrm{e} (c_\mathrm{e})$ is a constant, which is the case discussed in the next section.

\subsubsection{Constant electrolyte conductivity}\label{sec:constant_sigmae}
One assumption to simplify the electrolyte potential is to take that the electrolyte conductivity is a constant. Then, the integrals involving $\sigma_\mathrm{e}$ in (\ref{eq:TSPMe_T}c), (\ref{eq:SPMe_Phi}) and (\ref{eq:V_TSPMe}) can be computed analytically. This means a slight reduction in the computational cost, which is negligible for a laptop but can be important for simpler devices such as the processors in the battery management systems of the electric vehicles. 

Hence, the potentials in both electrode and electrolyte can be calculated using
\begin{subequations}
\begin{multline}
\Phi_{\mathrm{n}} = U_\mathrm{n} (\left.c_{\mathrm{n}}\right|_{r = R_\mathrm{n}}) - \frac{i_\mathrm{app}(2 L_\mathrm{n} - x) x}{2 L_\mathrm{n} \sigma_{\mathrm{n}}} + \frac{i_\mathrm{app} L_\mathrm{n}}{3 \sigma_{\mathrm{n}}} - \frac{i_\mathrm{app}}{6 \sigma_\mathrm{e}} \frac{L_\mathrm{n}}{\mathcal{B}_\mathrm{n}} \\
+ (1 - t^+) \frac{2 R T}{F} \frac{1}{L_\mathrm{n}} \int_0^{L_\mathrm{n}} \log \left( \frac{c_\mathrm{e}(x,t)}{c_\mathrm{e}(0,t)}\right) \dd x + \frac{2 R T}{F} \frac{1}{L_\mathrm{n}} \int_0^{L_\mathrm{n}} \arcsinh \left( \frac{i_\mathrm{app}}{a_\mathrm{n} L_\mathrm{n} j_\mathrm{n}} \right) \dd x,
\end{multline}
\begin{multline}
\Phi_{\mathrm{p}} = U_\mathrm{p} (\left.c_{\mathrm{p}}\right|_{r = R_\mathrm{p}}) + \frac{i_\mathrm{app}(2 (L - L_\mathrm{p}) - x) x}{2 L_\mathrm{p} \sigma_{\mathrm{p}}} - \frac{i_\mathrm{app} \left(2 L_\mathrm{p}^2 - 6 L L_\mathrm{p} + 3 L^2 \right)}{6 L_\mathrm{p} \sigma_{\mathrm{p}}} - \frac{i_\mathrm{app}}{6 \sigma_\mathrm{e}} \left(\frac{3 L_\mathrm{n}}{\mathcal{B}_\mathrm{n}} + \frac{6 L_\mathrm{s}}{\mathcal{B}_\mathrm{s}} + \frac{2 L_\mathrm{p}}{\mathcal{B}_\mathrm{p}} \right)\\
+ (1 - t^+) \frac{2 R T}{F} \frac{1}{L_\mathrm{p}} \int_{L-L_\mathrm{p}}^{L} \log \left( \frac{c_\mathrm{e}(x,t)}{c_\mathrm{e}(0,t)}\right) \dd x - \frac{2 R T}{F} \frac{1}{L_\mathrm{p}} \int_{L - L_\mathrm{p}}^{L} \arcsinh \left( \frac{i_\mathrm{app}}{a_\mathrm{p} L_\mathrm{p} j_\mathrm{p}} \right) \dd x,
\end{multline}
\begin{align}
\Phi_{\mathrm{e}} &= (1-t^+)\frac{2 R T}{F} \log \left( \frac{c_\mathrm{e}(x,t)}{c_\mathrm{e}(0,t)} \right) + \frac{i_\mathrm{app}}{2 \sigma_\mathrm{e}}
\begin{cases}
-\frac{x^2}{\mathcal{B}_\mathrm{n} L_\mathrm{n}}, & \quad \text{ if } 0 \leq x < L_\mathrm{n},\\
-\frac{2 (x - L_\mathrm{n})}{\mathcal{B}_\mathrm{s}} - \frac{L_\mathrm{n}}{\mathcal{B}_\mathrm{n}}, & \quad \text{ if } L_\mathrm{n} \leq x < L - L_\mathrm{p},\\
\frac{(L-x)^2}{\mathcal{B}_\mathrm{p} L_\mathrm{p}} - \left(\frac{L_\mathrm{n}}{\mathcal{B}_\mathrm{n}} + \frac{2 L_\mathrm{s}}{\mathcal{B}_\mathrm{s}} + \frac{L_\mathrm{p}}{\mathcal{B}_\mathrm{p}} \right), & \quad \text{ if } L - L_\mathrm{p} \leq x \leq L.
\end{cases}
\end{align}
\end{subequations}

Then, the Ohmic losses in the electrolyte can be written as
\begin{equation}
\Delta \Phi_\mathrm{e} = - \frac{i_\mathrm{app}}{3 \sigma_\mathrm{e}} \left( \frac{L_\mathrm{n}}{\mathcal{B}_\mathrm{n}} + \frac{3 L_\mathrm{s}}{\mathcal{B}_\mathrm{s}} + \frac{L_\mathrm{p}}{\mathcal{B}_\mathrm{p}} \right),
\end{equation}
and the electrolyte heat source term can be written as
\begin{equation}
Q_\mathrm{e} = - (1-t^+) \frac{2 R T}{F} \frac{i_\mathrm{app}}{L} \left( \frac{1}{L_\mathrm{p}} \int_{L-L_\mathrm{p}}^L \log \left( c_\mathrm{e} \right) \dd x - \frac{1}{L_\mathrm{n}} \int_{0}^{L_\mathrm{n}} \log \left( c_\mathrm{e} \right) \dd x \right) + \frac{i_\mathrm{app}^2}{3 L \sigma_\mathrm{e}} \left( \frac{L_\mathrm{n}}{\mathcal{B}_\mathrm{n}} + \frac{3 L_\mathrm{s}}{\mathcal{B}_\mathrm{s}} + \frac{L_\mathrm{p}}{\mathcal{B}_\mathrm{p}} \right).
\end{equation}

Note that, as mentioned in the previous section, these expressions for the potentials are very similar to those obtained in the limit where the variation of electrolyte concentration is small. Therefore, we can observe that the Ohmic losses in the electrolyte have the same form as those in \cite{Marquis2019}.

If the electrolyte conductivity is assumed to be high, then an additional simplification can be done by eliminating all the terms involving $\sigma_\mathrm{e}^{-1}$. Mathematically, this corresponds to the limit
\begin{equation}
    \Sigma_\mathrm{e} = \frac{R T_\mathrm{amb}}{F L i_0} \sigma_{\mathrm{e},\mathrm{typ}} \gg 1.
\end{equation}

\subsubsection{Fast electrode diffusion}\label{sec:fast_diffusion}
Finally, the last simplification that can be taken is fast electrode diffusion, which corresponds to the limit
\begin{equation}
    \mathcal{C}_k = \frac{R^2_k}{D_{k,\mathrm{typ}} t_0} \ll 1.
\end{equation}
In this case, the particle problem \eqref{eq:SPMe_cs} reduces to an ODE of the form
\begin{subequations}
\begin{align}
    \dv{c_k}{t} &= -\frac{J_k}{a_k R_k F}, \\
    c_k(0) &= c_{k,\mathrm{init}},
\end{align}
\end{subequations}
and we find several instances of this model in the literature, such as \cite{Richardson2012,Moyles2019asymptotic,Moyles2018experimental,Hennessy2019}. Note that this fast electrode diffusion model does not capture the relaxation effects that we observe when the current is switched off. These effects, mostly noticeable in the terminal voltage, are caused by the relaxation of the particle concentration to its steady state and therefore require a finite diffusion in the particles.

\section{Results}\label{sec:results}
After presenting the TSPMe model, we now validate it to show its relevance and applicability. To show the accuracy and performance of the reduced model, we start by comparing the TSPMe with the thermal DFN (TDFN). Then, to show the applicability of the TSPMe in real applications we compare it against experimental data on the LG M50 cell. This is a 21700 cylindrical cell, with a nominal capacity of 5 Ah, and with an NMC811 positive electrode and a graphite and SiO$_x$ negative electrode (see \cite{Chen2020} for further details). In both cases, we perform comparisons at different C-rates (C/2, 1C and 2C) and ambient temperatures (0~$^\circ \mathrm{C}$, 10~$^\circ \mathrm{C}$ and 25~$^\circ \mathrm{C}$). The code to reproduce the simulations is publicly available in an online repository (\url{https://github.com/brosaplanella/TEC-reduced-model}).

\subsection{Comparison between TSPMe and DFN}
We first look at the comparison between the reduced model (the TSPMe) and a thermal DFN model (using a lumped thermal model, see \cite{Marquis2020}) under different discharge conditions to assess how well does the TSPMe approximate the TDFN. Here we use a lumped thermal model coupled to the DFN for simplicity, as resolving the battery geometry would require a way heavier computational approach that is out of the scope of this work. This increase of computational cost is because, as detailed in \cite{Tranter2020}, we would have to define a model over three scales (battery, cell and particle, see Figure \ref{fig:multiscale}) with the corresponding battery geometry, which would also be challenging from the implementation point of view. The price we pay to use a lumped model for the TDFN is that we do not capture the thermal gradients in the battery. Validating the TSPMe model against a spatially resolved TDFN (like \cite{Tranter2020}) is an area for future work on this topic.

The models have been implemented in PyBaMM (Python Battery Mathematical Modelling) package, a Python-based open-source package to simulate battery models \cite{Sulzer2020}. The code and data are available online (see ``Data and code availability''). Both models are solved using a finite volumes method for the spatial discretisation using the same number of points in both methods for a given domain. For the results here, we used 30 points in the electrode particles and 20 points for each electrode and separator. With this number of points in the discretisation the model yields accurate results, and we do not observe a significant improvement in the solution when increasing the number of points in the mesh. The code to compare the different mesh sizes can be found in the online repository.

After spatial discretisation, we obtain a system of 122 ODEs for the TSPMe model and a system of 1262 ODEs and 100 algebraic equations for the TDFN model, which is an index one DAE \cite{Marquis2019}. Observe that, for the same mesh size, the TDFN model requires solving a system over 10 times larger than the TSPMe, and much more complex due to the algebraic constraints. In order to solve the models we used the SciPy ODE solver \cite{Virtanen2020} for the TSPMe, while for the TDFN we needed a DAE solver, so we used the CasADI solver \cite{Andersson2019}.

\begin{table
\centering
\begin{tabular}{| c c l c c c |}
\hline
\textbf{Symbol} & \textbf{Units} & \textbf{Description} & \textbf{Pos.} & \textbf{Sep.} & \textbf{Neg.} \\ \hline
$L_{k}$ & m & Thickness & $75.6\E{-6}$ & $12\E{-6}$ & $85.2\E{-6}$ \\
$R_{k}$ & m & Radius of electrode particles & $5.22\E{-6}$ & - & $5.86\E{-6}$ \\
$a_{k}$ & m$^{-1}$ & Particle surface area density & $3.82\E{5}$ & - & $3.84\E{5}$ \\
$D_{k}$ & $\mathrm{m}^2 \; \mathrm{s}^{-1}$ & Lithium diffusivity in particles & $4\E{-15}$ & - & $3.3\E{-14}$ \\
$\sigma_{k}$ & $\mathrm{S} \; \mathrm{m}^{-1}$ & Electrode conductivity & 0.18 & - & 215 \\
$c_{k,\mathrm{init}}$ & $\mathrm{mol} \; \mathrm{m}^{-3}$ & Initial particle concentration & 17038 & - & 29866 \\
$c_{k}^{\max}$ & $\mathrm{mol} \; \mathrm{m}^{-3}$ & Max. particle concentration & 63104 & - & 33133 \\
$U_k$ & V & Open circuit potential & Fig. \ref{fig:OCVs} & - & Fig. \ref{fig:OCVs} \\
$m_k$ & $\mathrm{A} \; \mathrm{m}^{-2} \left(\mathrm{mol} \; \mathrm{m}^{-3}\right)^{-1.5}$ & Reaction rate & $3.42\E{-6}$ & - & $6.48\E{-7}$ \\ \hline
$\varepsilon_{k}$ & - & Electrolyte volume fraction & 0.335 & 0.47 & 0.25 \\
$D_{\mathrm{e}}$ & $\mathrm{m}^2 \; \mathrm{s}^{-1}$ & Electrolyte diffusivity & \multicolumn{3}{c |}{see \cite{Chen2020}} \\
$\sigma_{\mathrm{e}}$ & $\mathrm{S} \; \mathrm{m}^{-1}$ & Electrolyte conductivity & \multicolumn{3}{c |}{see \cite{Chen2020}} \\
$t^+$ & - & Transfer number & \multicolumn{3}{c |}{0.2594} \\
$c_{\mathrm{e},\mathrm{init}}$ & $\mathrm{mol} \; \mathrm{m}^{-3}$ & Initial electrolyte concentration & \multicolumn{3}{c |}{1000} \\ \hline
$i_{\mathrm{app}}$ & $\mathrm{A} \; \mathrm{m}^{-2}$ & Applied current density & \multicolumn{3}{c |}{$48.69 C$}\\
$F$ & $\mathrm{C} \; \mathrm{mol}^{-1}$ & Faraday constant & \multicolumn{3}{c |}{96485}\\
$R$ & $\mathrm{J} \; \mathrm{K}^{-1} \; \mathrm{mol}^{-1}$ & Gas constant & \multicolumn{3}{c |}{8.314}\\
$a_\mathrm{cool}$ & $\mathrm{m}^{-1}$ & Cooling surface area density & \multicolumn{3}{c |}{219.42} \\
$L_\mathrm{batt}$ & m & Length scale of the battery* & \multicolumn{3}{c |}{$1\E{-2}$}\\ 
$T_\mathrm{amb}$ & K & Ambient temperature & \multicolumn{3}{c |}{298}\\ \hline
$\theta$ & $\mathrm{J} \; \mathrm{K}^{-1} \; \mathrm{m}^{-3}$ & Volumetric heat capacity & \multicolumn{3}{c |}{$2.85\E{6}$} \\
$\kappa$ & $\mathrm{W} \; \mathrm{m}^{-1} \; \mathrm{K}^{-1}$ & Thermal conductivity* & \multicolumn{3}{c |}{1.05} \\
$h$ & $\mathrm{W} \; \mathrm{m}^{-2} \; \mathrm{K}^{-1}$ & Heat exchange coefficient & \multicolumn{3}{c |}{20} \\ \hline
\end{tabular}
\caption{Dimensional parameters for the LG M50 cell. The electrochemical parameters are taken from \cite{Chen2020} and the thermal parameters are taken from \cite{Taheri2013}. The parameters marked with an asterisk indicate that they are not directly used in the TSPMe but have been used to evaluate the dimensionless grouping for its derivation (see \ref{sec:derivation_TSPMe}). The $C$ in the definition of $i_\mathrm{app}$ corresponds to the C-rate of the experiment. The open-circuit potentials of the electrodes are taken from \cite{Chen2020} as well, and they are plotted in Figure \ref{fig:OCVs}.}
\label{tab:parameter_values_LGM50}
\end{table}

For the simulations we used the parameter values shown in Table~\ref{tab:parameter_values_LGM50}. The electrochemical parameters are for the LG M50 cell  and have been taken from \cite{Chen2020}, while the thermal parameters ($\theta$, $\kappa$ and $h$) have been taken from \cite{Taheri2013}. For the thermal parameters, given that they are lumped parameters for the whole battery, we took from \cite{Taheri2013} the values for each part and calculated the lumped value using the thicknesses from \cite{Chen2020}. For the heat exchange coefficient, we took the intermediate value from the range discussed in \cite{Taheri2013}. Note that two of the parameters in Table~\ref{tab:parameter_values_LGM50} ($L_\mathrm{batt}$ and $\kappa$) do not directly appear in the TSPMe model, but their values have been used to calculate the dimensionless parameters in Table~\ref{tab:ND_parameter_values_LGM50}, which are needed in the asymptotic analysis presented in \ref{sec:derivation_TSPMe}. These parameters are identified with an asterisk in Table~\ref{tab:parameter_values_LGM50}. To capture the microstructure effects on the macroscopic transport, we use the Bruggeman correlation, so we set $\mathcal{B}_k = \varepsilon_k^{1.5}$. We also do not include the reversible heat generation effects $Q_\mathrm{rev}$ because the Peltier coefficient $\Pi_k = T \pdv{U_k}{T}$ needs to be consistent with the open-circuit potential $U_k$ and this would require a careful and exhaustive experimental analysis that is out of the scope of this work. This is an aspect that will be addressed in future work.

The comparison between the TDFN model and TSPMe is shown in Figures~\ref{fig:comp_models_25degC}-\ref{fig:comp_models_0degC}, with the error data in Table~\ref{tab:error_models} and the computational times in Table~\ref{tab:comp_time}. We observe that the TSPMe does a very good job representing the voltage and the cell temperature compared to the TDFN model, and that the discrepancy increases with the C-rate. This could be expected as, from the asymptotic analysis in \ref{sec:derivation_TSPMe}, we know that some of the assumptions on the parameter sizes break down if the applied current is too large. The data in Table~\ref{tab:error_models} shows a slight trend that the agreement between models improves as the temperature decreases, which could be caused because lower temperatures yield lower values of the thermal potential, which is one of the assumptions in the asymptotic analysis. 

\begin{figure
    \centering
    \includegraphics[scale=1]{img/comp_models_25degC.png}
    \caption{Comparison between the TDFN and TSPMe solutions at 25 $^\circ \mathrm{C}$ and different C-rates. The plots on the left column compare the terminal voltage and cell temperature, where the solid colour lines represent the TSPMe and the black dashed lines represent the TDFN. The plots on the right column show the absolute error between the two models for voltage and temperature, respectively.}
    \label{fig:comp_models_25degC}
\end{figure}


\begin{figure
    \centering
    \includegraphics[scale=1]{img/comp_models_10degC.png}
    \caption{Comparison between the TDFN and TSPMe solutions at 10 $^\circ \mathrm{C}$ and different C-rates. The plots on the left column compare the terminal voltage and cell temperature, where the solid colour lines represent the TSPMe and the black dashed lines represent the TDFN. The plots on the right column show the absolute error between the two models for voltage and temperature, respectively.}
    \label{fig:comp_models_10degC}
\end{figure}


\begin{figure
    \centering
    \includegraphics[scale=1]{img/comp_models_0degC.png}
    \caption{Comparison between the TDFN and TSPMe solutions at 0 $^\circ \mathrm{C}$ and different C-rates. The plots on the left column compare the terminal voltage and cell temperature, where the solid colour lines represent the TSPMe and the black dashed lines represent the TDFN. The plots on the right column show the absolute error between the two models for voltage and temperature, respectively.}
    \label{fig:comp_models_0degC}
\end{figure}

Meanwhile, the computational time to solve each model shows that the TSPMe is between approximately 10 times (2C) and 40 times (C/2) faster than the TDFN model. We also notice that, for the TSPMe, the computational time does not vary much with the C-rate, but for the TDFN model the computational time is inversely proportional to the C-rate. This is due to the way PyBaMM handles solver events (e.g. stopping the solver when the critical discharge voltage is reached). The parameters of the solver could be changed to reduce the solving time of the TDFN, but they would need to be adjusted in a case by case basis given that the solver is very sensitive to the parameter values and the experimental conditions. On the other hand, there is no significant effect of the experiment temperature on the computational time. The computational times have been calculated on a laptop with an Intel Core i7-7660U (2.50Ghz) processor and 16 GB RAM. To eliminate the artefacts caused by other background processes, each model has been solved 20 times and the values shown in Table~\ref{tab:comp_time} are the mean and standard deviation of the different runs for each C-rate and temperature. 

This huge difference in the computational speed is because the PDEs arising in the TSPMe are discretised into a system of ODEs which is well-behaved, while the discretised TDFN model gives a DAE system which is ill-posed and, therefore, solving this system requires very specific numerical schemes that are computationally expensive \cite{Gerdts2015}. Specialised techniques to speed up solving both systems of equations are available (like the reduction methods of residue grouping \cite{Smith2008} or balanced truncation \cite{Jun2015} for systems of ODEs). However, given that the TSPMe is intrinsically simpler than the TDFN, under equivalent conditions the former should always be faster. From this comparison, we conclude that the TSPMe provides a very good approximation to the thermal DFN model at a much lower computational cost. In the next section we assess if this model can be used to obtain meaningful predictions of experimental data.

\begin{table
    \begin{subtable}{.5\linewidth}
      \centering
        \begin{tabular}{r | c | c | c}
            & 0.5C & 1C & 2C \\ \hline
            25 $^\circ \mathrm{C}$ & 2.10 (5.87) & 5.59 (16.35) & 23.95 (63.61) \\
            10 $^\circ \mathrm{C}$ & 1.72 (5.10) & 4.97 (14.62) & 22.58 (60.71) \\
            0 $^\circ \mathrm{C}$ & 1.64 (4.98) & 4.82 (14.05) & 22.10 (59.15) \\
        \end{tabular}
        \caption{Voltage (mV)}
    \end{subtable}%
    \begin{subtable}{.5\linewidth}
      \centering
        \begin{tabular}{r | c | c | c}
            & 0.5C & 1C & 2C \\ \hline
            25 $^\circ \mathrm{C}$ & 0.03 (0.05) & 0.15 (0.29) & 1.14 (1.92) \\
            10 $^\circ \mathrm{C}$ & 0.02 (0.04) & 0.13 (0.24) & 1.07 (1.75) \\
            0 $^\circ \mathrm{C}$ & 0.02 (0.03) & 0.13 (0.23) & 1.04 (1.70) \\
        \end{tabular}
        \caption{Temperature ($^\circ \mathrm{C}$)}
    \end{subtable}
    \caption{Error between TSPMe and DFN models for different temperatures and C-rates. The first values are the root-mean-squared-error (RMSE) and the second values (in brackets) are the peak error.}
    \label{tab:error_models}    
\end{table}


\begin{table
    \begin{subtable}{.5\linewidth}
      \centering
        \begin{tabular}{r | c | c | c}
            & 0.5C & 1C & 2C \\ \hline
            25 $^\circ \mathrm{C}$ & $0.44 \pm 0.02$ & $0.43 \pm 0.03$ & $0.47 \pm 0.03$ \\
            10 $^\circ \mathrm{C}$ & $0.44 \pm 0.02$ & $0.44 \pm 0.04$ & $0.47 \pm 0.03$ \\
            0 $^\circ \mathrm{C}$ & $0.44 \pm 0.02$ & $0.43 \pm 0.03$ & $0.47 \pm 0.05$ \\
        \end{tabular}
        \caption{TSPMe}
    \end{subtable}%
    \begin{subtable}{.5\linewidth}
      \centering
        \begin{tabular}{r | c | c | c}
            & 0.5C & 1C & 2C \\ \hline
            25 $^\circ \mathrm{C}$ & $19.03 \pm 0.19$ & $9.25 \pm 0.23$ & $9.02 \pm 0.14$ \\
            10 $^\circ \mathrm{C}$ & $18.81 \pm 0.06$ & $9.01 \pm 0.04$ & $8.44 \pm 0.04$ \\
            0 $^\circ \mathrm{C}$ & $18.64 \pm 0.13$ & $9.09 \pm 0.09$ & $8.39 \pm 0.06$ \\
        \end{tabular}
        \caption{TDFN}
    \end{subtable} 
    \caption{Computational time (in seconds) for the TSPMe and the TDFN model. The times shown here are the solving times only (they do not include the system assembly time). To obtain reliable measurements each simulation has been run 20 times and the values shown here are the mean and standard deviation of the samples.}
    \label{tab:comp_time}
\end{table}


\subsection{Comparison with experimental data}
Now we compare the TSPMe model with experimental data to assess if the model can be used in real applications. The experimental data is for a commercial cell, the LG M50. For the rate tests, cells were placed inside an Espec thermal chamber and, for accuracy, thermocouples were used to record the temperature inside the chamber. Cell temperature was monitored on the cell mid surface and four cells were used for each C-rate and ambient temperature. One of these cells had two additional thermocouples: one on the positive tab and one on the negative tab. A 10 A Digatron battery cycler was used for the tests. For the temperature measurements, K-type thermocouples from RS components were used, which is accurate up to a standard deviation of $\pm 0.75 \%$. Three temperature settings 0 $^\circ \mathrm{C}$, 10 $^\circ \mathrm{C}$ and 25 $^\circ \mathrm{C}$ were used and for each cycle. The cycle consisted of a constant current (CC) charge step of C/3 to a cut-off voltage of 4.2 V, followed by a constant voltage (CV) step at 4.2 V until the charge current dropped to C/20, then a two-hour long rest period, before a constant current discharge step to 2.5 V and a final two-hour long rest before repeating it but for a different discharge C-rate. Different discharge C-rates were explored: C/10, C/2, 1C and 2C, running two cycles for each discharge rate. To reduce the data size, a variable time step combined with a voltage difference condition was used. For C/10 and C/2 discharge rates we used 60 s or 10 mV sampling and for 1C and 2C discharge rates we used 1 s or 10 mV sampling. Data was recorded once one of the two conditions was satisfied. Note that one of the cells cycled at 1C (cell 791) has not been included in the validation data because it gave faulty temperature measurements. However, the data for that cell is available in the online repository. We would like to highlight here that the purpose of this work is not to parameterise the thermal model for a specific cell, but to derive and validate a reduced thermal-electrochemical model. Therefore, in the simulations presented in this section we have used existing parameter sets and calibrated and validated them against the experimental data.

Still, there are some aspects to discuss about the parameters before we continue. We have used as a starting point the parameters shown in Table~\ref{tab:parameter_values_LGM50}, but some of them needed to be adjusted. The parameters related to the electrochemistry are taken from \cite{Chen2020}. These parameters are assumed to not depend on temperature, except for the reaction rates, which follow an Arrhenius relation. Note that in Table~\ref{tab:parameter_values_LGM50} the diffusion coefficient in the particles is assumed constant, however in real-life it depends on lithium concentration and temperature. Therefore, the diffusion coefficient used in the model is taken as an ``effective coefficient'', and adjusted for each experimental set-up (i.e. each combination of C-rate and temperature), so the constant value mimics the nonlinear dynamics for that particular set-up. That is, for each combination of C-rate and temperature, we change the value of $D_\mathrm{n}$ until we get a reasonable qualitative agreement with experimental data, leaving $D_\mathrm{p}$ unchanged. This assumption is discussed in detail in \cite{Chen2020}. In addition, for each temperature we have adjusted the initial concentration in the positive electrode, so the initial equilibrium potential matches the rest potential before starting the discharge.

The parameters for the thermal model have been taken from \cite{Taheri2013} and then tuned to match the data (in particular they were tuned for 1C at 25 $^\circ \mathrm{C}$) obtaining the values of $\theta = 2.32\E{6}$ $\mathrm{J} \; \mathrm{K}^{-1} \; \mathrm{m}^{-3}$ and $h = 16$ $\mathrm{W} \; \mathrm{m}^{-2} \; \mathrm{K}^{-1}$. Recall that we have neglected the reversible heat generation term $Q_\mathrm{rev}$. However, for the moderate C-rates in the experimental data its contribution should be relatively small. Finally, the value used for both the cell initial temperature and the ambient temperature (which as a modelling assumption we take them to be the same) is taken to be the average value of the final cell temperature of each experiment. This is because the measured temperature of the thermal chamber might not match the equilibrium temperature of the batteries due to experimental limitations in the chamber temperature controller. The parameters that were adjusted for the comparison with experimental data are shown in Table~\ref{tab:tuned_parameters}.

\begin{table
    \begin{subtable}{.35\linewidth}
      \centering
        \begin{tabular}{r | c | c | c}
            & 0.5C & 1C & 2C \\ \hline
            25 $^\circ \mathrm{C}$ & 0.9 & 2 & 6 \\
            10 $^\circ \mathrm{C}$ & 0.4 & 1 & 3 \\
            0 $^\circ \mathrm{C}$ & 0.22 & 0.55 & 1.5
        \end{tabular}
        \caption{$D_\mathrm{n}$ ($10^{-14}$ m$^2$ s$^{-1}$)}
    \end{subtable}%
    \begin{subtable}{.3\linewidth}
      \centering
        \begin{tabular}{r | c }
            & $c_{\mathrm{p},\mathrm{init}}$ \\ \hline
            25 $^\circ \mathrm{C}$ & 17150 \\
            10 $^\circ \mathrm{C}$ & 17750 \\
            0 $^\circ \mathrm{C}$ & 18150
        \end{tabular}
        \caption{$c_{\mathrm{p},\mathrm{init}}$  (mol m$^{-3}$)}
    \end{subtable}%
    \begin{subtable}{.35\linewidth}
      \centering
        \begin{tabular}{r | c | c | c}
            & 0.5C & 1C & 2C \\ \hline
            25 $^\circ \mathrm{C}$ & 24.45 & 24.68 & 24.30 \\
            10 $^\circ \mathrm{C}$ & 9.80 & 10.10 & 9.60 \\
            0 $^\circ \mathrm{C}$ & 0.02 & 0.35 & -0.30 \\
        \end{tabular}
        \caption{$T_\mathrm{amb}$ ($^\circ \mathrm{C}$)}
    \end{subtable}
    \caption{Values for the parameters that have been tuned for the comparison with experimental data.}
    \label{tab:tuned_parameters}    
\end{table}



The comparison between TSPMe, TDFN and experimental data for voltage and temperature is shown in Figures~\ref{fig:comp_exp_25degC}-\ref{fig:comp_exp_0degC}, with the error metrics for TSPMe shown in Table~\ref{tab:error_data}. The first thing we notice is that, in all cases, the TDFN and TSPMe perform very similarly, therefore the errors with experimental data do not arise from the model reduction process. Looking at Figure~\ref{fig:comp_exp_25degC}, we observe a very good agreement of the voltage at 25 $^\circ \mathrm{C}$ at all C-rates. This could be expected, given that the parameters in \cite{Chen2020} were measured at room temperature. On the other hand, we observe that temperature shows very good agreement with data during the relaxation, which shows that the tuned value of the heat exchange coefficient is a good estimation. The discrepancies observed during the discharge phase, especially at C/2, point out that the discrepancy arises from the heat source terms. It is possible that the issue arises from the irreversible heat source term. If that is the case, given that from \eqref{eq:TSPMe_T} we know that the irreversible heating scales like $i_\mathrm{app} \log (i_\mathrm{app})$ while the heating due to Ohmic losses scale like $i_\mathrm{app}^2$, as the C-rate increases the Ohmic heating dominates and the discrepancy due to the irreversible heat becomes less apparent. The discrepancy could also arise from the reversible heat source term which has not been included in the simulations. In both cases, further work is required to validate these hypotheses.

Looking at the voltage comparison at lower temperatures, we notice that they show worse agreement as the temperature decreases. This is reasonable given that, again, the parameters were measured at room temperature. And even the reaction rates, which have a temperature dependence, were measured in the range between 25 $^\circ \mathrm{C}$ and 60 $^\circ \mathrm{C}$ and therefore the values at 10 $^\circ \mathrm{C}$ and 0 $^\circ \mathrm{C}$ are an extrapolation. Moreover, for temperature, we notice a slight increase in the goodness of fit as the ambient temperature decreases, which is not what we expected. This could also be a side-effect of a discrepancy in the heat source terms, as discussed previously.

\begin{figure}[!tb]
    \centering
    \includegraphics[scale=1]{img/comp_exp_25degC.png}
    \caption{Comparison between TSPMe, TDFN and experimental data at 25 $^\circ \mathrm{C}$. The colour solid lines represent the experimental data for the different cells studied (with the cell number matching that of the dataset), the black solid line represents the TSPMe and the grey dotted line represents the TDFN model. The black cross represents the initial equilibrium potential of the model.}
    \label{fig:comp_exp_25degC}
\end{figure}

\begin{figure}[!htb]
    \centering
    \includegraphics[scale=1]{img/comp_exp_10degC.png}
    \caption{Comparison between TSPMe, TDFN and experimental data at 10 $^\circ \mathrm{C}$. The colour solid lines represent the experimental data for the different cells studied (with the cell number matching that of the dataset), the black solid line represents the TSPMe and the grey dotted line represents the TDFN model. The black cross represents the initial equilibrium potential of the model.}
    \label{fig:comp_exp_10degC}
\end{figure}

\begin{figure}[!htb]
    \centering
    \includegraphics[scale=1]{img/comp_exp_0degC.png}
    \caption{Comparison between TSPMe, TDFN and experimental data at 0 $^\circ \mathrm{C}$. The colour solid lines represent the experimental data for the different cells studied (with the cell number matching that of the dataset), the black solid line represents the TSPMe and the grey dotted line represents the TDFN model. The black cross represents the initial equilibrium potential of the model.}
    \label{fig:comp_exp_0degC}
\end{figure}

\begin{table
    \begin{subtable}{.55\linewidth}
      \centering
        \begin{tabular}{r | c | c | c}
            & 0.5C & 1C & 2C \\ \hline
            25 $^\circ \mathrm{C}$ & 72.99 (0.97) & 49.94 (0.96) & 35.99 (0.96) \\
            10 $^\circ \mathrm{C}$ & 116.32 (0.89) & 73.89 (0.87) & 61.83 (0.80) \\
            0 $^\circ \mathrm{C}$ & 99.39 (0.91) & 102.41 (0.65) & 100.81 (0.30) \\
        \end{tabular}
        \caption{Voltage (mV)}
    \end{subtable}%
    \begin{subtable}{.45\linewidth}
      \centering
        \begin{tabular}{r | c | c | c}
            & 0.5C & 1C & 2C \\ \hline
            25 $^\circ \mathrm{C}$ & 0.75 (0.67) & 1.46 (0.83) & 2.05 (0.96) \\
            10 $^\circ \mathrm{C}$ & 0.98 (0.67) & 1.50 (0.88) & 1.56 (0.98) \\
            0 $^\circ \mathrm{C}$ & 1.09 (0.72) & 1.46 (0.92) & 1.11 (0.99) \\
        \end{tabular}
        \caption{Temperature ($^\circ \mathrm{C}$)}
    \end{subtable}
    \caption{Error between TSPMe and experimental data for different temperatures and C-rates. The first values are the root-mean-squared-error (RMSE) and the second values (in brackets) are the coefficient of determination $R^2$. Note that the latter is dimensionless and the closer it is to one, the better the fit is.}
    \label{tab:error_data}    
\end{table}

After comparing the TSPMe and TDFN against experimental data we conclude that the TSPMe can provide an accuracy similar to the TDFN but at a much lower computational cost. The discrepancies observed between the models and the data do not arise from the model reduction, as they are observed in the TDFN as well, so they most likely come from the parameter set. It is well-known that an accurate parameter set is required to obtain reliable predictions, and the parameter set from \cite{Chen2020} does not take into account the temperature dependence of the parameters. Therefore, future work is required to obtain a more accurate parameter set for thermal-electrochemical models so we can validate the model against the LG M50 experimental data.

\section{Conclusions}
In this work, we have derived a reduced thermal-electrochemical model that can be used to simulate the behaviour of lithium-ion batteries in a fast and accurate way (compared to the thermal Doyle-Fuller-Newman model). The reduced model, which we refer to as Thermal Single Particle Model with electrolyte (TSPMe), is a SPM-type model (Single Particle Model) but, as opposed to many of the SPM-type models found in the literature, the TSPMe has been derived systematically using asymptotic techniques. This systematic method allows us to determine in advance the range of validity of the model and the accuracy of the reduced model.

Before presenting the TSPMe, we discussed the concept of the Single Particle Model (SPM), both for pure electrochemical models and coupled thermal-electrochemical models (Section \ref{sec:SPM_discussion}). From that analysis we concluded that the fundamental feature of SPM-type models is that they split the model into two steps: the first step involves solving a system of partial differential equations (PDEs) to determine the intercalated lithium and electrolyte ion concentrations, while the second step involves finding any other variable of interest from closed form expressions. This is a crucial distinction because the computationally expensive part when solving the model is solving the PDEs. Moreover, when comparing the different SPM-type models in the literature we find that the main differences between them are in the second step, thus the system of PDEs is the same across most of the models. The same idea naturally extends when a thermal-electrochemical model is considered. The reduction of the thermal DFN model to the TSPMe was based on only two assumptions: the deviations from the open-circuit potential are small and the heat transfer of the battery is limited by the heat exchange with the environment. With only these assumptions, we can already perform the model reduction, as any other assumption usually taken in other reduced models came as a consequence of the assumptions and the governing equations.

In terms of performance, we found that the TSPMe approximates very well the thermal DFN model, especially at low C-rates, but it is significantly faster. At high C-rates (2C) it is over ten times faster while at low C-rates (C/2) the TSPMe is over forty times faster. Moreover, despite all the limitations of the model (due to the assumptions taken to derive it) and the parameters used in the simulations (which were estimated at room temperature), the TSPMe showed good agreement with experimental data on a commercial cell. This good accuracy and low computational time make the TSPMe model (with other models of the same family) a very good candidate for battery design and control.

There are two main lines to explore as possible extensions of this work. On one side, a full thermal parameterisation, including thermal parameters and thermal dependence of electrochemical parameters in the relevant temperature range (0 $^\circ \mathrm{C}$ to 25 $^\circ \mathrm{C}$), is required for a better comparison between the TSPMe and experimental data. This will allow us to rule out discrepancies due to the parameters rather than the model, and thus obtain a more critical validation of the model. On the other side, the model should also be validated against models resolving the full jellyroll structure of the cylindrical battery (e.g. \cite{Tranter2020}) to assess the discrepancies between the TSPMe, which assumes a lumped thermal model, and the fully resolved model.

\section*{Acknowledgements}
This work is supported by the Faraday Institution [EP/S003053/1 grant number FIRG003] and the HVM Catapult project [number 8164]. FBP and WDW would like to thank Prof Emma Kendrick and Kieran O'Regan, from the University of Birmingham, for the useful discussions.

\section*{Data and code availability}
The code and the experimental datasets are available on the repository:

\url{https://github.com/brosaplanella/TEC-reduced-model} (DOI: doi.org/10.5281/zenodo.4085227).

\section*{CRediT author statement}
\textbf{Ferran Brosa Planella:} Conceptualization, Methodology, Software, Formal analysis, Data Curation, Writing – Original Draft. \textbf{Muhammad Sheikh:} Investigation. \textbf{W. Dhammika Widanage:} Conceptualization, Writing – Review \& Editing, Supervision, Funding acquisition.

\bibliographystyle{unsrt}
