Porous Electrode Model with Particle Stress Effects for Li(Ni1/3Co1/3Mn1/3)O2 Electrode


 Inaccuracies in some equations have been identified. The corrections are listed below. Note that the errors do not affect any results. The conclusions are still valid.

Lithium transition metal oxides are currently the focus of battery manufacturers as attractive positive intercalation electrode materials in lithium-ion batteries (LIB). [1][2][3] After LiCoO 2 (LCO) was firstly introduced by Goodenough, 4 its commercial success led to the development and research studies of other metal oxides like LiNiO 2 (LNO), LiNi 0.8 Co 0.15 Al 0.05 O 2 (NCA), LiMnO 2 (LMO), Li(Ni 0.5 Mn 0.5 )O 2 (NMO) and Li(Ni x Mn y Co z )O 2 (NMC). 2 Among them, NMC materials have been under the spotlight because of its potential for high energy density applications in portable electronic devices as well as in the expanding market of electric vehicles. 5,6 A determining factor in the electrochemical performance of LIB is the lithium insertion mechanism within the metal oxides of the positive electrode. In order to elucidate the internal processes of ion migration, diffusion as well as intercalation kinetics, physics-based models like the volume-averaged model 7-10 and the single particle model 11 have been largely preferred. A physics-based model offers not only a better understanding of the internal processes; it also gives the possibility to optimize material properties or operating condition. Nevertheless, model parameterization and the choice of the physical parameters remain a big challenge. Due to the lack of experimental data or manufacturer's specification, most works rely on intrinsic transport and kinetic parameters from literature sources, or guessed values.
One reliable way to obtain and quantify the diffusional and charge transfer properties is to combine a physics-based model with experimental study like galvanostatic intermittent titration technique (GITT). To characterize porous electrode kinetics, different model techniques can be employed. For example, the use of GITT analytical theory 12 has been largely utilized. 10,[13][14][15] However, care must be taken with the assumption of uniform current distribution and the definition of the electrochemically active area. 7,13 The single particle model, on the other hand, allows the investigation of the diffusion of lithium throughout the particle bulk but neglects the porosity of the particle agglomerate which consists of many primary particles. 11 A particular phenomenon reported in using GITT coupled with a physics-based model to study porous electrode is the long time constant during relaxation. 9,16 GITT model studies of lithium manganese-rich NMC positive electrode from Dees et al. 9 showed that the material did not relax to full equilibrium. The reason for that is the difference in the equilibrium potentials of the two active electrode components giving rise to a more dynamic (de)lithiation interaction. 9 For LiNi 0.5044 Co 0.1986 Mn 0.2970 O 2 z E-mail: jyko@kth.se (NCM523) material, Verma et al. 13 proposed a model formulation for the active area for primary/secondary particle agglomerate materials to simulate GITT and half-cell performance. As were shown in these studies, the different physical behaviors exhibited by different electrode materials need to be captured by the model to accurately extract mass transport and kinetic parameters.
Although several modifications have been done to the physicsbased models coupled with GITT, none of the modeling studies has, to our knowledge, included stresses in electrode materials. However, many authors have reported stress-related behavior acting on electrode particles. Particle-level strains may be induced by both external loads during production and internal electrochemical processes during cycling. 17 During charging or discharging, lithium is removed (deintercalation) or inserted into (intercalation) the host lattices of the active material particles, causing lattice parameter variation, resulting in expansion or contraction of the particles. 18,19 Such volumetrical strains can also impose stresses on the surrounding active, filler or carbon particles. 20 During intercalation, LiMn 2 O 4 gives a 6.5% volume change due to a change in lattice parameters. 21 For NMC materials, lithium slab space increases with delithiation 14,22 and the lattice parameter c is shown to increase by 1.5%; 23 whereas the lattice parameter c increases by 2.5% for LCO. 23 As a result, stresses arise in the electrode materials. 17,21,24 An effect of stress generation is its influence on the electrochemical potential of the solid species participating in the electrode reaction, and thus on the open-circuit potential of an electrode, 20,25-27 the lithiation and delithiation kinetics, and lithium transport. A further issue which has not been widely addressed is the dynamic evolution of stress-related (chemo)mechanical properties of electrode particles at different state-of-charge (SOC). Due to the scarcity of data available for complete mechanical characterization throughout a wide SOC range, stress behavior during battery cycling has been mostly investigated assuming constant mechanical properties, disregarding of the dependency on the lithium concentration.
In this study, a physics-based porous electrode stress model is developed to investigate the GITT potential responses of Li(Ni 1/3 Mn 1/3 Co 1/3 )O 2 (NMC111) during charging at the beginning of life (BOL). It is an extended pseudo-two-dimensional (P2D) model based on Newman et al. 28 to extract transport, kinetic, stress and thermodynamic properties across different electrode voltages and temperature. In particular, a stress proportionality factor which includes the electrochemo-mechanical properties is derived, stemming from tensile stress that occurs on the electrode particle surface due to volume change, impacting the free energy for lithium intercalation. 29 The stress affects not only the equilibrium potential, but it also drives the lithium diffusion in the solid bulk. In this work, we focus on analyzing both the processes of polarization to the end of relaxation during a whole GITT transient and on examining the interfacial behavior responsible for the long times needed for NMC111 to relax to its equilibrium potential during GITT pulses. In contrast to the standard model based on the works on Newman et al. 28

Experimental
Materials.-All measurements were performed in a two-electrode pouch cell set-up, half-cell configuration (NMC111/Li metal). The separator was Whatman glass fibers microfilter, grade GF/A. The electrolyte was 1 M LiPF6 in EC:DEC 1:1% in volume (BASF, Selectilyte LP40). All cells were assembled in an inert-atmosphere glove box (LC-1, LC Technology Solutions Inc.). NMC111 electrodes were purchased from ElectrodesandMore, and the specifications are reported in Table I. The electrodes were used as received. Figure 1 shows the scanning electron microscope (SEM) of pristine NMC111 (x3.0k magnification), collected by means of a Hitachi S-4800 tabletop SEM.
Electrochemical testing.-Formation.-The Li/NMC111 cells were cycled between +2.7 V and +4.3 V (vs. Li/Li + ) in Constant Current (CC) -Constant Voltage (CV) mode for 4 cycles. 28 mA/g NMC were applied for CC charge and discharge, and 14 mA/g NMC was used as cutoff criterion during CV (at the end of charge only, i.e. +4.3 V). Formation was performed using either Solartron Analytical 1287A (Ametek) or Biologic multi-channel VMP-3 (BioLogic). The tem-  GITT.-GITT was performed after formation. Two low-current pulses were used, namely +0.064 mA (8.15 mA g −1 NMC or C/20) and +0.13 mA (16.6 mA g −1 NMC or C/10). Between each pulse, relaxation was allowed for 3 and 2 hours respectively, to allow the potential to relax to its steady-state value. The pulses were repeated until the potential cutoff of +4.3 V (vs. Li/Li + ) was reached. GITT measurements were performed in a temperature chamber at 10°C, 25°C and 40°C by means of a Solartron Analytical 1287A (Ametek).

Model Development
The porous electrode stress model is an extended pseudo-twodimensional (P2D) physicochemical model based on Newman and Tiedemann. 28 The model is developed to fit each GITT charging response of the NMC111 porous electrode based on linearization of equilibrium potential around the initial lithium concentration before a low-current pulse of +0.064 mA is applied and the incorporation of particle stress effect. The P2D model, as depicted in Figure 2, has one dimension extending in the cross-cell direction between a lithium metal counter electrode and the positive current collector, and another pseudo dimension extending into the solid active material particles. At the 1D counter electrode-to-current collector level, two domains are used, representing the NMC porous composite electrode and the separator, with the metallic planar lithium counter electrode being represented by boundary condition at the outer separator boundary.
Concentrated solution theory is used to describe the ion transport through the electrolyte. The volume-averaged transport equations account for the porous electrode and separator. The governing equations and the corresponding boundary conditions are presented in Table II. The electrolyte mass and charge transport for concentrated binary electrolyte are described by Equation 1 and Equation 2. The current expressions in Equation 3 account for the overall charge balance. Equation 4 and Equation 5 describe the solid phase charge balance in the electronically conducting matrix. The pseudo dimension at the particle level describes the mass transport or solid lithium diffusion across the radius of spherical active material particles (Equation 6). The two dimensions are coupled at the surface of the intercalation particles through Butler-Volmer kinetics (Equation 7). The model takes into account the non-uniform current distribution across the electrode thickness and assumes that the active material is in single solid solution phase. The model input parameters are listed in Table III. The electrolyte properties were taken from the work of Lundgren et al. 30 The Newman-based model in Table II is extended to study the stress effect during charging based on the analogy between thermal expansion in solid mechanics and diffusion-induced stress (DIS). 17,21,[31][32][33] For stressed solids, the equilibrium interfacial potential depends on the entropy contribution, taking into account the concentration effect, Table II. Newman-based model equations for lithium-NMC111 electrode.

Equations Equation No. Boundary conditions
Porous Electrode and Separator level: Electrolyte mass transport [5] Particle level: Bulk active material mass transport [9] and on the stress experienced by the solid lattice, 26,27 according to where c s is the solid lithium concentration. The first term on the right hand side, E conc , denotes the equilibrium potential in the stress-free state, while the second term denotes the equilibrium potential change under stress, where μ(σ) is the change in chemical potential due to stress σ in the particle. [11] where E 0 ref is constant, a θ is the activity of vacancy in the lattice and a Liθ is the activity of lithium in the lattice. During the seven GITT perturbations done in the range of 3.7−4.1 V vs Li/Li + at each temperature of 10°C, 25°C and 40°C, a change in lithium activity and thus its concentration is induced, causing a shift in the equilibrium potential from each of the initial equilibrium potential E 0 . E 0 is obtained experimentally from GITT and used as the initial reference potential   in a stress-and gradient-free state before each perturbation. Figure 3 shows the experimental E 0 plotted against the lithium fraction inside the active material, x pos , each time before perturbation . x pos is defined as c s,s /c s,max , where c s,s is the lithium concentration on the particle surface and c s,max is the maximum solid bulk concentration. The equilibrium potential which changes along the GITT transient in the stress-free state can then be described as [12] where x pos,0 is the initial lithium fraction before the current pulse is applied corresponding to each E 0. Assuming a substantially small change of E conc with x pos curve, the change in E conc with x pos can be linearized. Thus, Equation 12 is approximated as [14] where dE conc /dx pos emulates the slope of the OCP curve. Note that dE conc /dx pos is a fitting parameter.

E conc is given by the Nernst equation as
To derive the expression for σ, we consider the stress developed on the surface of a spherical particle of radius R which affects the free energy for lithium intercalation, as shown in Equation 10. The conventional stress-strain relationships in spherical coordinates, 17,21,31,32 for the radial and tangential components, are where r , θ , σ r and σ θ are the strains and stresses in the radial and tangential direction, E is the Young's modulus, ν is the Poisson's ratio and the partial molar volume. A further assumption is that the elastic properties are constant for the small range of lithium concentration change during each GITT perturbation and relaxation. Due to spherical symmetry, the radial and tangential strains, in the infinitesimal formulation of deformation are given by where u is the radial displacement. Since species diffusion in solids is a much slower process than the rate of elastic deformation, mechanical equilibrium is treated as a static equilibrium problem. 21,34 The equation for static mechanical equilibrium in the bulk of a sphere is given by Combining the four equations above with the condition of the stress to remain finite at r = 0 results in the following expressions for the tangential stress 33 where c s,ave (r) is the average concentration in the spherical volume of radius r within the particle of radius R, defined as: The radial stress σ r (R) is zero following the conventional boundary condition at free surfaces. 17,21,31,32 The tangential stress on the particle surface, r = R, takes the form The stress-dependent shift in chemical potential is expressed as As the and solid-mechanics parameters (E and ν) are assumed to be constant over each GITT pulse and the concentration difference c s,ave (t)-c s,s (t) is normalized by the maximum bulk concentration c s,max , the second term in Equation 10 results in [23] where the stress proportionality factor γ stress = E 2 3(1−v) · cs,max F includes the electrochemo-mechanical properties and is a fitting parameter.
The modified equilibrium potential for a solid active material under stress, E eq , in Equation 10 during each GITT transient now takes the form  [24] with c s,0 = c s,s (R,0). Hence, E eq does not only take into account the contribution of the local concentration changes on the particle surface but also the average concentration changes of the particle. At the same time, it depends on the thermodynamic as well as the electrochemomechanical properties. As lithium is being removed from the particle during the GITT charging process, the tangential stress is tensile at the particle surface (positive by definition). Therefore, the E eq is expected to be reduced under the effect of tensile stress. 35 Stress interacts with diffusion in a two-way coupling in which diffusion induces stress (DIS), which in turn affects diffusion (stressenhanced diffusion). 17,21,[31][32][33]36 In this study, the fitted effective intercalation diffusion coefficient D s,eff described by the Fick's second law in Equation 6 includes the impact of gradients of stress as a driving force for diffusion inside the particle: Equation 25 is taken from the derivation of a Fickian-type flux equation based on the coupling of concentration and stress presented in Appendix. The D s,eff is assumed constant inside the particle as the particle concentration and stress gradients are not expected to change significantly. As shown in Equations 24 and 25, the stress impacts both the solid lithium transport and on the insertion kinetics on the particle surface.
The number of parameters extracted by fitting the GITT transients are kept as few as possible to restrict the number of degrees of freedoms. In particular, four parameters are being optimized: the intercalation diffusion coefficient D s,eff , the kinetic rate constant for the electrochemical reaction k, the slope of the OCP curve dE conc /dx pos and the stress proportionality factor ϒ stress . All four parameters are determined for one full GITT transient, i.e. from the start of applying a polarization to the end of the relaxation.
The coupled time-dependent non-linear partial differential equations are solved in COMSOL Multiphysics 5.3.1 and parameter estimation is carried out using the Levenberg-Marquardt optimization solver. The simulated potential response was fitted to the experimental galvanostatic pulses in the GITT experiments. Model optimization and parameterization are carried out by employing the least square fitting method, i.e. by minimizing the following function χ: where α represents the fitting parameters, ϕ s,exp and ϕ s,sim are the experimental and simulated electrode potential values. Figure 4 shows examples of typical GITT potential responses. At 0 min, i.e. as soon as a low-current pulse of +0.064 mA (C/20) is applied, a potential jump can be seen which is then followed by a slower potential increase. The instantaneous response of the potential corresponds to the ohmic drop of the cell, and at this time-scale the initial polarization. Thereafter, the ten-minute non-instantaneous response of potential increase indicates the process as lithium is being extracted from the electrode, slowly building up a concentration gradient in the electrolyte as well as within the solid bulk particles. When the applied current is interrupted after 10 min, the cell is allowed to relax back to the equilibrium condition (OCP), thus the potential decays gradually, corresponding to the relaxation of the concentration gradients. The slowly changing voltage is caused by the slow diffusion processes in the cell. During this process, the concentration of lithium within the particles is gradually redistributed until the concentration of lithium becomes uniform and reaches an equilibrium state.

Results and Discussion
Firstly we compare the fitted experimental GITT potential response with both the standard Newman porous electrode model 28 and the extended model. One example of the comparison for a GITT pulse at around 3.73 V and 3.91 V at 25°C is shown in Figures 4a and 4b. Referring to the example for the GITT transient curve fits at 3.73 V in Figure 4a, both the standard and the extended models are able to fit the GITT experimental data well. However, at a higher voltage of 3.91 V in Figure 4b, discrepancies can be clearly seen between the two models. As can be seen, the slow relaxation behavior right after the pulse is stopped cannot be captured properly by the standard Newman model. The reason for this is that the lithium concentration gradient relaxes too quickly. Thus, there is a need to use an extended model for the fitting of the GITT curves. As shown in both Figures 4a and 4b, both the potential curves can be fitted relatively well with the extended model, justifying the need of including an additional parameter to compare the intercalation behavior during polarization and relaxation for the wide range of voltage studied here. Figure 5 shows every fitted GITT potential experimental response at 10°C, 25°C and 40°C. All the GITT curves have been fitted with the same extended stress model to create a consistent and comparable parameter set. Hence, the capability of describing and comparing the same physical processes can be carried out at every electrode potential with a uniform set of four physical parameters (D s,eff , k, dE conc /dx pos , ϒ stress ). The extended model fits the experimental data well at all potentials and temperatures. For all three temperatures, the behavior of a slow relaxation has been observed, predominantly at higher poten- tials. The slow relaxation behavior in NMC111 is very similar to the GITT voltage relaxation profiles reported by Dees et al. 9 for a blended electrode due to the active material's multiphase region. However, NMC111 behaves as an ideal solid solution in the standard operating voltage window 37 and is therefore not expected to exhibit such multiphase behavior. For higher temperatures in combination with higher voltages, the occurrence of the slow relaxation to equilibrium potentials during the three-hour open-circuit stand becomes more apparent. The model predicts that the relaxation to the equilibrium potential can take much longer than the allocated three hours in this work, especially when the potentials are higher. The long-term relaxation process during open-circuit stand is in alignment with the findings by Gerschler et al. 38 where they measured the OCP for up to 100 hours. The observed behavior that the voltage equalizes slower at higher temperatures is also in agreement with their work. 38 While it is common to carry out OCP measurements before any battery cell testing, allowing the cell to rest for up to tens or hundreds of hours at each potential cutoff is exhaustive and seldom done. However, these results indicate that one may have to wait long times in order to be sure that the voltage equilibrium is reached, particularly when the potential or temperature is higher.
In our extended model, we consider solely the effect caused by the nature of NMC111 material and exclude any side reactions that might occur. Factors such as the degradation effects at higher temperature, potentials, and the influence of counter electrode are neglected. In the experiments, the GITT potential responses of the NMC111 electrode were measured versus a lithium metal electrode, acting both acts as counter and reference electrode. We are aware that for this setup there may occur relaxation phenomena stemming from the lithium counter electrode. 7,11 However, contributions from the counter electrode to the relaxation of the measured potential do not seem likely in this study, since such response behavior should have been seen for all potentials, and not only for certain potentials which we observed here. It is further considered unlikely that the slower relaxation detected at higher potentials, between approximately 3.8 V and 4 V stems from increased side reactions as the electrode potentials still remains relatively low.
Along the GITT potential responses, the polarization and relaxation dynamics involve two mechanisms at two levels simultaneously: a) the inter-particle kinetic lithium exchange to reach a common equilibrium potential at the NMC particle surfaces and b) the intra-particle lithium diffusion to flatten out the concentration gradients. The processes of polarization and equilibration are described by the effective solid lithium transport (D s,eff ), electrochemical kinetic rate (k), material thermodynamic (dE conc /dx pos ) and surface stress (ϒ stress ) parameters. These parameters are extracted and investigated across the voltage range at different temperatures. Figures 6a and 6b show the D s,eff and D as a function of potentials at 10°C, 25°C and 40°C. From Equation 25, the approximation of the lithium diffusion coefficient D can be obtained. As can be seen in Figure 6, D s,eff and D are generally higher when lithium fraction x pos ≈ 0.94 to x pos ≈ 0.85 or in the potential region of ∼3.70 V to 3.75 V during initial state of charging. They decreases thereafter when the charging proceeds. As lithium is further depleted from x pos ≈ 0.80 onwards or when the potential continues increasing from ∼3.76 V, D s,eff appears to be quite steady at a value of ∼10 −14 cm 2 /s as seen in Figure 6a. As shown in Figure 6b, D continues decreasing when lithium is further depleted from x pos ≈ 0.80 onwards. Faster diffusion at lower potentials leads to a relatively faster voltage relaxation, which can be seen in Figure 4 and Figure 5. As the potential increases, the diffusion becomes slower, which contributes to the slower voltage relaxation discussed earlier. In our model, the lithium transport is governed by gradients of the chemical potential gradient, i.e. gradients in the concentration as well as the mechanical stress. Hence, D s,eff not only depends on T and SOC but also on the chemo-mechanical properties. 39 Slower diffusion indicates influence from higher Young's Modulus E and/or partial molar volume , particularly for GITT at higher voltages. From Figure 6b, D represents the diffusivity driven by the concentration gradient, but dismissing the stress effect. The dependence of D with x pos is larger across the delithiation process compared to D s,eff , implying that stress becomes more apparent at higher voltages, enhancing the lithium transport. D s,eff varies minimally with temperature, which agrees to Cui et al. 14 showing less than an order of magnitude change across a vast range of temperatures. We note that although it is more common to report a tendency of faster diffusion with increasing temperature, 14,40 the opposite behavior of slightly lower D s,eff at higher temperatures is not previously unreported in scientific literature. 15 Overall, D s,eff at these three temperatures expanding across a voltage window of operation from 3.7-4.1 V vary within the order of magnitude of around 10 −14 -10 −13 cm 2 /s.  15 This could be attributed to a lot of reasons, among them are the different type of experimental and/or techniques and quantitative means of investigating the solid-phase diffusion. Besides, computation of diffusivity also depends on the assumed diffusion length l. Most work relies on the assumption of diffusion occurring in the secondary particles, 11,14,42 while we assume it to be in the primary particles, whereas others adopted the electrode thickness 41 as characteristic length. The assumption of adopting primary particles is based on the fact that the primary particles are well separated from each other and do not appear to form dense agglomerates, as shown in the SEM image in Figure 1. Due to the different scale of diffusion lengths involved, a comparison using the characteristic diffusion time constant τ would be more reliable and meaningful. In this study, τ lies around 2.5·10 4 -2.5·10 5 s. This is similar to ∼2.5·10 5 s in the work of Cui et al. 14 who used analytical model for the GITT data but lower than 2.25·10 2 -2.25·10 4 s in the work of Wu et al. 11 who adopted single particle model to describe open circuit relaxation. Apart from the choice of diffusion length, another concern lies in how the parameterization in each study is done. Having that said, direct comparison of either the D s,eff or time constant extracted using different techniques can be rendered meaningless if one does not make clear of the underlying model assumptions and limitations. The main distinction of our model with the other models is the coupling of a porous electrode model to analyze the physical phenomena at the electrode and particle levels. Besides, the whole GITT pulse, i.e. from polarization to relaxation, is being fitted by the developed model so that both the electrochemical reaction and diffusion are accounted for simultaneously.
The fitted electrochemical kinetic rate constant k from the model in Figure 7 shows a continuous increase as x pos decreases from 0.95 to 0.55 (or as potential increases from 3.7 V to 4.1 V), resulting from a lower activation energy barrier for the lithium deintercalation reaction as the charging process proceeds. As temperature increases from 10°C to 25°C and 40°C, the kinetic rate constants gradually increases. The data are fitted to a polynomial equation (Equation 27) using temperature-independent coefficients.
The fitted isothermal polynomial equations are shown in black line in Figure 7 and the coefficients for the polynomial fit in Table IV. The activation energies obtained through the Arrhenius equation k = C exp(−E a /RT) where C is a constant shows an increase of E a by a factor of around 30% when x pos drops from ∼0.95 to ∼0.55. Figures 8a and 8b show how dE conc /dx pos and ϒ stress develop with x pos . Both parameters are fitted here as a local average quantity to approximate the extent of the equilibrium potential variation with the particle lithium concentration. dE conc /dx pos remains relatively low from around 0.2 to 0.6 for x pos = 0.95 to x pos = 0.73 or at voltages 3.70 V to ∼3.79 V. As charging continues, dE conc /dx pos increases significantly. Almost similar values can be seen for temperatures 10°C, 25°C and 40°C, indicating that we would expect similar SOC dependencies of  the OCP curves for the three temperatures. ϒ stress is from around 0 to 0.5 for x pos = 0.95 to 0.73 or at voltages 3.70 V to ∼3.79 V. As charging proceeds further from 3.80 V, ϒ stress increases rapidly. ϒ stress at 25°C and 40°C are quite comparable whereas it appears to be lower for 10°C. As ϒ stress is coupled to the electrochemo-mechanical properties, larger values which can be seen at higher voltages imply an increase of the stiffness of the material, or increased expansion upon lithium intercalation-deintercalation, thus contributing to the buildup of surface stress on the NMC particle. Stress has been associated with solid particle volume change which is linked to variation in lattice parameters. Therefore, the increase in stress-related properties could be influenced by an increase in (001) plane distance, which is the total thickness of transition metal oxide slab and the lithium slab space in the bulk, during delithiation, as being reported from an XRD analysis. 14 The change of ϒ stress with potentials is especially important to be noted, as most literature 17,21,34 assumes constant solid mechanics properties when modeling particle-level stresses. For example, stress evolution depending on the lithiation rate 44 and a dramatic change of Young's modulus in layered structures with the degree of lithiation has been reported. 45 From the expression for ϒ stress , one can presume an increase in the Young's modulus or/and partial molar volume of the guest lithium with potential increase. From Figure 8, one can conclude that the thermodynamic and mechanical stress effects are more SOC-dependent than temperature-dependent.  Figure 9 shows the temporal variation of the fitted electrode potential response E s , the potential changes due to diffusion and stress, E conc and E stress , and reaction overpotential η, at ∼3.73 V, ∼3.83 V and ∼3.91 V for 10°C, 25°C and 40°C. E conc minus E stress equals the total change in the equilibrium potential E eq (c s ,σ) from which the system is perturbed during the GITT pulse. For different electrode voltages or temperatures, the performance limiting mechanism varies. In the case of low electrode voltage and low temperature, there is no to little stress contribution; the large overpotential is mainly dominated by the equal contribution from the activation overpotential and diffusion. As the temperature or the voltage increases, diffusion limitation and a prominent contribution from the particle stress term can be seen. At higher electrode potentials, the lower diffusivity and the higher SOC-dependent thermodynamic or stress factor (dE conc /dx pos or ϒ stress ) are responsible for the lithium transport limitation, while the activation overpotential decreases due to better reaction kinetics (higher k). The contribution from stress E stress is slightly lower than the pure diffusive process E conc . However, both E conc and E stress contribute to the polarization and slow potential relaxation seen during GITT. As temperature increases, the overpotential becomes more dominated by diffusion and stress. Better kinetics at higher temperature reduces the activation overpotential. Slightly poorer diffusion at higher temperature increases E conc . However, the stress term E stress that increases with temperature drives diffusion, contributing to the overall lower overpotential during polarization. During the relaxation process, due to slow diffusion, the buildup of large concentration gradients in both E conc and E stress requires more time to subside. In a nutshell, as electrode voltage and temperature increases, electrode performance is mainly limited by slow lithium transport. The effect of SOC on the lithium transport limitation is much more apparent than the temperature effect.
The combination of diffusivity, kinetics, thermodynamic and stress factor has an impact on the lithium concentration distribution that af-fects the speed of the relaxation process during GITT. As a result of slower interfacial kinetics and non-limiting diffusive behavior at lower potentials, lithium distribution in the spherical particle is smoother, leading to faster relaxation during open circuit. As electrode voltage increases, steeper concentration gradient arises due to the increasing limitation due to diffusion and stress. This leads to the aforementioned longer time constant for the relaxation process. To further illustrate the effect of the physical parameters on the different distribution of lithium, the spatial lithium concentration within the NMC particle is compared in Figure 10 at ∼3.73 V and ∼3.91 V, each at 25°C and 40°C respectively. The particle radius is normalized and the lithium concentration is normalized to the corresponding initial concentration. During charging (de-intercalation) in the first ten minutes, the lithium distribution gradually builds up, forming a gradient which drives the lithium transport along the particle radius from the center to the surface. The lithium concentration at the surface corresponds to the potential response of the electrode obtained from the GITT. Figures 10a  and 10c show that the lithium concentration distribution at ∼3.73 V is more homogenous with smaller gradient near the surface as a result of the good diffusive behavior. At ∼3.9 V, the diffusion occurs only in the particle's near-surface region during the 10-min charging as a results of slow diffusivity and increased stress properties (higher ϒ stress ), leading to steeper concentration gradients build up near the surface, as seen in Figures 10b and 10d. Consequently, the removal of the steep concentration gradients during GITT relaxation after which the current is interrupted takes a longer time.

Conclusions
A pseudo-two-dimensional porous electrode model incorporating a particle surface stress proportionality factor is developed. The equilibrium potential is based on linearization of the potential around the initial lithium concentration in a stress-and gradient-free state and the chemical potential change due to stress experienced by the solid lattice. The stress affects not only the equilibrium potential, but it also drives the lithium diffusion in the solid bulk. As the Newman-based standard model fails to describe sufficiently the dynamic performance of NMC111, the extended model with stress effect proves to be an improved version in successfully fitting the processes of polarization to the end of relaxation during a whole GITT transient. The GITT potential response curves are fitted well for NMC111 over a voltage range of 3.7-4.1 V vs Li/Li + at 10°C, 25°C and 40°C. Across the voltage range, D s,eff varies within an order of magnitude at around 10 −14 -10 −13 cm 2 /s. The strong dependency of the extracted k, dE conc /dx pos and ϒ stress with lithium fraction suggests the necessity of using varying kinetic, thermodynamic and electrochemo-mechanical properties to study intercalation materials. The fitting procedure shows minimal variance of D s,eff and dE conc /dx pos with temperature. From the model, low stress contribution is observed at low electrode voltage and low temperature; the overpotential dominated by the equal contribution from the activation overpotential and diffusion. As the temperature or/and the voltage increases, the limitation from the diffusion and particle surface stress are profound. The lithium transport behavior in solid NMC bulk as a result of the combined physical intercalation properties is further analyzed and compared. At higher potentials, a large lithium concentration gradient requires long relaxation time during GITT; at lower potentials, the lithium distribution is more homogenous leading to faster relaxation. The slow relaxation phenomenon is more prominent at higher temperatures. Therefore, this gives the cautious note that longer time could be needed to reach the equilibrium state at higher potentials and higher temperatures, which is particularly useful when doing measurements for the OCP curve. The generality of the derived stress proportionality factor in this model offers the possibility to be applied in other variants of NMC or insertion electrodes. It also serves as a groundwork for solid mechanics' analysis of stress during battery cycling to study electrode durability.