A Coupled Electrochemical-Thermal Failure Model for Predicting the Thermal Runaway Behavior of Lithium-Ion Batteries

Thermal runaway is always a troublesome problem that hinders the safe application of high energy lithium-ion batteries. There is an urgent need to interpret the voltage and temperature changes and their underlying mechanisms during thermal runaway, in order to guide the safe design of a battery system. This paper is dedicated to building a coupled electrochemical-thermal model that can well predict the voltage drop and temperature increase during thermal runaway. The model can capture the underlying mechanism of 1 (cid:2) the capacity degradation under high temperature; 2 (cid:2) the internal short circuit caused by the thermal failure of the separator; and 3 (cid:2) the chemical reactions of the cell components that release heat under extreme temperature. The model is validated using by experimental data, therefore the modeling analysis has high ﬁdelity. We employ the model to analyze 1) the capacity degradation under extreme temperature; 2) the inﬂuence of the SEI decomposition and regeneration on the thermal runaway behavior; 3) the heat generation by internal short circuit in the thermal runaway process. The discussions presented here help extend the usage of lithium-ion batteries at extreme high temperature ( > 80 ◦ C), and guide the safe design of lithium-ion batteries with less hazard level during thermal runaway.

Lithium-ion batteries are favored to be used in today's electrochemical energy storage systems, given their extended cycle life and high energy density. 1 Higher energy density of lithium-ion batteries is required for applications with limited space, such as electric vehicles, 2 space craft, 3 smartphones, 4 etc. However, higher energy density sacrifices the safety performance, 5 which is raising more concerns in recent years. [6][7][8] The safety problems always occurs with fire, 9,10 explosion 11 and toxicity, 12 threatening the human lives and properties. A battery management system (BMS) should reduce the hazards of battery failure by monitoring the battery status. [13][14][15] The common signals that a BMS relies on for fault diagnosis are voltage and temperature. 16 A model that can simultaneously predict the voltage drop and temperature rise during battery failure will benefit the design of safer battery systems. 17 Such a model should be capable of capturing the underlying coupled electrochemical-thermal mechanisms of battery failure at extreme temperature.
The coupled electrochemical-thermal failure of a lithium-ion batteries at extreme temperature includes some key processes: 1 the capacity degradation under high temperature; 18 2 the internal short circuit (ISC) caused by the failure of separator; 19,20 3 the chemical reactions of the cell components that release heat under extreme temperature, 21,22 finally triggering thermal runaway(TR). 23-26 1 and 2 are more of electrochemical processes, and 3 is more of a chemical process. As all of the failure processes are triggered by temperature rise, the electrochemical-thermal model is highly coupled by temperature-related functions, e.g., heat generation by ISC, degradation caused by high temperature etc.
There are many available degradation models for predicting the degradation of lithium ion batteries. Empirical models have capacity as its only output, 27 which is usually predicted by empirical equations. 28 However, an empirical model does not consider correlated degradation mechanisms, e.g. loss of active materials (LAM), loss of lithium inventory (LLI), and ohmic resistance increase (ORI) etc., therefore it is hard for further establishing connections to degradation induced by high temperature. Upon this, a mechanistic model can help, because it can capture the degradation mechanism (LAM and LLI) based on the first principles. 29, 30 Christensen and Newman 31,32 first proposed a mechanistic model whose parameters directly link with the LAM, LLI * Electrochemical Society Member. z E-mail: hexm@tsinghua.edu.cn; ouymg@tsinghua.edu.cn and ORI of a lithium ion battery. The dynamics of LAM, LLI and ORI are described by Arrhenius type Equations. The growth of the solid electrolyte interface (SEI) film, which causes LLI and ORI, has been modeled in the literature [33][34][35] and relationship compatible with the mechanistic model have been identified. 36 As the SEI decomposition and regeneration 37,38 are important chemical processes pertinent to battery failure at high temperature, the modeling approaches need to be modified to describe the interactions between the electrochemicalthermal processes. ISC, from which the electric energy of the cell can be released 39,40 and the thermal runaway can be initiated, 41,42 , plays an important role in the battery failure process. ISC usually occurs when the separator fails to obstruct the contact of the cathode and anode due to various kinds of reasons, including mechanical crush, 43,44 dendrite growth, [45][46][47] separator shrinkage, 48 etc. The heat generation during ISC is critical for a coupled electrochemical-thermal model. Zavalis 49 et al. have built a detailed 2D electrochemical model that can simulate varies kinds of ISC. The distribution of ISC current and heat generation can be presented in 2D. Santhanagopalan 50 et al. have modeled four kinds of ISC, pointing out that the joule heat is critical to initiate a temperature rise. Further studies started to use an equivalent resistance to control the degree of ISC in simulation, [51][52][53] and can fit the experimental data well. 54 Recently, Coman 17 et al. pioneered the use of an efficiency factor to describe the proportion of electric energy that is converted into heat, considering the venting process, making it possible for further establishment of an electrochemical-thermal coupled battery failure model. TR, which releases energy instantaneously, is a key prediction target of an electrochemical-thermal coupled model. Dahn's group 55,56 pioneered in establishing a battery TR model using chemical kinetics that are described by Arrhenius type Equations. The kinetic parameters can be acquired from calorimetric tests. 57, 58 Spotnitz and Franklin 21 summarized the reaction kinetics for several kinds of cell materials, and demonstrated the validity of the Dahn's TR model. The Dahn's TR model has further been modified for 3D TR simulation, 59 and TR propagation within a battery pack 42,60,61 in the following years. Coman et al. 62,63 amended the venting process on the Dahn's TR model. Another recent trend is to combine the TR model with an electrochemical model together to simultaneously predict the voltage and temperature during TR. Lee 64   the influence of the design variables on the features of battery thermal runaway. Ping 66 et al. built a 2D electro-thermal coupled model to predict the TR, which is triggered by the heat generated during cycling under adiabatic conditions. Both those models can simultaneously predict the battery voltage and temperature, however, they did not consider the influence of the high temperature on the electrochemical behavior, and could not simulate the degradation caused by exposing the cell components to high temperature. In other words, they implement one-way coupling by electrochemical→thermal, rather than thermal→electrochemical. Therefore the Step 1 in the thermal failure still cannot be well modeled to the best knowledge of the authors. This paper aims to establish a coupled electrochemical-thermal model that can predict the voltage and temperature of the lithiumion battery during thermal runaway tests. Coupled electrochemicalthermal dynamics that can describe 1 the capacity degradation under high temperature; 2 the internal short circuit (ISC) caused by the thermal failure of separator; and 3 the chemical reactions of the cell components that release heat under extreme temperature, are built by specific mathematical functions with physical meanings. The voltage drop and temperature rise predicted by the coupled electrochemicalthermal model can fit well with the experimental data, whereas the predicted degradation also matches the test results. The analysis section discusses the capacity degradation under extreme high temperature exposure, the influence of the SEI regeneration on the total heat generation, and the critical conditions when an ISC triggers TR.

Theoretical
The electrochemical model.-The mechanistic model defined in Ref. 36 is used to simulate the electrochemical behavior of the lithiumion battery. The model assumes that the open-circuit-voltage of a full cell is the difference of the cathode potential U ca (y) and the anode potential U an (x), as shown in Figure 1. The output voltage of the model should be calculated with additional voltage drop caused by ohmic potential loss, and that caused by entropy change at high temperature, as calculated in Eq. 1: where V mdl is the voltage output of the model, U ca (y) is the cathode potential, U an (x) the anode potential. I is the current. I = 5A for discharge test, whereas I is the current for ISC during high temperature test, because there is no external load. R cell is the average ohmic resistance of the battery throughout the discharge process, T is the temperature, 0001V is the entropic coefficient to fit the voltage depleting for T<80 • C.
U ca (y) and U an (x) vary with the state-of-charge (SOC) of the cathode and the anode, y and x, respectively. The functional relations for U ca and U an are measured from half-cell test, as in Ref. 67 and shown in Figure 1. Note that y and x are the SOC in the half-cell test, not the y in Li y (Ni α Co β Mn γ )O 2 or the x in Li x C 6 . At the beginning of discharge, the cathode is at its low SOC (y→0) and the anode in its high SOC (x→1). As the discharge process goes on, y and x change proportionally to the current I according to Eqs. 2 and 3: where Q ca and Q an are the available capacity for the cathode and anode, respectively. η is the columbic efficiency, η = 1 especially for lithium-ion batteries. y 0 and x 0 denote the initial value for y| t = 0 and x| t = 0 in the discharge process, respectively. y end and x end denote the value for y and x at the termination of discharge, respectively. The values of {y 0 , x 0 , Q ca , Q an , R cell }, as shown in Figure 1, are identified using genetic algorithm with a cost function considering root mean square error, as in Ref. 36. Some of the parameters can be referred from Ref. 67 with minor revisions to fit the experimental data of a specific cell used in this paper. Note that there is no load across the cell during the validation tests using accelerating rate calorimetry(ARC), therefore the electrochemical-thermal effect of heat generation is not considered in the electrochemical model. The parameters {y, x, Q ca , Q an , R cell } in the electrochemical model directly link to the degradation forms of LLI, LAM and ORI. The correlated dynamics are further built as follows.
The degradation model under high temperature.-The coupled electrochemical-thermal processes related with SEI decomposition and regeneration.-In 2014, Tanaka and Bessler 38 pointed out that the decomposition and regeneration of SEI occur simultaneously at high temperature for lithium-ion cells. A similar phenomenon was discussed by Richard and Dahn 37 fifteen years ago. Here we reformulate the correlated mechanisms in Figure 2, trying to make the complex process more clear. Figure 2a shows the original state at the surface of the graphite anode, with compact SEI protecting the intercalated lithium from electrolyte.
1 SEI decomposition Figure 2b shows that once the temperature rises higher than the decomposition temperature (T onset,SEI ), the decomposition of SEI initiates with the reaction kinetics calculated by Arrhenius Equation as Eq. 4: where dc d SEI dt refers to the decomposition rate of SEI, A SEI is the preexponential factor, c SEI is the normalized concentration of SEI, E a,SEI is the activation energy, R 0 = 8.314J · mol −1 · K −1 is the molar gas constant, T is temperature in the model. The parameters related with Arrhenius Equations are collected in Table II. 2 Li reacts with electrolyte at the anode interface The intercalated lithium in the anode reacts with the electrolyte, when the SEI breaks down, as shown in Figure 2b. Such a reaction generates heat, of which the amount can be measured by DSC. 58,68 Figure 2. Coupled electrochemical-thermal processes that cause SEI decomposition and regeneration at extreme high temperature, usually higher than 60 • C.
The reaction is modeled by a modified Arrhenius Equation as in Eq. 5: where dc d anode dt means the rate for the reaction between the intercalated lithium and the electrolyte, A anode is the pre-exponential factor, c anode is the normalized concentration of full reaction between the anode and the electrolyte, E a,anode is the activation energy. The parameters related with Arrhenius Equations are listed in Table II 60 is a correlation factor used to modify the reaction rate, taking into account that the thicker the SEI layer is, the slower the reaction will be.
The reaction between the intercalated lithium and the electrolyte is like the "formation" process to form SEI layer on the surface of anode during manufacturing. The formation of new SEI layer is called as the regeneration of SEI, as shown in Figure 2c. 3

Loss of lithium inventory
The SEI regeneration at high temperature will definitely cause LLI. κ LLI , the rate of LLI, is regarded as proportional to dc d anode dt as in Eq. 6: anode dt [6] where K LLI = 2.7 is used to fit the heat generation data. Furthermore, x LLI , which is the total LLI caused by SEI regeneration can be modeled by Eq. 7: Consequently, the x 0 in Eq. 3 will decrease by x LLI from its original x 0,0 , as in Eq. 8. The net change of the c SEI , which also represents the thickness of the SEI layer, is defined as κ SEI in Eq. 10. The c SEI can thus be calculated by Eq. 11. [10] c SEI = c SEI,0 − t 0 κ SEI dτ [11] 5 Resistance growth The growth of SEI layer will cause an increase in the cell resistance. The increasing rate of the cell resistance caused by SEI regeneration is defined in Eq. 12: 36 [12] where κ R is the increasing rate of the cell resistance, and K R = 0.074 links κ R and κ LLI . Furthermore, the increment in the cell resistance ( R SEI ) can be integrated by Eq. 13.
Resistance growth of the cell.-The increase of cell resistance not only includes the increment caused by SEI growth ( R SEI ), but also contains the increment caused by separator melting ( R sep ) and electrolyte vaporization ( R ele ). The R cell can thus be modeled by Eq. 14: [14] where R cell,0 is the initial value of R cell at 25 • C, ξ(T) is function that compensates the resistance decrease caused by temperature rise. According to the test results as shown in Figure 3, ξ(T) is given by Eq. 15:  where T ξ = 523K is the constant that sets a reference temperature, b ξ = 0.83 is an offset factor that is calibrated from the linear interpolation in Figure 3. R sep and R ele have physical meanings as shown in Figure 4. R sep does not increase until the temperature reaches the melting point of the PE base (the separator of the cell has PE base with ceramic coating). The R sep drops after the temperature reaches 200 • C or higher, because the separator starts to get thinner and thinner as temperature rises, as shown in Figure 4a. R ele starts to rise at around 100 • C, because the DMC (boiling point 91 • C) in the ternary electrolyte (1M LiPF 6 & DMC:DEC:EC = 1:1:1) starts to vaporize, 69 as shown in Figure 5. The vaporization of solvent leads to swelling of the cell, thereby increasing the cell resistance. The swelling stops at around 130 • C, which is higher than the boiling point of DEC (boiling point 128 • C). R sep and R ele are set as lookup table in the model to fit the experimental results, as shown in Figure 4.
Loss of active material at the anode and cathode.-The LAM at the anode and cathode is regarded to be dominated by the leakage of the ternary electrolyte, which brings deactivation of some regions on the porous electrodes, as shown in Figure 5. The leakage of electrolyte, often occurring between 100 • C-130 • C as shown in Figure 5b, can be inferred from the temperature results, as shown in Figure 5d. The LAM of the anode and cathode are described by Eq. 16 and 17: Q an = c LAM,an · Q an,0 [17] where c LAM,ca and c LAM,an are the normalized concentration indicating the degree of LAM at high temperature for cathode and anode respectively. Both c LAM,ca and c LAM,an drop from 1 to a lower level larger than 0. Q ca,0 = 27.3Ah and Q an,0 = 32.2Ah are the original capacity of the cell before LAM. Assume that the LAM of the anode is mainly controlled by the electrolyte leakage, using a variable named c LAM,ele in the model, as in Eq. 18. The LAM of the cathode is not only controlled by c LAM,ele , but also influenced by the dissolution of transition metal, 30,70 which is marked by c diss LAM,ca as in Eq. 19. c LAM,an = c LAM,ele [18] c LAM,ca = c diss LAM,ca · c LAM,ele [19]

Variable
Update function Initial value The decreasing rate of c diss LAM,ca can be modeled by the Arrhenius type Equation in Eq. 20: where κ diss LAM,ca = dc diss LAM,ca dt is the dissolution rate, A LAM,ca is the preexponential factor, E a,LAM,ca is the activation energy. The parameters related with Arrhenius Equations are listed in Table I.
Eq. 21 is established to simulate the LAM caused by the leakage of electrolyte: where κ LAM,ele = dc LAM,ele dt is the rate of vaporization, A LAM,ele is the preexponential factor, E a,LAM,ele is the activation energy. The parameters related with Arrhenius Equations are listed in Table I. The offset 1/3 in Eq. 21 can regulate the final value of c LAM,ele to be larger than 1/3, because EC does not gasify until the temperature reaches its boiling point at 248 • C. 69 Finally, Table I lists the correlated dynamics that are used to model the LAM at the cathode and anode. Figure  6a displays an equivalent circuit representation of the cell during ISC. The short circuit current can be calculated using Eq. 22:

Internal short circuit under extreme high temperature.-
where V mdl is the voltage calculated by Eq. 1, R ISC is the equivalent resistance for the ISC. Figure 6b shows the value of R ISC that is set in the model. The variation of R ISC encounters four stages: 1) the self-discharge at high temperature between 50-140 • C; 2) the increase caused by the melting of PE base of the separator; 3) the decrease caused by the thickness decreasing of the separator; 4) hard ISC at a critical temperature T * ISC . The critical R ISC that reflects the degree of hard ISC is marked as R * ISC as shown in Figure 6b.
.01 in this model. Heat generation due to the ISC is modeled by Joule heating as shown in Eq. 23: where Q ISC is the heat generation power of ISC, ζ = 0.28 is the efficacy coefficient, which is proposed in Ref. 17 to fix the non-ohmic heat generation, I is the ISC current as in Eq. 22, R cell the resistance of the cell as in Eq. 14, R ISC is the equivalent ISC resistance as in Figure 6b.
Thermal model.-The thermal model is built to predict the temperature rise during the TR process. 60 The temperature T is modeled using Eq. 24. The derivative of T, i.e., dT dt , satisfies Eq. 25 according to the Energy Balance. In Eq. 25, Q is the total power of heat generation, M = 720g is the mass of the cell; C p = 1100J · kg −1 · K −1 is the specific heat capacity of the cell. The specific heat capacity C p was measured before TR test using ARC, and similar values can be seen in Ref. 71.
Q is the sum of all kinds of heat generation/dissipation terms, as in shown Eq. 26: [26] where Q chem represents the heat generation power by chemical reactions, Q ISC represents the heat release power generated by ISC, and Q dis represents the heat transfer/dissipation power into the environment. The ARC can provide an adiabatic environment during test, therefore we can set Q dis = 0 in the model. Q chem can be calculated using Eq. 27, where Q z denotes the power in watts dissipated as heat by an arbitrary reaction z. z can be {SEI, anode, separator, cathode,1, cathode,2, electrolyte} as listed in Table II. The decomposition of The heat generation rate Q z is proportional to the decomposition rate ( dc d z dt ) of the normalized concentration (c z ) of reactant in reaction z, as shown in Eq. 28. h z is the enthalpy of the chemical reaction z; m z is the total mass of the reactants inside the cell; T onset,z is the onset temperature of the reaction z. The determinant condition, T>T onset,z within the bracket implies that the reaction z only happens when the T is higher than T onset,z .
c z can be integrated from its derivative dcz dt , as shown in Eq. 29, where c z,0 represents the initial value of c z , κ z is a proportionality factor, dcz dt equals the difference between the decomposition rate dc d z (t) dt and the regeneration rate dc g z (t) dt of reactant z, as shown in Eq. 30.
The decomposition rate dc d z (t) dt conforms with Arrhenius Equation as shown in Eq. 31: [31] where A z is the frequency factor; n z,1 and n z,2 are the orders for reaction z; E a,z is the activation energy; R 0 = 8.314J · mol −1 · K −1 is the ideal gas constant; g z is the correction term of specific reaction. The values of the parameters that are related to chemical kinetics and have z as subscripts have been listed in Table II.
The Q Ca+An in Eq. 27 denotes the rapid oxidation-reduction reaction between the cathode and anode after a TR is triggered. The rapid oxidation-reduction reaction resembles the reaction between the fuel and oxygen in a combustion reaction. 74 The exact mechanism of the rapid oxidation-reduction reaction during TR is not clear, because the heat release rate is too fast. Figure 7 provides hypothesis for the reaction mechanism involving the rapid oxidation-reduction reaction. If the separator is intact, the decomposition reactions are confined within a limited spaces, which we call a reaction system (marked as SYS in Figure 7). The SYS ele ca , SYS ele sep , and SYS ele an react independently reflecting in the Q z in Eq. 27. We have sufficient knowledge on the  reaction kinetics of Q z , because the reaction system is well-defined, and the heat generation can be well measured using DSC. However, once the separator collapses at extremely high temperature, the cathode and anode will come into contact with each other, forming a complex reaction system comprised of SYS ele ca + SYS ele an . The reaction of SYS ele ca + SYS ele an is quite rapid as defined in Ref. 58, with an onset temperature of T onset,Ca+An ). Liu 75 et al. believe that the oxidationreduction reaction is caused by "chemical-crosstalk" with migration of oxide material across the separator, making the case even more complex. This is a new area of active researches, and here we choose to model the system based on the assumptions shown in Figure 7 and Ref. 58.
Q Ca+An can then be calculated using Eq. 32: where H Ca+An is the total energy released to heat the battery during the rapid oxidation-reduction reaction; t is the average time for energy release, here we choose t = 10s for the model based on experimental observations; and T onset,Ca+An is set as 260 • C for the cell, because the cell used in our experiments has a PE-based ceramic coated separator with a collapse temperature of 260 • C as reported in Ref. 76. The total chemical energy released during TR process, H, satisfies the energy balance Equation 33: [33] where T represents the total temperature rise measured by ARC test; H chem is the total energy released by all of the chemical reactions, which is calculated by Eq. 34; H ISC is the total thermal energy release by ISC, which is calculated using Eq. 35. H Ca+An is set to 308000J in order to fit the T in our experiments.
Experimental TR tests are conducted using ARC to validate the coupled electrochemical-thermal model. The ARC (Figure 8a) is manufactured by Thermal Hazard Technology. A common ARC test follows the heat-wait-seek method. 37,74 The ARC with extended-volume (EV-ARC) can hold the cell with large format in the tests. The 25Ah cell employed for the ARC test is a rechargeable LiNi x Co y Mn z O 2 polymer battery manufactured by AE Energy Co. Ltd. It has LiNi x Co y Mn z O 2 as its cathode and graphite as its anode. The separator is polyethylene (PE) based, with a ceramic coating. The detailed material composition of the battery can be found in Ref. 67.
The voltage and the temperature were recorded during the test, in order to validate the coupled electrochemical-thermal model. The experimental results are shown in Figure 8b. A unique test, which is called the TR test with early termination, 67 is set up to study the degradation behavior of lithium-ion battery exposed to "near" runaway conditions. In the ARC test with early termination, the heating process was stopped once the temperature reached a pre-defined value (T cool ), as shown in Figure 8b. The cell sample was pulled out of the chamber and cooled down to ambient temperature. T cool = {80 • C, 90 • C, 100 • C, 110 • C, 120 • C, 130 • C, 140 • C, 160 • C} for the ARC test with early termination. The cells were discharged at 5A (C/5 rate) after cooling to ambient temperature. The voltage and capacity of the discharge process are used to check the performance of coupled electrochemical-thermal model in simulating the capacity degradation under extreme temperature. Further details on the ARC results can be found in Ref. 67 and. 74

Results and Discussion
Model validation.-The simulation is conducted using Simulink 8.7 with code in Matlab 2016a. The solver is chosen as ode23s (stiff/Mod. Rosenbrock) with variable steps. The relative tolerance is 10 −8 during simulation. The numerical solutions of Arrhenius Equations can be referred to Ref. 21. Figure 9 shows the verification results for the coupled electrochemical-thermal model. The simulation results can fit the ARC test results for the T-t, V-t and V-T profile, as shown in Figures 9a-9c, indicating the coupled electrochemical-thermal model can predict the evolution of the voltage and the temperature well. Figure 9d shows that the electrochemical-thermal coupled model can also fit the dT/dt results well. The model has high accuracy in simulating the dT/dt results, especially for the temperatures between 80-300 • C, as shown in Figure 9e, indicating that the model can capture the dynamics of heat generation caused by chemical reactions.   Figure 10a shows the different heat generation rates that contribute to the total heat generation in Eq. 26. The proportion of the heat generation by ISC (Q ISC ) in the total heat generation (Q Total ) is similar with that was reported in Ref. 77. Figure 10b shows the normalized concentrations of different reactants. Most of the reactions occur sequentially with concentrations dropping from 1 to 0. A special case is that the c SEI,0 = 0.15 for the SEI decomposition. The increment in c SEI can be clearly seen in Figure 10b for 60 • C<T<90 • C, indicating the speed of SEI regeneration is faster than that of SEI decomposition at that time. Figure 10c presents the cumulative heat with a unit of Joule for all the heat generation terms, allowing a more rapid appreciation of the reactions that drive the temperature rise. Moreover, the percentage of individual heat generation terms in the total heat generation is summarized in Figure 10d. Interestingly, the heat generation by ISC only occupies 2% of the total heat during TR, indicating that ISC is not the major heat sources of TR.
The discharge curves of the electrochemical model fit well with that from the experiments, for all T cool between 80 • C and 130 • C, as shown in Figure 11a. For the T cool >130 • C, the cell capacity approximately drops to zero, therefore the results for T cool >130 • C are not shown. The capacity retention rate of the model also fits well with the experimental data, as shown in Figure 11b, indicating that the variables in the electrochemical model can capture the real capacity degradation mechanisms at high temperature. The critical variables {y, x, Q ca , Q an , R cell } are shown in Figures 11c-11e, respectively. It can be seen that the LLI starts from 90 • C due to SEI regeneration, whereas the LAM starts from 100 • C due to electrolyte leakage.
We add remarks here to clarify the relationship between the concentrations that are used for electrochemical decomposition and that for thermal failure. This paper considers the "thermal→electrochemical" coupling that models the influences of the high temperature exposure on the battery capacity degradation by Eqs. 16-21. The logic flow for the "thermal→electrochemical" cou-pling should be 1) c z →Q z by Eqs. 28-31; 2) Q z →T by Eqs. 24-27; 3) T→c LAM by Eqs. 18-21. Therefore, the concentrations of chemical reactions (c z ) affects the evolution of the concentrations of electrochemical decompositions (c LAM ). A close loop was not established to describe the feedback effect of c LAM on the c z , because we need more data that can help to build the relationship. Correlated issues might be interesting for further study.
The good prediction of the coupled electrochemical-thermal behavior convinces us that the model has good validity, and can be used for further analysis to discuss: 1) the capacity degradation under extreme temperature; 2) the regeneration of SEI and its influence on the TR behavior; 3) the influence of the ISC on the total heat generation during TR; and 4) parameters that may influence the critical features of TR.

Modeling analysis.-The degradation under storage at extreme
temperature.-The coupled electrochemical-thermal model helps study the degradation under storage at extreme temperatures. The cell is assumed to be suddenly exposed to an isothermal environment, like a thermal chamber, or another extreme environment. We wish to determine how long the cell can tolerate the high temperature storage, in order to explore the tolerance of the lithium-ion cells exposed to extreme conditions. The heat dissipation term, Q dis , should be modified by Eq. 36: [36] where h conv = 5 W · m −2 · K −1 is the convection coefficient for natural convection, A S is the total area of the cell, T storage is the storage temperature. T storage is set as {80 • C, 90 • C, 100 • C, 110 • C, 120 • C, 130 • C} in the simulations. The cell does not go into TR for all the T storage , because the heat generation is not sufficient at that temperature range, indicating that the cell has good thermal stability and can pass a hot box test.  Figure 12 shows the capacity degradation for the cell stored at T storage = 120 • C. The cell can maintain a capacity larger than 80% within 4.6h, as the orange dots shown in Figure 12a. If the electrolyte leakage does not occur, which means that the R ele = 0 and κ LAM,ele = 0 in the model, the capacity degradation will be delayed, and the capacity retention of 80% can last up to 10h. Figures 12b-12g compare the {y, x, Q ca , Q an , R cell } for the high temperature storage with or without electrolyte leakage. The changes in the stoichiometric coefficient y and x look similar for the two cases, as shown in Figures  12b and 12e. According to Eqs. 16 and 17, c LAM,ca and c LAM,an are the normalized Q ca and Q an . The capacity retention rate of Q ca and Q an can be reflected by c LAM,ca and c LAM,an , respectively. For the cell with electrolyte leakage, both c LAM,ca and c LAM,an decrease much faster than that without electrolyte leakage, whereas the R cell increases faster than that without electrolyte leakage. Therefore, the major degradation mechanism for the high temperature storage with electrolyte leakage is the LAM and ORI. For the cell without electrolyte leakage, the value of x in Eq. 3 drops from 0.86 to 0.72, 83.7% of its original value, whereas the c LAM,ca drops to 81.6%. Hence, the major degradation mechanism for high temperature storage without electrolyte leakage is a combination of LAM and LLI. Limited ORI also contributes to the capacity degradation without electrolyte leakage, but not as much as that with electrolyte leakage.
The duration for the capacity degradation to 80% for the cells stored at different T storage are compared in Figure 13. There is slight difference in the retention time to 80% capacity for the cell with and without electrolyte leakage, for T storage ≤100 • C, because leakage only starts at temperature higher than 100 • C for the cell used in experiments. For T storage >100 • C, the retention time to 80% capacity for cell with electrolyte decreases 100% faster than that without electrolyte leakage, indicating that avoiding electrolyte leakage can significantly improve the capacity retention time for the cell exposed to high temperatures. This modeling analysis helps understand the capacity loss mechanism at extremely high temperatures, and help extend the usage of the lithium-ion batteries at extreme working conditions.
The regeneration of SEI and its influence on the thermal runaway features.-Modeling analysis helps reveal the influence of the SEI regeneration on the TR features of lithium-ion batteries. K g SEI , as defined in Eq. 9, controls the speed of SEI regeneration. Figure 14a shows the influence of the K g SEI on dT/dt during TR. The dT/dt can maintain at a specific level by the feedback gain K g SEI , indicating stable heat generation by the SEI decomposition and regeneration reaction. 7 dT/dt ≈ 0.027 • C · min −1 for the calibrated model with K g SEI = 6. dT/dt rises from 0.022 • C · min −1 to 0.034 • C · min −1 , when K g SEI rises from 2 to 10. The magnitude of stable heat generation caused by SEI decomposition and regeneration may vary for different kinds of carbon based   anode materials, 68 therefore K g SEI is a critical parameter in the model, especially for the simulation between 60-120 • C. Figure 14b shows that a larger increment of c SEI can be observed for the model with larger K g SEI , indicating a faster SEI regeneration in the reactions. However, c SEI will always approach zero before the temperature reaches 150 • C, indicating that at T>150 • C, the rate of SEI decomposition is much faster than that of the SEI regeneration. Figure 14c shows that the TR time will be shortened for larger values ofK g SEI , because the value of dT/dt increases. Figure 14d displays the voltage profile during TR. Interestingly, the voltage for cell with smaller K g SEI will be lower than that with larger K g SEI . This is because that the self-discharge lasts longer for a cell with smaller K g SEI . Modeling analysis also helps reveal the influence of the SEI regeneration on the capacity degradation during TR. K LLI , as defined in Eq. 6, controls the rate of decrease of x, which is an important indicator of LLI. The value of K LLI has little influence on the properties of heat generation, therefore we do not compare the dT/dt curves here for different K LLI . Figure 15a illustrates the decreasing trend of x for K LLI from 2.3 to 3.1. A larger K LLI will definitely accelerates the LLI at the anode. The downward slope at approximately 250 • C indicates an ISC after the collapse of separator. Figure 15b shows the difference in the voltage drop for the model with different K LLI . A larger K LLI brings faster consumption of lithium inventory, leading to a faster voltage drop. Therefore K LLI is a critical parameter for predicting the LLI in the model during thermal runaway.  Figure 16a shows the triggering time of TR can be significantly decreased for the case with T * ISC = 135 • C. While the trigger time of TR are quite close for other T * ISC , because the rate of temperature rise at high temperatures is already high. The dT/dt curves in Figure 16b clearly show the influence of the ISC on the heat generation rate during TR (note that dT/dt is proportional to the total heat generation power Q, as shown in Eq. 25. The speed up effect on the dT/dt curves for T * ISC = 135 • C is obvious, compared to that for T * ISC = 210 • C. H ISC , which is the total heat released during an ISC as defined in Eq. 35, has been marked in Figure 16b. T ISC is in proportion with H ISC as in Eq. 37: [37] Figure 16b displays that the lower T * ISC is, the larger H ISC will be, and the higher temperature rise T ISC will be, therefore delaying the ISC to higher temperature is useful to mitigating battery TR. Our observations in the modeling analysis not only match the experimental results in our recent paper, 78 but also helps interpret the mechanisms of failure for some data sets with different T * ISC . Figures 16c-16f provides some ARC data that can validate the proposed coupled electrochemical-thermal model. Figures 16c and  16d are the ARC test results for the 25Ah lithium-ion cell with same materials, but manufactured in a different batch. The ceramic coating on the separator of the cell is not uniform, therefore an ISC will occur at 135 • C, when the shrinkage starts for the PE base. All the T-t, V-t, dT/dt curves look similar with the simulation case with T * ISC = 135 o C, as in Figures 16a and 16b. The rapid heat release during TR does not occur until the temperature rises to 250 • C or higher. Figures 16e and 16f provide another case for a 3.3Ah pouch cell with similar cell chemistries. The ISC occurs around 205 • C, giving rise to a small peak in the dT/dt curve (circled in Figure 16f), similar to that shown in Figure 16b. Figures 16c-16f further convince us that the coupled electrochemical-thermal model has high fidelity in capturing the underlying mechanisms during TR.
The H ISC = 22260J for T * ISC = 135 • C only occupies approximately 7% of the total electric energy that is charged into the cell. Since the efficacy coefficient for electric energy release is set as ζ = 28% in Eq. 23, these results imply that the electric energy has not been fully released during ISC. The low rate of energy release during ISC may be attributed to the capacity degradation caused by electrolyte leakage. Hence, simulations were conducted to compare the ISC behavior for a cell with and without electrolyte leakage. The results are shown in Figure 17. Although the ISC occurs at similar times, TR oc-curs immediately after ISC is triggered for the cell without electrolyte leakage, as shown in Figure 17a. The maximum temperature during TR rises from 869.1 • C to 907.5 • C (nearly 40 • C increase), for the cell without electrolyte leakage. Figure 17b shows that the ISC can last longer for the cell without electrolyte leakage, therefore the H ISC increases to 53443J (16% of the charged electrical energy). The ISC does not end until the rapid oxidation-reduction occurs, therefore there is still some electric energy left in the cell. The ISC process is forced to be shut down in the model, because the rapid oxidation-reduction will break the circuit that can forms stable ISC current. In summary, the cell becomes more dangerous if there is no electrolyte leakage before an ISC occurs. From this perspective, releasing the electrolyte in time before the occurrence of an ISC will help reduce the hazard during TR. However, the exhaust flammable electrolyte may be ignited outside the cell, arousing other safety problems, and will be discussed in our future work. 79 The influence of the degree of ISC (R * ISC ) on the TR behavior has been discussed by modeling analysis. The temperature curves (T-t and dT/dt) look similar for R * ISC = 0.01 and R * ISC = 0.001 , as shown in Figure 18, indicating the heat generation rate for different R * ISC can be similar. Figure 18b reveals the mechanism underlying this phenomenon: because a smaller R * ISC leads to lower level of voltage after an ISC is triggered. A lower level of voltage denotes a lower level of the heat generation power Q ISC . Hence the R * ISC has little influence on the TR behavior, as long as it is sufficiently small to reflect a hard ISC.
Moreover, we try to study the influence of the efficacy coefficient ζ on the TR behavior using the model. Figure 19a shows that if the value of ζ increases to 100%, the TR will be much severe, and the TR will be triggered much more rapid once T reaches T * ISC = 135 • C. The maximum temperature during TR will further increase from 907.5 • C to 937.5 • C, due to the increase in ζ. Figure 19b further confirms that a more severe hazard might be brought by the increase in ζ. In conclusion, the parameter ζ has a big influence on the TR behavior of lithium-ion cells. The accurate measurement of ζ still requires further study. Ref. 17 suggests that ζ should include the effects of venting, which is also interesting for further investigation.
Finally, the results of modeling analysis have been concluded by Figure 20. The lower value of the T * ISC is, the larger influence will be on the TR features, because the capacity degradation will significantly reduce the cell capacity at higher temperatures. The leakage of electrolyte can cause up to 40 • C difference in the maximum temperature at TR for the cell used in this study, whereas the ζ can cause 30 • C     difference. Figure 20 tells us that if the degree of ISC can be controlled by the electrolyte concentration and the heat release efficiency of joule heating, we can reduce the maximum temperature of the cell by up to 100 • C. According to Figures 17-19, the occurrence of an ISC is a critical factor to trigger TR, therefore managing the triggering temperature of ISC will help postpone the TR to higher level. However, the total heat generation by ISC is not the major heat source that heat the cell to 800 • C or higher. Figure 21 summarizes the percentage of the heat generation by ISC in the total heat generation of TR, for R * ISC = 0.01 ,ς = 0.28, no matter there is electrolyte leakage or not. The percentage of the heat generation by ISC rises when T * ISC decreases, indicating the triggering condition might be severer if massive ISC occurs at lower temperature. Less than 10% of the total heat is released by ISC, confirmed by our modeling analysis. The thermal energy released by chemical reactions covers the majority of the total energy released during TR. The rapid oxidation-reduction reaction is the main reaction that determines the maximum temperature, and the mechanism requires further study.
Effect of critical temperatures on the thermal runaway features.-Further modeling analysis are conducted to discuss the critical temperatures that may significantly influence the TR features in this section. The collapse of the separator sets a key tipping point for the energy release style during TR, as shown in Figure 7. Before the collapse of separator, the reactions occur within their own systems. However, after the collapse of separator, the cathode (oxidizer) and anode (re- Figure 21. The percentage of the heat generation by ISC in the total heat generation of TR, R * ISC = 0.01 ,ς = 0.28. ducer) have chances to electrically contact each other, and rapid reactions will start at high temperatures. The critical temperatures that represent the critical point during TR include: 1) T sep,mlt , the melting point of the base of the separator, T sep,mlt = 120∼140 • C for PE, and T sep,mlt = 160∼180 • C for PP; 2) T sep,clp , the collapse temperature of the separator, T sep,clp = T sep,mlt for a single-layer separator, while for a multilayer separator T sep,clp can be up to 200 • C or higher; 3) T * ISC , the temperature for ISC, to the best knowledge of the authors, most of the test results points that T * ISC = T sep,clp . However, under the cases when the cell is under mechanical abuse, T * ISC can be lower than T sep,clp ; 4) T onset,Ca+An , the onset temperature of the rapid oxidation-reduction reaction between the cathode and anode, as shown in Figure 7 and Figures 22d and 22e.
The critical behavior of TR can be interpreted if the relationship between the four critical temperatures is clear. Generally, the relationship among T sep,mlt , T sep,clp and T * ISC is clear, as shown in Figures  22a-22c. However, the relationship between T onset,Ca+An and T * ISC is still questionable. In Ref. 73 the authors believe that T onset,Ca+An can be lower than T * ISC , because the oxygen and other oxidizers that are released by the cathode can migrate from the cathode to the anode, trigging TR before ISC, as shown in Figure 22d, here we called Case 1 for the style of TR. On the contrary, in this paper, we have T onset,Ca+An >T * ISC , which is called as Case 2, as shown in Figure 22e. Case 2 has two branches, one is that the heat released by ISC can further trigger TR, which is called the Case 2.1, whereas another one is that the heat released by ISC cannot further trigger TR, which is called the Case 2.2. We emphasize Case 2.2, because if the conditions that lead to Case 2.2 are well understood, then chances that the TR is triggered by ISC will be minimized.
Modeling analysis are conducted for the Case 1, Case 2.1 and Case 2.2, respectively. In Case 1, the T onset,Ca+An = 130 • C, referring to the onset temperature for oxygen release as in Ref. 80, and T * ISC = 170 • C>T onset,Ca+An . In Case 2.1, T * ISC = 170 • C<T onset,Ca+An = 260 • C, whereas in Case 2.2, T * ISC = 170 • C<T onset,Ca+An = 600 • C. Figure 23 displays the simulation results for the three cases. Figure 23a shows that the TR is triggered immediately when the temperature reaches T onset,Ca+An = 130 • C for Case 1, there will be no ISC because the circuit for short has been broken in advance. Interestingly, there is less TR for Case 2.2, although a fast temperature rise is still triggered at T * ISC = 170 • C, similar to Case 2.1. The maximum temperature for Case 2.2 is 472.6 • C, which is much lower than that for Case 2.1 (861.6 • C), because the thermal energy of the rapid oxidation-reduction reaction has not been released. Eq. 38 defines the condition that when the fierce oxidation-reduction reaction will not be initiated: coupling makes the modeling analysis on the heat generation by ISC more reliable, and proving that two-way coupling based on physical equations is possible in TR simulations. 2. New physical based functions are established, e.g. Eqs. 6, 9, 10, 14, although some of the parameters are still empirical. Those physical based equations can guide further research on quantify specific side reactions with parameters involved in two-way coupling. 3. Successfully quantify the proportion of the ISC in the total heat generation during TR, with proper experimental validation. The quantified results indicate that the ISC is not the major heat source during TR. Maybe one day we would observe that there will be battery TR without ISC, as predicted by the modeling analysis.
We also give our suggestions on the future work, to further improve the capability of the electrochemical-thermal coupled battery failure model. We encourage further research to focus on the "electrochemical-thermal" way of coupling, which can solidify the two-way coupling of the model, making it closer to reflect the real physical/chemical processes during TR. And next significant progress may be fulfilled in the simulation considering the influence of battery state-of-charge on the TR behavior: 1. More functions that can describe the coupling way of "electrochemical→thermal" should be built, in order to simulate the influence of battery state-of-charge on the TR behaviors. Based on some unsystematic results, we believe that the proposed ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 207.241.231.83 Downloaded on 2019-04-28 to IP model may be easily extended for cells with state-of-charge larger than 75%. However, more data are required to come to a more solid conclusion. 2. The influence of state-of-charge on the TR behavior should be established by setting new equations between c LAM in the electrochemical model and c ca , can in the thermal model. More data are needed to quantify the electrochemical-thermal coupled relationship.

Summary
A coupled electrochemical-thermal TR model for lithium-ion batteries has been built to simultaneously predict the voltage and temperature during TR. The model can fit the experimental data well; it captures well the underlying dynamics of 1 the capacity degradation at high temperatures; 2 the ISC caused by the thermal failure of separator; and 3 the chemical reactions of the cell components that release heat under extreme temperature. The modeling analysis that discusses capacity degradation under extreme temperatures can be used to extend the usage of lithium-ion batteries at extreme high temperature. Our modeling analysis also studied the influence of the SEI decomposition and regeneration on the TR behavior. The discussion on the critical parameters relating with the ISC process benefits the safety design of lithium-ion battery with lower hazard levels during TR.