Multi-Scale Thermo-Electrochemical Modeling of Performance and Aging of a LiFePO 4 /Graphite Lithium-Ion Cell

Lithium-ion batteries show a complex thermo-electrochemical performance and aging behavior. This paper presents a modeling and simulation framework that is able to describe both multi-scale heat and mass transport and complex electrochemical reaction mechanisms. The transport model is based on a 1D + 1D + 1D (pseudo-3D or P3D) multi-scale approach for intra-particle lithium diffusion, electrode-pair mass and charge transport, and cell-level heat transport, coupled via boundary conditions and homogenization approaches. The electrochemistry model is based on the use of the open-source chemical kinetics code CAN- TERA, allowing ﬂexible multi-phase electrochemistry to describe both main and side reactions such as SEI formation. A model of gas-phase pressure buildup inside the cell upon aging is added. We parameterize the model to reﬂect the performance and aging behavior of a lithium iron phosphate (LiFePO 4 , LFP)/graphite (LiC 6 ) 26650 battery cell. Performance (0.1–10 C dis-charge/charge at 25, 40 and 60 ◦ C) and calendaric aging experimental data (500 days at 30 ◦ C and 45 ◦ C and different SOC) from literature can be successfully reproduced. The predicted internal cell states (concentrations, potential, temperature, pressure, internal resistances) are shown and discussed. The model is able to capture the nonlinear feedback between performance, aging, and temperature.©TheAuthor(s) 2016. (CC http://creativecommons.org/licenses/by/4.0/), any work [DOI: 10.1149/2.0761702jes] All rights reserved. Manuscript

Mathematical modeling and numerical simulation have become standard techniques in lithium-ion battery research and developmentfrom the atomistic scale up to the system scale. [1][2][3][4][5] Historically, most lithium-ion cell models are based on the work of John Newman and coworkers who developed a one-dimensional model of electrochemistry and mass and charge transport in porous electrodes on the ∼100 μm scale, 6 which was later extended by transport in the active materials particles on the ∼1 μm scale, giving rise to 1D + 1D or "pseudo-2D" (P2D) models. 7,8 This type of model is widely used today. Extensions include solid electrolyte interphase (SEI) formation, 9 aging mechanisms, 10,11 impedance simulations, 12 and multi-phase chemistry in lithium-air 13,14 and lithium-sulfur 15 cells.
Temperature has a strong influence on the performance and lifetime of a lithium-ion battery. A straightforward approach has been to include heat sources and heat conductivity to the standard P2D type models. [16][17][18] However, temperature gradients typically occur on a higher scale, that is, the mm and cm scale of a single cell, as compared to electrode scale described by typical P2D models. Therefore, model extensions to the cell scale are necessary. Consequently, scale-coupling methods have been developed that describe both, electrochemistry on the electrode-pair scale, and heat transport on the cell scale more efficiently. Scale coupling usually uses independent computational domains on the various scales coupled through adequate boundary conditions. As a result, models with various dimensionalities have been presented, for example, 3D + 1D + 1D (cell scale + electrode-pair scale + particle scale), [19][20][21][22][23] 3D (cell scale), 23 2D + 1D (cell and electrode-pair scale + particle scale), 24 2D + 1D (cell scale + electrode-pair scale), 25 3D + 2D + 1D (cell scale + electrode-pair scale + particle scale), 26 and 1D + 1D + 1D (cell scale + electrodepair scale + particle scale). 27 These studies reveal significant temperature gradients inside the cell, emphasizing the requirement for spatially resolved thermal models.
Limited battery lifetime is today a major cost driver, in particular for stationary energy storage applications where long lifetimes are required. 28,29 Aging mechanisms are complex and consist of coupled chemical, electrochemical and mechanical processes. [30][31][32][33][34][35] Aging mechanisms have been integrated in P2D-type models, for example, SEI formation, 9,10,[36][37][38] lithium plating, 11,39,40 capacity loss due to break and repair of the SEI during cycling, 41,42 or active material delamination. 43,44 Due to the nonlinear temperature dependence of ag- * Electrochemical Society Member. z E-mail: christian.kupper@hs-offenburg.de ing mechanisms, temperature gradients inside the cell lead to spatial variations in aging. 45 Thermal models on multiple dimensions are computationally expensive and often require parallel computing algorithms. Aging mechanisms are generally complex and require specialized modeling approaches. Still, both, thermal effects and degradation reactions dominantly influence macroscopic cell behavior. Furthermore, they are strongly and nonlinearly coupled in both directions: aging is highly temperature-dependent, while fast aging reactions cause significant heat release, known as thermal runaway. The coupling occurs over multiple scales-ageing reaction on the nanoscale (particle), performance on the microscale (electrode), and heat transport on a macroscale (cell). A comprehensive understanding of lithium-ion battery performance, safety and lifetime requires therefore both, the prediction of thermal effects, and the prediction of aging effects, under consideration and coupling of the involved scales. Furthermore, the computational expense should be sufficiently moderate to allow systematic dynamic studies over long simulation times for lifetime prediction. In this article we present a multi-scale model of a lithium-ion battery cell that meets these requirements. The model includes the following features: a) 1D + 1D + 1D multi-scale model of particle scale, electrode-pair scale, and cell scale to capture all relevant transport processes under moderate computational effort; b) flexible multi-phase electrochemistry to describe both main and side reactions, including secondary-phase formation, based on the use of the open-source chemical kinetics code CANTERA; 46 c) thermal model capturing the feedback between performance, aging and temperature; d) implementation using a stiff solver that allows both short-term simulation of strong gradients (e.g., during thermal runaway) and long-term simulations (e.g., during calendaric aging). We parameterize and test the model with the example of a cylindrical 26650 lithium iron phosphate (LiFePO 4 , LFP)/graphite cell including SEI formation as aging mechanism.

Modeling and Simulation Approach
1D + 1D + 1D multi-scale transport.-The computational domain of the multi-scale model is shown schematically in Figure 1. It describes transport processes on three distinct scales, each of which is modeled in one dimension, resulting in an overall 1D + 1D + 1D (pseudo-3D or P3D) model. Macroscopic (cell) scale.-On the macroscopic scale (here: x dimension, centimeter scale, cf. Figure 1), we describe heat transport  99 due to conduction using energy balance according to Please refer to the List of Symbols for definition of all symbols and units used throughout this article. We assume that heat conduction takes place along a single through-plane dimension of the cell. In cylindrical cells, this is the radial direction, thus div(λgrad T ) = 1 x ∂ ∂ x (xλ ∂ T ∂ x ). In flat cells (pouch cells), this is the through-cell direction, thus div(λgrad T ) = ∂ ∂ x (λ ∂ T ∂ x ). We therefore explicitly neglect heat conduction in the axial (in-plane) direction. Although axial heat conduction is known to be more effective due to anisotropic heat conductivity, this assumption is at least partially justified by the fact that typical cell aspect ratios are such that the radial (through-plane) distances are much shorter than axial (in-plane) distances. The full resolution of the cell temperature and assessment of the validity of the 1D assumption for a specific cell geometry would require 3D thermal transport models. 47 Boundary conditions at the cell center (x = 0) and at the cell/ambient interface (x = d cell /2, assuming radial symmetry) are respectively. The heat sourceq V in the energy balance Equation 1 results from physicochemical processes on the mesoscopic (electrodepair) scale, as will be discussed in the next section, and requires specific attention when upscaling from electrode-pair to cell scale, as will be described in the Upscaling section.
Mesoscopic (electrode-pair) scale.-On the mesoscopic scale (here: y dimension, hundred micron scale, cf. Figure 1), we describe mass and charge transport in the liquid electrolyte as well as charge transport in the electronic phase in through-plane direction of the electrode pair. In the liquid electrolyte, species conservation is given by and charge conservation follows from assumed local electroneutrality, In these two equations, the three right-hand side terms represent transport fluxes of species i, source term due to (electro-)chemical reactions, and source term due to double layer charging/discharging, respectively. Assuming the presence of an electric double layer at the electrode/electrolyte interface, we describe the electric double layer current as dt , [6] where, the electric-potential difference between liquid electrolyte and solid conducting phase is given as 48 Using Eq. 6, we cast the charge neutrality condition (5) into the form of a differential equation (see Appendix Charge neutrality and double layer capacitance), which we use as governing potential equation. Convection inside the liquid electrolyte is neglected, thus the species fluxes are given by where, the first and second terms on the right-hand side represent diffusion and migration, respectively. Following porous electrode theory, effective transport coefficients are calculated from bulk properties by correcting for porosity and tortuosity, 49 The frequently-used Bruggeman approximation, D eff i = ε elyt 1.5 D i , is a special case of this general expression, based on the assumption that the tortuosity factor τ = τ 2 = ε −0.5 elyt . 50 The general flux Equation 9 can host both, diluted solution theory (DST) (Nernst-Planck type model) and concentrated solution theory (CST), depending on the choice of diffusion and migration coefficients. 13 In the limit of diluted solution, Concentrated solution theory for a binary electrolyte can be cast into the form of Eq. 9 by selecting (cf. Appendix Transport coefficients from concentrated solution theory) where, the index i refers to the two ions (e.g., Li + and PF 6 -). The diffusion coefficient D 0 , electrolyte conductivity σ, transference numbers t 0 − = 1 − t 0 + , and lumped activity parameter ν = (1 − t 0 + )(1 + ∂ ln f∓ ∂ ln c ) are the common parameters applied in lithium-ion battery models with concentrated solution theory. 51 Typical lithium-ion battery electrolytes are concentrated solutions, therefore CST provides a more accurate picture of species transport. However, the parameterization of CST requires specialized and timeconsuming experiments. 51 In many cases, parameters may be unavailable, for example because the electrolyte composition of a commercial cell is unknown. Here, DST can be advantageous because the only parameters are the diffusing species' diffusion coefficients which can be relatively easily derived from cell impedance experiments. Using a correctly-parameterized DST model may give more accurate simulation results than a CST model parameterized to a different type of electrolyte. Also, the DST framework in Eqs. 9-11 can host electrolytes containing more than two ionic species, for example, polysulfides in lithium/sulfur cells. 15 Here, CST is even more difficult to A306 Journal of The Electrochemical Society, 164 (2) A304-A320 (2017) parameterize. Therefore, it is advantageous to have a flexible model framework that can host both, CST and DST.
We furthermore assume that the electronic resistance within the electrode is negligible compared to the ionic resistance in the liquid electrolyte, resulting in a spatially constant φ elde throughout each electrode. This assumption is valid for typical lithium-ion battery electrode materials. Graphite electrode conductivities are two orders of magnitude higher than that of the electrolyte. [52][53][54] Although pure lithium iron phosphate is a very bad electron conductor, 55-57 thanks to carbon coating or conductive additives used in commercial LFP batteries the electrode is likely to have one order of magnitude higher conductivity than the electrolyte. 51,[58][59][60][61][62] Boundary conditions for the species conservation are J i = 0 at both electrode/current collector interfaces. For the ionic potential, the boundary conditions are d( φ)/dx = 0 at both electrode/current collector interfaces. Electrode-pair voltage follows from where, we assume an additional potential drop due to the electronic resistance R cc of the current collection system. We set φ elde,ca = 0V as potential reference. The externally applied current density i is given as Heat production takes place within the electrodes. We assume that, on the electrode-pair scale, temperature gradients are small with respect to gradients on the cell scale. Therefore, the electrode pair itself is modeled isothermally, yet with time-dependent temperature as given by the cell model (cf. Upscaling section). The area-specific heat is given by integrating all local heat sources over the electrode-pair length,q A = L EP 0 (q chem (y) +q ohm (y)) dy + R cc i 2 , [16] where, reversible and irreversible heating due to all chemical and electrochemical reactions is given by (see Appendix Heat source terms) q chem = Nr n=1 r n A V n − H n + Fν e,n φ n , [17] and ohmic (Joule) heating in the liquid electrolyte is given bẏ [18] The last term in Eq. 16 represents ohmic heating due to resistance of the current collection system. In summary, we have presented governing equations for species concentrations (Eq. 4) and electric potential (Eq. 8) on the electrodepair scale, as well as the required constitutive equations. We still require a closing relationship between faradic current and electricpotential difference, i V F = f ( φ eff ). This relationship follows from electrochemistry, as described in the Electrochemistry and multiphase chemistry section.
Microscopic (particle) scale.-On the microscopic scale (here: z dimension, micrometer scale, cf. Figure 1), we describe diffusive mass transport inside the active materials particles. We assume spherical particles and Fickian diffusion. The mass conservation for lithium inside the active material (Li,AM) is given by The solid-state diffusion coefficient may depend on the concentration of intercalated lithium c Li,AM . Boundary conditions are J Li,AM = 0 in particle center and J Li,AM = r P 3ε AMṡ V Li,AM at particle surface. The geometric factor r P 3ε AM is derived in Appendix Geometric factor. In some cases we do not want to model solid-state diffusion on the microscale, for example, if current rates are low (and we want to save computational time); for a model-based sensitivity analysis (where simulations runs with and without solid-state diffusion are compared); or if the active material is showing a phase-change behavior (such as lithium iron phosphate 63 ) and therefore requires a different treatment of the microscale. Mass conservation of bulk lithium is then given by ∂ε AM c Li,AM ∂t =ṡ V Li,AM , [20] which can be used instead of Eq. 19.
Upscaling.-Each of the three scales described above represents transport processes modeled in one dimension. This section describes the coupling strategy between the scales.
Between macroscale (cell) and mesoscale (electrode pair) we use a homogenization approach, as shown in Figure 2. We use a reduced number of electrode-pair models N EP (typically, 1. . . 10) which we distribute along the x direction. Each electrode-pair model m thus represents the behavior of a hollow cylinder V m of the cell. The cell scale passes the local (on the x-scale) temperature T m down to the electrode-pair scale. In turn, the representative electrode-pair models pass their heat sourcesq V m up to the cell scale, where they are used locally (on the x scale) within the radial sections represented by the electrode-pair models.
Lithium-ion battery cells consist of stacked or wound sheets of electrode pairs with a total area in the m 2 range (e.g., 0.171 m 2 for the 26650 cell investigated here). This area is referred to as "active electrode area" A e . 11 Using this parameter, the volumetric (on the x scale) heatq V m is calculated from the area-specific heatq A m (Eq. 16) of electrode pair m according tȯ The total cell current is given as weighted sum over the representative electrode pairs, while the cell voltage is identical for all electrode pairs due to their parallel connection (assuming negligible potential differences throughout the current collection system, due to relatively small resistance of the metal foils of about 10 μ compared to the full cell resistance in the order of > 1000 μ 64 ), [23] Number and size of representative hollow cylinders are typically selected as to resolve the highest temperature gradients. We have empirically observed that even in the limiting case of one single volume section (i.e., one electrode pair represents the complete cell), the thermal cell behavior can be well described. This is due to the fact that, the spatial discretization of the heat transport Equation 1 (indicated as solid line in Figure 2) can be chosen independently from the number of electrode-pair models (indicated as solid points in Figure 2). The latter example would lead to a temperature gradient under a constant heat source along the x-scale.
Between mesoscale (electrode pair) and microscale (particle) we apply the "standard" pseudo-2D (P2D) approach, 7,8 that is, a single particle is modeled at every grid point of the electrodes.
Electrochemistry and multi-phase chemistry.-We use a generalized multi-phase chemistry framework that was described in detail in Ref. 13 and that is based on the use of the open-source chemical kinetics code CANTERA 46 (cf. Simulation methodology section). Briefly, the electrode pair ( Figure 1, y scale) is assumed to consist of up to seven layers (positive and negative electrode, positive and negative current collectors, separator, positive and negative gas reservoirs). Each layer may host an arbitrary number of bulk phases (solid, liquid, or gaseous) characterized by their respective volume fractions ε, where each bulk phase may host an arbitrary number of chemical species. Each layer may furthermore host an arbitrary number of interfaces, characterized by their volume-specific surface area A V , where each interface may host an arbitrary number of interfacial reactions between adjacent bulk phases.
The following Eqs. 24-26 represent the implementation of chemical kinetics in CANTERA. The rate of a single interfacial reaction follows from mass-action kinetics, 65,66 [24] where the concentrations c i refer to the concentrations at the electrode/electrolyte interface as given by the transport models. The forward rate constant is given by transition state theory through an extended Arrhenius expression, 65,67 which consists of a thermally-activated part and a potential-dependent part. The reverse rate constant follows from thermodynamic consistency according to where, c 0 i is the standard concentration (reference state 66 ) of species i. CANTERA uses the following definition of standard concentrations: for bulk liquid and solid phases (ConstDensityThermo class), c 0 i = ρ/M (total concentration in the phase); for the gas phase (IdealGas class), c 0 i = pgas RT (total concentration in the phase). The last term in Eq. 26 is required for unit consistency in Eq. 24. Note the interfacial reaction rate r has units of mol m 2 s , while the forward rate preexponential factor k 0 f has units of mol In an ideal intercalation material, the equilibrium half-cell potential depends on the (concentration-independent) standard Gibbs energy G 0 and the concentrations of the involved species according to the Nernst equation, [27] which can be derived from Eqs. 24-26 by setting r = 0 and φ = φ eq . A real intercalation material shows additional concentration dependencies that can be described using a (concentration-dependent) excess Gibbs energy G E , 18,65,68 In our simulation approach, we take this into account by using a concentration-dependent Gibbs reaction energy G ( 26. In CANTERA, the Gibbs reaction energy is calculated from the species molar enthalpies and entropies, [29] We parameterize h Li (c i ) and s Li (c i ) of intercalated lithium species from experimental half-cell potentials assuming reference values for the vacancies h V = 0, s V = 0. The effective potential difference between electrode and electrolyte φ eff occurring in Eq. 25 follows from its nominal value (Eq. 7) by correcting for a film resistance, representing, for example, the SEI (cf. Calendaric aging section). The overall (volumetric) bulk species source terms entering the mass balance Eq. 4 are obtained by summing over all interfacial reactions, scaled to the specific area of the respective interface, The faradic current entering the charge balance Eq. 8 is obtained in analogy from the source terms of electrons, F ν e,n r n A V n [32] The interfacial reactions cause a net mass transfer between bulk phases, leading to a change of volume fractions. Therefore, we introduce additional conservation equations for all bulk phases, 13 Note that the present electrochemistry framework does not use the Butler-Volmer equation. Instead, the more fundamental transition state theory, Eqs. 24-26, is used. In order to cast this theory into Butler-Volmer type kinetics, we additionally multiply the forward rate constant k f by a factor exp(−α f G 0 RT ) (see Appendix Mass action kinetics and Butler-Volmer formulation).
In summary, this section has presented complete expressions describing electrochemical thermodynamics and kinetics. It can generally host an arbitrary number of charge-transfer and non-chargetransfer reactions. The required model parameters are the forward rate preexponential factor k 0 f , activation energy E act,f , and symmetry factor α f (in case of charge-transfer reactions) of all reactions, as well as molar enthalpies h i and molar entropies s i of all species. The potential-difference φ enters from the meso-scale transport model.
Gas pressure inside the cell.-SEI formation during aging causes gas release, hence, pressure increase inside the cell housing. While a pouch cell can (partially) accommodate pressure increase by volume expansion, in the cylindrical cell studied here the housing is rigid, thus only pressure rise needs to be addressed. In order to quantify this effect, the void spaces of the cell housing are modeled as 0D gas reservoir. 69 The governing equations are the continuity equation, A308 Journal of The Electrochemical Society, 164 (2) A304-A320 (2017) species conservation, [35] and ideal gas law, The gas-phase species source termṠ gas i follows from integration over the complete electrode volume on mesoscale and macroscale, similar to the electrical current (Eqs. 15 and 22), according tȯ The gas pressure and composition is assumed homogenous throughout the cell, including the gas-filled porosity inside the anode.
Simulation methodology.-The governing equations and relationships presented above (except the chemistry, see below) are implemented in the in-house multiphysics software package DENIS (Detailed Electrochemistry and Numerical Impedance Simulation). 13,48 DENIS is a C/C ++ code with modular structure. Its main functionality is to cast the model equations into the form of a differential-algebraic equation system, based on conveniently editable text input files. The solution vector y in the particular case presented in the results has a dimension of 1928. We apply the finite-volume method to discretize the partial differential equations for T (Eq. 1), c i (Eq. 4), φ (Eq. 8) and c Li,AM (Eq. 19) using non-equidistant non-adaptive discretization. The number of grid points is 10, 31 and 14 in x, y, and z dimensions, respectively. We use N EP = 5 electrode-pair models distributed unevenly on the x scale, representing radial sections of 0.5, 4, 4, 4 and 0.5 mm of the cell.
The DAE system (38) is numerically solved using the implicit timeadaptive solver LIMEX (version 4.3A). 70,71 The equation system can be either solved for I cell when E cell is given as independent variable (potentiostatic simulation), or vice versa (galvanostatic simulation). The complete electrochemistry (thermodynamics and kinetics) is evaluated using the chemical kinetics software CANTERA (version 2.2). 46 In particular, CANTERA implements Eqs. 24-26 for reaction mechanisms of arbitrary complexity. The reaction mechanisms as well as all thermodynamic and kinetic data are provided to CAN-TERA in the form of conveniently editable text input files, which are available from the authors upon request. CANTERA is a C ++ code which we interface from DENIS to obtain the reaction rates r needed in Eqs. 17, 31, 32, 37, and reaction enthalpies H needed in Eq. 17. For all bulk phases we use CANTERA's ConstDensi-tyThermo ("incompressible_solid") model, except for the gas phase for which we use the IdealGas model. For all interfacial reactions the InterfaceKinetics model, as described by Eqs. 24-26, is used. The stoichiometry-dependent thermodynamics of intercalated lithium species is implemented via an additional DENIS/CANTERA feedback loop.
We use MATLAB (version 2016a) for controlling all DENIS simulations via S-function interfaces, as well as for data evaluation and visualization. Wall-clock computational time is typically 4 hours on a 3.4 GHz Pentium processor desktop computer (single-core simulation) for one discharge curve as shown, for example, in Figure 5. This comparatively long computational time is due to the fact the code was designed for flexibility 13,48 (including full coupling to CANTERA), not numerical efficiency.

Results and Discussion
Parameterization.-While the modeling approach presented above is generic in the sense that it can host different types, geometries and chemistries of lithium-ion cells, we here specifically parameterize the model in order to represent a commercially-available cylindrical 26650-size high-power LFP/graphite cell. The large number of required parameters poses a particular challenge to physically-based battery modeling. 72 We use the following parameterization approach: (1) Use literature values on the identical or similar cells and components for an initial set of model parameters; (2) use macroscopic experimental data, in particular electrical and thermal behavior upon charge and discharge, to re-parameterize selected performance-sensitive parameters, in particular, rate coefficients (preexponential factors the charge-transfer reactions) and thermal parameters (heat transfer coefficient); (3) use macroscopic aging data from literature (capacity as function of calendaric aging time) to parameterize the rate coefficients of the aging reaction (pre-exponential factor, activation energy); (4) test the model against macroscopic experimental data over a wide range of conditions (charge/discharge rates from 0.1 C to 10 C, temperatures, SOC). Parameter identification was performed by manually varying values (automated fitting was not possible due to high computational time of the simulation). The resulting full set of model parameters is given in Table I (macroscopic  cell scale), Table II (mesoscopic electrode-pair scale), and Table III (microscopic particle scale). In addition, Figure 3 and Figure 4 show the thermodynamics and the transport data of intercalated lithium as function of intercalation stoichiometry. A particular strength of the present modeling approach is the ability to host complex multi-step electrochemical reaction mechanisms. The thermodynamic data of all used species are given in Table V, and the reaction equation and kinetic coefficients are given in Table VI. LFP is known to be a two-phase material with charge/discharge equilibrium-potential hysteresis. 73,74 In the present model we neglect the hysteresis and assume that charge and discharge thermodynamics are described by the same data ( Figure 3). More sophisticated models that are able to capture the equilibrium-potential hysteresis have been published 75 and may be integrated into the present framework in future, but are out of scope of the present study. As LFP shows a two-phase behavior during charge and discharge, the model of spherical diffusion used for the graphite anode is not applicable here 56,76,77 and requires a more complex mathematical description like the coreshell model used in Dargaville et al. 78 or the multi-particle model of Farkhondeh et al. 75 Again, these may be integrated into the present framework in future, but are out of scope of the present study. Therefore, bulk diffusion is not modeled at the cathode, recognizing that the model starts to become less reliable for high C-rates (>5 C), as will also be shown below.
The present framework is able to host both, DST and CST transport models of the electrolyte (cf. Mesoscopic (electrode-pair) scale section). While we have investigated both models, the simulation results show rather small differences (results not shown). As the CST parameters are not known for the electrolyte of the investigated cells, the results shown in the remainder of this section are based on DST.
It should be noted that a number of other parameters have to be considered as "preliminary", either because the true materials and  composition of the commercial cells used for this study are not known, or because there is no experimental data available for these materials. In this case, available data from similar materials and composition were used. By fitting selected model parameters to experimental data on the investigated cells, deviations of "preliminary" parameters from their true values and deviations due to simplified modeling assumptions (two-phase behavior, DST) are implicitly included in the fits. Note that fitted parameters may still not be unique, in particular when parameters show similar sensitivities to cell behavior (e.g., preexponential factors of cathode and anode charge-transfer reaction). Still, our parameterization approach is justified by the fact that the thermal and electrical cell performance and aging behavior can be successfully reproduced over a wide range of operating conditions, as will be shown in the following.   3.6 V with C/10 final current. The model is able to quantitatively reproduce the experimental electrical and thermal behavior over the complete range of investigated conditions. Going into detail, there are a number of aspects where simulation and experiment deviate, as discussed in the following. (a) The simulated cell capacity shows lesser dependence on C-rate than the experiments. We believe that this is an artefact of the experiment, as the data at 0.1 C and 1 C were carried out   with a different individual cell than the data at 5 C and 10C; note the model parameters were fitted to the data at 0.1 C and 1 C. (b) At 5 C and 10 C, simulations show an increasingly different voltage than the experiments. Overpotential is underpredicted in case of discharge and overpredicted in case of charge. This is likely a result from neglecting micro-scale (particle-scale) transport in the cathode and phase-change effects in the present simulations (cf. Parameterization section). (c) As consequence of this, the temperature increase is underpredicted in case of discharge and overpredicted in case of charge at high C-rates. (d) The details of the shape of the discharge and charge curves are not fully reproduced by the model. Likely, the applied thermodynamic data and electrode balances do not fully correspond to those of the experimental cells.
Reducing the remaining disagreement between simulation and experiment would certainly be feasible by revising the parametric base. However, this would not contribute to the mechanistic understanding of the cell behavior that will be discussed below, and is therefore out of scope of the present work. Note also that empirical equivalent circuit-based models typically show a much better agreement with experiments than physically-based models, which however cannot contribute to understanding of the internal chemistry and transport mechanisms.
Internal cell states.-In this section we discuss the internal cell states during constant-current discharge at 5 C and 25 • C ambient temperature of a fresh (non-aged) cell. This condition was chosen as representative example for the multi-scale modeling and simulation approach. All results therefore correspond to one single data set shown in Figure 5 and Figure 6.
On the macroscopic (cell) scale, temperature profiles are shown in Figure 7 for 16 equidistant time steps between beginning (t = 0 s) and end (t = 720 s) of discharge. Temperature shows a gradient between center (x = 0) and surface (x = 13 mm) of the cell. While the surface temperature increases from 25 • C to 40.3 • C during discharge, the core temperature increases to 43.5 • C. The temperature difference between surface and core remains relatively low but still significant. Erhard et al. 26 have simulated and measured temperature differences within an LFP/graphite cell. Their highest C-rate (2 C) showed a maximum core temperature of about 40 • C and a maximum temperature difference to the surface of about 5 • C. They have emphasized the strong influence of this relatively small temperature difference on local cell properties like voltage or current density.
The spatiotemporal behavior of the internal states on the mesoscopic (electrode pair) scale is shown in Figure 8. These profiles were taken at the center of the cell (macroscale x = 0). The concentration of Li + and PF 6 ions in the electrolyte is shown in Figure 8a. Note the concentrations of Li + and PF 6 are identical, as enforced by charge neutrality (cf. Mesoscopic (electrode-pair) scale section). They both show a gradient between negative electrode (high concentration) and positive electrode (low concentration), consistent with the location of formation and consumption of Li + , respectively. The gradient becomes smaller for progressing discharge. This is due to the increas-  ing cell temperature (cf. Figure 7), increasing the ion diffusivity of the electrolyte. Figure 8b shows the normalized electric-potential distribution in the electrolyte. Its absolute value (not shown) changes by around 1 V during discharge due to the changing electrode potential. For Figure 8b, we have normalized the potential to the value at separator center to make the small potential gradient (ca. 10 mV) visible. The potential drop is decreasing with progressing discharge because of the increasing cell temperature, increasing the ionic conductivity of the electrolyte. The potential is increasing from the positive to the negative electrode. This drives a migration flux of Li + from the negative to the positive electrode, which adds to the diffusive flux due to the concentration gradient (Figure 8a) in the same direction. For PF 6 -, migration flux is in the opposite direction, i.e., from positive to negative electrode. It is cancelled by diffusive flux, resulting in a net zero flux for the PF 6 ion. The combined fluxes satisfy both, charge neutrality and net Li + transport. Figure 8c shows the lithium stoichiometry in the active materials. At the positive electrode, it increases (lithium intercalation), while at the negative electrode it decreases (lithium de-intercalation) during discharge. Importantly, the data show a considerable spatial gradient of the stoichiometry-in other words, a spatial distribution of the local SOC. This is more pronounced at the positive electrode compared to the negative electrode.
Finally, on the microscopic (particle) scale, the distribution of lithium stoichiometry in the graphite anode particle is shown in Figure  9. These data were taken at the center of the cell (x = 0 mm) for a representative particle close to the anode/separator interface (y = 100 μm). The stoichiometry continuously decreases during discharge, starting at the particle surface (z = 0 μm). The wave-like shape of the profiles is a result from the non-constant solid-state diffusion coefficient (cf. Figure 4). The diffusion coefficient has a minimum at a stoichiometry of ca. 0.3, corresponding to the strongest gradient in the stoichiometry profiles.
Calendaric aging.-The results shown in the previous section were obtained with a fresh (non-aged) cell. In the present section, we use the model to investigate calendaric aging. We assume that aging is due to SEI formation at the anode/electrolyte interface according to   [39] This reaction takes place in parallel to the main charge-transfer reaction (interalation/deintercalation). Note that SEI formation was also included in all simulations shown in the previous sections, however did not affect the results shown there due the short time scales of simulation.
The film resistance entering Eq. 30 is given by (see Appendix SEI film resistance and formation rate) SEI film growth is known to be diffusion-limited, resulting in a growth rate proportional to 1/δ SEI (see Appendix SEI film resistance and formation rate). We include this into the model by scaling the preexponential factor of the SEI formation reaction with a factor of 1/δ SEI (cf. Table VI). The details of the complex multi-step SEI formation mechanism 36 are not resolved here, but are assumed to be implicitly included in the fit of the rate coefficients to experimental ageing data. Additionally the typical dry-out of the cell during aging due to electrolyte consumption is considered by multiplying the specific area of the anode active material surface with the aging factor ε elyt (t) ε elyt,init (cf .  Table VI).
In order to simulate calendaric aging, all parameters were initialized and the electrode stoichiometries were set to 100% SOC, representing a fresh fully-charged cell. The cell was virtually aged  by carrying out a transient simulation with zero-current boundary condition for different times up to 500 days. After aging, the capacity of the virtual cell was determined by first fully charging (1 C, CCCV, 3.6 V cutoff, C/10 break, 1 h resting time) and then fully discharging (1 C, CC, 2.0 V). The same types of simulations were carried out for different initial SOC and ambient temperature. Note the wall-clock time for 500 days aging simulation is only 10 min. Figure 10 shows simulated (this work) and experimental (Grolleau et al. 84 ) state of health (SOH) as function of aging time for three different states of charge at 30 • C ambient temperature. Here, the SOH is defined as . [41] Within the scatter of the experimental data, the simulation is able to reproduce capacity loss both qualitatively and quantitatively for 30 • C. Capacity loss is up to 5% for 500 days calendaric aging under the investigated conditions. At 45 • C, the experiments show a stronger linearity than the simulation, which follow a √ t behavior as expected from diffusion-limited film growth. 85,86 The data indicate that additional degradation mechanisms not considered in the model become relevant at elevated temperature, for example, cathode decomposition due to precipitation of iron particles. 87 It should be noted that only two fit parameters of the model are associated with aging, that is, pre-exponential factor and activation energy of the SEI formation reaction (cf. Table VI). The SOC dependence results from the simulated internal potential distribution (cf. Figure 8b) driving the chargetransfer SEI formation reaction. Note that the fitted activation energy (106 kJ/mol) is higher than that reported before (43 kJ/mol). 88,89 The origin for this difference may be associated with the different types of investigated cells and different chemical ageing models.
In order to further investigate aging behavior, we have simulated the internal resistance of a fresh cell and an aged cell (500 days, 65% SOC, 30 • C). The internal resistance of the virtual cell was determined following a typical experimental protocol: [90][91][92] Perform complete charge/discharge cycle to determine the reference capacity for SOC and SOH definition; full charge using a CCCV protocol; discharge with 1 C rate for 15 min, corresponding 25% of the SOC; rest for 15 min; discharge at 0.1 C rate for 30 s; current step to 1 C and discharge for 30 s; rest for 15 min. The internal resistance was determined from the current and voltage values before and 3 s after the step according to [42] In addition, we have used the simulated internal cell states to quantify the contributions of the different cell components to the internal resistance, in particular, anode, cathode, separator, and current collector. Results of this analysis are shown in Figure 11, where a fresh cell and an aged cell are compared. The internal resistance of the fresh cell shows a distinct SOC dependence well-known from experiments, 30,92 with a minimum at intermediate SOC and maxima for both full and empty cells. This shape is a result from anode and cathode, while separator and current collector show SOC-independent behavior. The aged cell shows an increased internal resistance due to SEI film resistance (cf. Eq. 30 and Eq. 40). Interestingly, internal resistance increase is pronounced at 0% SOC, while at the same time internal resistance of the cathode decreases. This is caused by electrode imbalancing during aging. 45 In our simulations the stoichiometry range at the cathode (x in Li x PO 4 ) changes from 0.012-0.990 (fresh cell) to 0.012-0.915 (aged cell) and that of the anode (x in Li x C 6 ) from 0.010-0.590 (fresh cell) to 0.0064-0.562 (aged cell). Activation overpotential and therefore electrode internal resistance strongly increase toward high and low intercalation stoichiometries. 93 An SOC of 0% corresponds to the upper x value in Li x PO 4 , which is decreasing upon ageing, and to the lower x value in Li x C 6 , which is also decreasing upon ageing, therefore causing the observed variation in electrode resistances.
Apart from the macroscopic electrical behavior in terms of capacity and internal resistance, the physically-based model allows insight into internal cell states that are difficult to access experimentally. Figure  12 shows the volume fractions of the phases of the negative electrode as function of time for a 500-day calendaric aging simulation. SEI volume fraction continuously increases from its initial value of 1.0% to 4.68%. The latter value corresponds to a film thickness of 75.6 nm. At the same time, electrolyte volume fraction decreases by 5.11% due to its decomposition (cf. Reaction 39). The reaction product SEI has an overall higher density than the reactant electrolyte (cf. Table  IV). The resulting "void" volume increases by 1.69%, which can be directly interpreted as drying-out of the cell. The "void" volume in the electrode is actually filled with the formed gaseous ethylene. The gas also fills the "void" volume of the cell beyond the wound electrode pairs (cf. Table IV). As a result of the formed gas, the pressure inside the cell increases. Pressure as function of time is shown on the right axis of Figure 12. The simulation predicts a pressure increase to over 6 bar after 500 days aging time. This value is in the same order of magnitude compared to values from literature where cells with graphite anode and similar electrolyte reached 3 to 5.5 bars at 60 • C and 100% SOC after just 50 days of aging. 31 The high value might be due to one or several shortcomings in the model: (a) the initial void volume fraction inside the cell and/or the electrodes might have been underestimated; (b) the mechanical expansion of the cell housing and/or compression of the cell components due to increasing gasphase pressure is not considered in the model; (c) potential follow-up reactions of the gaseous ethylene are not considered.
The multi-scale model is able to reproduce all expected physicochemical changes inside the cell at least qualitatively, including   capacity loss, internal resistance increase, electrolyte decomposition, SEI formation, and pressure buildup. Quantitative conclusions are subject to the uncertainty of the model parameters (cf. Parameterization section). Further validation, in particular of the internal cell behavior, is beyond the scope of the present study, but will be subject of future investigations.

Conclusions
Lithium-ion batteries show a complex thermo-electrochemical performance and aging behavior. We have presented a modeling and simulation framework that is able to describe both multi-scale heat and mass transport and complex electrochemical reaction mechanisms.
The transport model is based on a 1D + 1D + 1D (pseudo-3D or P3D) multi-scale approach. Heat transport in the radial cell direction (1D, macroscale) is modeled as conductive process. Mass and charge transport on the electrode-pair scale (1D, mesoscale) is modeled as diffusion and migration. The charge neutrality condition is cast into a time-dependent partial differential equation by including double layer charging/discharge, allowing stable numerical simulation. Intraparticle transport of lithium atoms (1D, microscale) is modeled as Fickian diffusion with concentration-dependent diffusion coefficient. A 0D model of the void cell volume was added, allowing to describe gas-phase species concentration and pressure buildup during aging.
The electrochemistry model is based on the use of the open-source chemical kinetics code CANTERA, which is coupled to the transport model via the chemistry source terms. In this approach, charge-transfer reaction kinetics are not modeled with Butler-Volmer equations, but with more fundamental relationships based on mass-action kinetics and transition state theory. The framework allows to describe multireaction multi-phase thermoelectrochemistry. In the present study we couple the main reactions (intercalation at the particle surface) to SEI formation as parallel side reaction at the anode.
We have parameterized the model to reflect the performance and aging behavior of an LFP/graphite 26650 battery cell. The model was demonstrated against cell voltage and surface temperature (0.1-10 C discharge/charge at 25, 40 and 60 • C) as well as calendaric aging experimental data (500 days at 30 • C and 45 • C and different SOC) from literature. The predicted internal cell states (concentrations, potential, temperature, pressure, internal resistance) during 5C discharge and during calendaric aging were shown and discussed. The model is able to capture the nonlinear feedback between performance, aging, and temperature.
Future work will address the further parameterization and validation of the model using an extended experimental data base; and the application of the model to various scenarios, including sensitivity analyses, thermal runaway, and lifetime prediction.

Acknowledgments
CK acknowledges financial support from the state of Baden-Württemberg through the cooperative graduate college DENE (Dezentrale Erneuerbare Energiesysteme). This work has been partially funded by the German Ministry of Education and Research (BMBF, grant no. 03FH013PX3). The authors thank Michael Danzer and Harry Döring (Zentrum für Sonnenenergie und Wasserstoffforschung, ZSW, Ulm) for providing experimental data of lithium-ion battery cells as part of a project founded by the Volkswagen Foundation. WB thanks Timo Danner, Michael Hedwig, Birger Horstmann, Christian Hellwig, Manik Mayur, Wolfgang Mielke, Sundar Sathyamoorthy, Svenja Spitznagel, and Nanako Tanaka for their various contributions to model development, implementation, and parameterization.

Appendix
Charge neutrality and double layer capacitance.-In the following we describe the reformulation of the charge-neutrality Equation 5 into a differential equation for the local potential difference, Eq. 8. We assume that the electric current due to double layer charging and discharging is described as (Eq. 6) dt . [A1] Mass and charge conservation within the double layer requires that the electronic double layer current entering the electronic phase (solid conductive matrix) is counterbalanced by specific adsorption/desorption of ions from/into the ionic phase (liquid electrolyte). This results in an additional species source term in the continuity Equation 4, In the present study, we assign the double layer current to one single ionic species, that is, i = Li + . The faradic current i V F follows from the production rate of electrons in the electronic phase due to charge-transfer reactionsṡ V e , [A3] Due to overall charge neutrality of all (electro) chemical reactions, Inserting Eqs. A1-A4 into the charge neutrality Equation 5 yields a new governing equation for the electric potential difference φ, This is the governing potential equation that we solve for. It has the advantage of being a differential equation instead of an algebraic Equation 5, allowing a more straightforward numerical implementation and numerically more stable simulations.
Transport coefficients from concentrated solution theory.-The electrolyte transport model described in the Mesoscopic (electrode-pair) scale section can host CST as shown in the following. We compare here to the CST model of a binary electrolyte as originally developed by Newman and co-workers, where the solvent is used as reference species and its velocity is taken to be zero (i.e., convection is neglected). 8,94 Transference number is considered with respect to solvent velocity. For a binary electrolyte described by concentrated solution theory, Eq. 9 reads Using expressions (12)- (13) and assuming effective transport parameters, the transport coefficients are expressed by 13 where, the index i refers to the two ions (e.g., Li + and PF 6 -) denoted by subscripts + and − in the following, and Journal of The Electrochemical Society, 164 (2) A304-A320 (2017) Inserting A7-A9 into Eq. A6 yields [A13] We first derive the governing equation for the ionic potential. The ionic current density i l (in A/m 2 ) in the electrolyte is given by 66 [A14] In our case, Inserting Eq. A13 and using Eqs. A10-A12 yields [A16] Finally inserting into Eq. A9 gives the potential equation, This is the standard expression used in literature for a concentrated binary electrolyte solution. 11 We next derive the governing equation for the ionic concentration. We start from Eq. 4 for the cation (in the following, c = c Li + and J = J Li + ), Inserting Eq. A13, using Eq. A11 and A16, and assuming concentration-independent t + gives This is the standard literature expression for a concentrated binary electrolyte solution. 95,96 If we recognize thatṡ Li + = 1 z + F ∂i l ∂t , the right-hand side modifies to, which is the standard literature expression for a concentrated binary electrolyte solution with source term in a porous electrode. 11 Note that, in the absence of a source term (e.g., in the separator), ∂i l ∂ y = 0 and so the second term on the right vanishes. We have therefore shown that our formalism introduced in Eqs. 4,9,12,13 is identical to the standard literature formalism, Eqs. A17 and A19.
Heat source terms.-We apply here a spatially-resolved (on the mesoscopic scale) description of heat sources. The total energy released by an interfacial chemical reaction is given by its reaction enthalpy, and, in case of a charge-transfer reaction, is distributed among thermal and electrical energy release. The released electrical energy is given by the work of the electron during its movement over the electrode/electrolyte potential step according tȯ Therefore,q We can further reformulate this expression by recognizing that i V F = Fṡ V e = Fν e r A V , yieldingq [A24] Generally, an arbitrary number of chemical and electrochemical reactions can take place at different interfaces within the electrode. The released heat follows as sum over all reactions,q heat = Nr n=1 r n A V n − H n + Fν e,n φ , which is the expression used in our model (Eq. 17). It is valid for both electrochemical and thermochemical reactions, as for the latter, ν e = 0. Eq. A24 implicitly includes both reversible and irreversible heat contributions. We will next show that it can explicitly be cast into a function of these contributions. The activation overpotential is defined as, where, the equilibrium half-cell potential difference follows from thermodynamics according to We convert Eq. A24 to a function of current according tȯ Inserting Eqs. A26 and A27 and recognizing ν e = −z for a charge-transfer reaction written in standard reduction direction yieldṡ With the temperature derivative of Eq. A27 this can be further converted tȯ In this expression, the first and second terms on the right-hand side represent reversible and irreversible heat, respectively. Eq. A30 is often used in thermal models of lithium-ion batteries. [16][17][18]20,23 In our physically-based modeling approach, we use instead the implicit form of Eq. A24.
Geometric factor.-Volume-specific interfacial area and volume fraction of an active material particle in a composite electrode are defined as where, V is the electrode volume. Assuming a spherical particle, volume and surface area are given by A AM,P = 4πr 2 P , respectively. Inserting Eqs. A32-A34 into Eq. A31 yields an expression for the specific area as function of volume fraction, Upon intercalation, the volumetric source termṡ V Li,AM is converted to a surface flux via Inserting Eq. A35 into Eq. A36 yields the geometric factor introduced in the Microscopic (particle) scale section.

Mass action kinetics and Butler-Volmer formulation.-The
Butler-Volmer equation can be derived from transition state theory 65,66 by inserting Eqs. 25-26 into Eq. 24, substituting the absolute potential difference φ eff with the activation overpotential, and converting the reaction rate to current via i = z Fr. [A38] Extensive algebraic manipulation results in the Butler-Volmer form with the exchange current density and the exchange current density factor Eq. A41 shows that i 00 is a function of G, which, in lithium-ion batteries, is a function of SOC (cf. Electrochemistry and multi-phase chemistry section). In order to ensure SOC independent i 00 , we multiply k 0 f with a factor of exp(−α f G RT ), as implemented in CANTERA (exchange_current_density_formulation in the interfaceKinetics class).

SEI film resistance and formation rate.-
The film resistance in Eq. 30 depends on the thickness of the SEI layer. Starting from the area-specific resistance of the SEI layer according to the volume-specific film resistance is obtained by dividing with the specific surface area, Inserting Eq. A35 and A42 into Eq. A43 yields In the following an explicit relationship between SEI resistance and SEI volume fraction is developed. In analogy to Eq. A33, the volume of the SEI layer on one particle is equal to which, solved for the layer thickness, yields The particle-related volumes are related to the macroscopic volume fractions via Inserting Eqs. A34 and A47 into Eq. A46 yields Finally, Eq. A48 can be cast into Eq. A44, yielding the formula of the volume fraction dependent SEI resistance, For diffusion-limited film growth, the film formation rate r is proportional to the diffusion flux (Fick's first law) of the limiting species (here: electron 36 ) with concentration c, [A50]   Eq. 10 * Units of mol, m and s depending on reaction stoichiometry, see Electrochemistry and multi-phase chemistry section.