Physicochemically-informed continuum level model of a solid electrolyte interphase growth in Li-ion batteries

Contents lists available at ScienceDirect Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour Physicochemically-informed continuum level model of a solid electrolyte interphase growth in Li-ion batteries Klemen Zelič a,b, Meysam Esmaeilpour c, Saibal Jana c, Igor Mele a,b, Wolfgang Wenzel c, Tomaž Katrašnik a,b,∗ aUniversity of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, Ljubljana, Slovenia bNational Institute of Chemistry, Hajdrihova 19, Ljubljana, Slovenia cKarlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, Eggenstein-Leopoldshafen, Germany HIGHLIGHTS •A novel continuum model for SEI growth in Li-ion batteries, derived from a KMC model. •Innovative upscaling method preserves physicochemical consistency with a KMC model. •The model accurately predicts SEI properties such as thickness and porosity. •Demonstrates good agreement with both KMC simulations and experimental data. •Computationally efficient, allowing integration into continuum level battery models. ARTICLE INFO Keywords: Li-ion batteries Solid electrolyte interphase Molecular dynamics Kinetic Monte Carlo Continuum model Scale-bridging ABSTRACT Despite extensive research, understanding the SEI’s formation mechanism, structure, and its impact on battery performance remains challenging due to its complexity. To enable model-based design studies and to enhance understanding and prediction of the macroscopically observable consequences of SEI layer on battery performance and safety, continuum models featuring high level of prediction capability are needed. This objective of this paper is to resolve this challenge through an innovative physicochemically-informed continuum level model derived using a scale-bridging methodology, which, for the first time, enables highly consistent transfer of detailed KMC level based governing equations and reactions rates to the physicochemically-informed continuum level model. This was made possible by the innovative methodology relying on identification of rate-limiting reactions, deriving dynamic equations, and implementing dimensionality reduction. The resulting continuum model accurately replicates KMC results and experimental results while significantly reducing computational complexity. Furthermore, it, for the first time, enables distinguishing between ‘bad’, ‘good’, and ‘inorganic’ SEI growth scenarios on the continuum scale, offering valuable insights into electrode/electrolyte interface design. Due to its computational efficiency and scalability the proposed model can be integrated into higher-scale battery models, making possible advanced virtual performance, degradation and safety assessments with higher level of prediction capability. 1. Introduction The performance and safety of lithium-ion batteries (LIBs) are closely related to the formation and stability of the solid electrolyte interphase (SEI) on the surface of the anode active material [1–3]. The importance of this film, formed due to the reduction reaction of the electrolyte on the anode surface, was recognized many decades ago [4]. Despite extensive research efforts, the formation mechanism, structure, and impact of SEI on battery performance are still not fully understood ∗Corresponding author at: University of Ljubljana, Faculty of Mechanical Engineering, Aškerčeva 6, Ljubljana, Slovenia. E-mail address: [email protected] (T. Katrašnik). due to its high complexity [5–9]. Experimentally, the heterogeneity and evolving nature of the SEI make it difficult to isolate and study specific mechanisms without interference from other concurrent processes. This is compounded by limitations in current analytical techniques that may not fully capture the fast dynamics or the chemical diversity of the SEI layer. Theoretically, modeling these interactions requires assumptions and simplifications that can affect the accuracy of predictions. The lack of comprehensive, high-fidelity experimental data to validate these https://doi.org/10.1016/j.jpowsour.2024.235814 Received 12 July 2024; Received in revised form 17 October 2024; Accepted 6 November 2024 Journal of Power Sources 627 (2025) 235814 Available online 19 November 2024 0378-7753/© 2024 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/bync-nd/4.0/ ). K. Zelič et al. models further complicates theoretical approaches. A detailed discussion of the limitations of experimental and modeling techniques for SEI investigation is described in the Ref. [3]. These challenges significantly influence performance and technological application of batteries being one of the most important electrochemical energy storage technologies. In addition to extensive experimental investigations [5,7,10], an important part of the current knowledge on the formation mechanism of SEI has been obtained through theoretical studies based on modeling and simulation [11–14]. In particular, for the standard carbonatebased electrolytes, many studies based on atomistic simulations report underlying phenomena of SEI formation [15–18]. Recently, the coauthors of this work published an atomistically informed kinetic Monte Carlo (KMC) model [19], which provides interesting insight into key aspects of SEI characteristics, including thickness, porosity, and volume fraction. Through the analysis of these observables, the authors have categorized the final SEI thickness into three distinct groups based on their parameter study: ‘inorganic SEI’, ‘good SEI’ and ‘bad SEI’. According to [19], a ‘bad SEI’ forms as a very thin, porous, and frequently discontinuous organic layer, whereas a ‘good SEI’ develops into a thick, continuous organic layer. The term ‘inorganic SEI’ describes the formation of solely a compact inorganic layer. The descriptors ‘good SEI’ and ‘bad SEI’ are intentionally used to reflect their respective impacts on battery performance and longevity — positively in the case of ‘good SEI’ and negatively for ‘bad SEI’. The KMC model presented in [19] thus serves as an invaluable investigative tool for exploring a range of initial conditions of the electrolyte and electrode that determine the SEI layer type which will develop under specific starting conditions. Although, atomistic simulations are indispensable for revealing basic SEI formation phenomena, to study, understand and predict the macroscopically observable consequences of SEI layer on battery performance and safety, a continuum models are needed, e.g. [1,20]. Very good explanation on how continuum models can improve predictions compared to existing atomistic models is given in Refs. [21,22]. Continuum level models can feature different levels of detail and modeling basis, however, continuum SEI models are usually based on phenomenological approaches fitted to the electrochemical experiments (i.e., charge–discharge curves, electrical impedance spectroscopy (EIS), galvanostatic intermittent titration technique (GITT)) [12,23,24]. This is a shortcoming, as such types of models are lacking of prediction capability and they, in general, do not properly respond to variations in applied materials and their compositions, which is a serious deficiency for supporting model based design studies of batteries. In the literature, there are few reports of continuum models that are inspired by lower scale ansatzes, which are tested to obtain good agreement with experimental results [6,9]. Hence, modeling ansatzes were insightfully postulated by lower scale phenomena, while the complete model was not derived via mechanistic scale-bridging from lower scales. Several approaches can be used to address the scale bridging problems focusing on transferring knowledge from lower scale models to higher (continuum) scale models, particularly in the complex context of the Solid Electrolyte Interphase (SEI) layer formation in lithiumion batteries. The direct coupling of lower scale models to continuum models is currently hindered by the high complexity of SEI formation mechanisms and the multitude of chemical reactions involved, resulting in prohibitively high computational complexity and time. To tackle this challenge, one viable approach involves employing machine learning as an intermediate step between the lower and upper scale [25]. While this technique ensures low computational times, it does not offer an cause-and-effect insight into physical phenomena occurring at lower scales, hence, hindering obtaining detailed insight into potential mitigation strategies. Consequently, the development of a continuum level model, which is derived using a physicochemically consistent scale bridging from lower scales represents a challenge in maximizing trade-off between preservation of the mechanistic background of lower scale models and reaching computational complexity and times, which enable simulating continuum level models over longer spans of the operating time of the battery. The state of the art in this field is described in detail in the Refs. [3,26]. Despite recent advances, it can be succinctly concluded that this critical gap in achieving a physicochemically consistent scale bridging between discrete models, such as atomistic scale models, and continuum models still persist. To resolve this challenge, we present an innovative scale-bridging methodology, which, for the first time, enables highly consistent transfer of detailed KMC level based governing equations and reactions rates to the continuum scale, with only a minor loss in physicochemical consistency but with a very significant increase in computational efficiency. The proposed methodology consists of three steps: •Thermodynamic derivation of the dynamic equations for the system, •Identification and isolation of the rate limiting reactions and the predominant species in the system, •Dimensionality reduction. Resulting innovative continuum scale model, which preserves consistency to lower scales, was first calibrated to the KMC model and later bench-marked to the other sets of results generated by the KMC model to first demonstrate its accuracy and, in particular, prediction capability. In addition to the excellent agreement between the proposed model and the KMC model used for benchmarking, the simulation results also show strong alignment with published state-of-the-art experimental data. However, an even more important merit of the model arises from its physicochemical consistency with the KMC model, which is reflected in its unprecedented prediction capability of simulating formation of ‘inorganic SEI’, ‘good SEI’ and ‘bad SEI’, during calendar aging of cells, which was, to the best of authors knowledge, not yet demonstrated on the continuum level in a predictive manner. This characterizes the developed model as a suitable model for supporting early stage virtual design studies of battery performance and safety. 2. Methods Applied methodological steps are presented in Fig. 1. In this figure analytic derivation, which is one of the pillars of the scale bridging methodology, comprises: thermodynamic derivation of the dynamic equations for the system, identification and isolation of the rate limiting reactions and the predominant species in the system and dimensionality reduction. This pillar yields differential equations for species transport and reactions. This is a key pillar of the innovative scale bridging methodology, which ensures physiochemical consistency with atomistic models. In addition, calibration of the continuum scale model, which is the second pillar of the scale bridging methodology, comprises the optimization of model parameters with the aim to achieve the original KMC model by the use of the differential evolution algorithm. Both constituent parts of the scale bridging methodology (gray arrows on Fig. 1) result in the versatile continuum model which was solved by finite volume method to obtain the continuum model simulation results (rightmost cell of 1). These results are subsequently compared to the calibration and validation set of data generated by original KMC model in order to demonstrate the accuracy and prediction capability of the derived continuum model. Detailed description of each of these parts is presented in following subsections. 2.1. Derivation of continuum model equations Four distinct physico-chemical processes inherent to the system were analytically upscaled from the discrete KMC notation to continuum equations. These processes encompassed the migration of diffusing species, chemical reactions, surface adsorption and nucleation. The following subsections describe the derivation of continuum equations from KMC model for all listed phenomena, which ensures the physical consistency of the presented model with the lower scale. Journal of Power Sources 627 (2025) 235814 2 K. Zelič et al. Fig. 1. Schematic representation illustrating the scale - bridging methodology employed for the transformation of the original KMC model [19] to the continuum scale. Fig. 2. Schematic representation of the gas lattice model. Square mesh cells depict available empty spaces for diffusing species. Yellow and purple spheres symbolize individual particles of two diffusing species. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 2.1.1. Migration of diffusing species In the original KMC model [19], the behavior of diffusing species is characterized by what is known as a gas lattice model [27]. This model simulates the random motion of individual particles within the ensemble, where particles undergo random hops from their current lattice site to neighboring sites. The activation energy, dictating the height of the potential barrier between two lattice sites, determines the probability of a particle’s hopping event at each time step of the simulation. The system’s collective dynamics (i.e. the diffusion) is obtained by tracking the trajectories of each particle on the lattice. Fig. 2provides a schematic representation of two diffusing species within the context of the gas lattice model. To transition from the gas lattice model to the continuum model, we applied the regular solution theory [28–30]. Ultimate objective of this scale bridging step was to describe transport phenomena within the domain depicted in Fig. 2using two molar concentration fields instead of tracking the exact positions of each particle within the ensemble. According to Boltzmann’s definition, the total entropy of the system represented in Fig. 2is proportional to the logarithm of the number of all possible configurations at fixed particle number 𝑆= −𝑘𝐵log (𝑁 𝑁1𝑁2(𝑁−𝑁1−𝑁2)).(1) Here, 𝑆represents the system’s entropy, 𝑘𝐵is the Boltzmann constant, 𝑁denotes the total number of lattice sites, 𝑁1indicates the number of lattice sites occupied by the first species’ particles, and 𝑁2signifies the number of lattice sites occupied by the second species’ particles. In practical applications, 𝑁,𝑁1and 𝑁2assume large values, as each particle on the lattice represents a single molecule of the Li-ion battery electrolyte. Consequently, Stirling’s formula allows us to transform Eq. (1) into: 𝑆= −𝑘𝐵𝑁[log(𝑁) −𝑁1 𝑁log (𝑁1 𝑁)−𝑁2 𝑁log (𝑁2 𝑁) +𝑁−𝑁1−𝑁2 𝑁log (𝑁−𝑁1−𝑁2 𝑁)].(2) Note that the obtained expression is general and remains valid regardless of the model’s dimensionality. The same entropy expression would be derived whether starting from a 1D, 2D, or 3D lattice. We introduce non-dimensional molar concentrations, 𝑐1=𝑁1∕𝑁 and 𝑐2=𝑁2∕𝑁, for each species. This transformation yielded the following formula for the free energy density of the system: 𝑓= −𝑇 𝑠=𝑅𝑇 [𝑐1log(𝑐1) −𝑐2log(𝑐2) − (1 −𝑐1−𝑐2) log (1 −𝑐1−𝑐2)], (3) where 𝑅represents the gas constant, 𝑇stands for temperature, and 𝑠 denotes entropy 𝑆per mole of lattice sites. The volume integration of the derived free energy density in Eq. (3) yielded the total energy density of the system, F. Minimizing F represents the driving force propelling the system toward equilibrium. The chemical potential 𝜇of the system, which governs the diffusion process, was computed from F using the Lagrange variations principle [31]. The resulting chemical potential dependency on concentration is described by: 𝜇=𝑅𝑇 [log 𝑐1+ log 𝑐2− log(1 −𝑐1−𝑐2)].(4) This obtained chemical potential conforms to the three-component regular solution [32]. By introducing the derived 𝜇values into the Cahn–Hilliard equation [33], we established a system of two partial differential equations that describe the transport of two diffusing/interacting species within: 𝜕 𝑐𝑖 𝜕 𝑡= ∇(𝑀𝑖𝑐𝑖∇𝜇)+ ∇𝑗 .(5) Journal of Power Sources 627 (2025) 235814 3 K. Zelič et al. Here, 𝑖ranges over one and two for the first and second species, 𝑡 denotes time, 𝑀𝑖represents the mobility of the 𝑖th species and ∇𝑗term represent the divergence of any additional molar fluxes or sources in the system. The final forms of these equations, explicitly incorporating 𝜇, will be presented in Section 3. 2.1.2. Description of reactions In the original KMC model, reaction rates were described using activation energies that determine the energy barriers and the probability of two colliding particles (reactants) to undergo a reaction. The probability of particle collision was inherently present in the gas lattice (KMC) model, as every particle trajectory was tracked throughout the simulation. To transform this discrete reaction model into the continuum model, we needed to describe two critical aspects that determined the reaction rate: the probability of collision and the probability of reaction. This translated into the commonly used proportionality of the reaction rate to the reaction rate constant (𝑘, representing the probability of reaction) and the reactants’ concentration (𝑐𝑖, representing the probability of collision). 𝑟=𝑘𝑃𝑡∏ 𝑖 𝑐𝑛𝑖 𝑖(6) Additional factor denoted 𝑃𝑡= exp(−𝑑∕𝜆)was added in all reactions that explicitly include electrons as a reacting species. 𝑑denotes distance between surface and site of the reaction and 𝜆is characteristic length. 𝑃𝑡that takes in to account probability for electron tunneling to the site of the reaction [34]. KMC model use identical description of tunneling probability. Since the probability of collisions in the KMC (gas lattice) model was linearly proportional to the concentration of particles, all reactions were modeled as zero-order reactions. Consequently, the reaction rate is linearly proportional to the concentration of reactants and inversely proportional to the concentration of products. Eq. (6), coupled with the fact that all reactions in the system were zero-order, led to the formulation of a reaction flux 𝑗 𝑟= ∇⋅𝑗 ,(7) which could be seamlessly integrated with the previously obtained Cahn–Hilliard equation (Eq. (5)). 2.1.3. Surface adsorption In continuum model, SEI film growth is modeled by the reaction equation at the surface, which transforms diffusing species in to the organic SEI. Reaction equation in most basic form reads: 𝑟=𝑘𝐴𝑐 .(8) Here 𝑟denotes reaction rate, 𝑘is reaction rate constant, 𝐴is surface area where reaction takes place and 𝑐is concentration of reactant species. Note that the surface where reaction happens is reactive and it enters in to the reaction as reactant. 𝐴as a measure of quantity of such a surface, which can in this case be interpreted as amount of this reactant. In a porous SEI film, this is a total surface area of porous media. Since porosity and the amount of SEI film vary with the distance from graphite surface, it is convenient to introduce the space dependent volumetric density of surface area: 𝑎=𝑑 𝐴 𝑑 𝑉,(9) where 𝑑 𝑉is infinitesimal volume around the point in space where 𝑎is defined. This newly introduced parameter can further be linked to the porosity, which will be a dependent variable in the final set of equation and one of the main results of simulation: 1 −𝜀=𝜁 𝑎. (10) Here 𝜁is a new parameter in the proposed model, which corresponds to the size of average feature in the porous structure of SEI. Incorporating previous equations in Eq. (8) yields: 𝑟=∫𝑘 𝜁(1 −𝜀)𝑐 𝑑 𝑉 .(11) Table 1 List of the reactions and species involved in SEI growth according to Ref. [19]. Entry Reactants Products 1 Electrode (surf) + EC-Li+C2H4OCOOLi 2 Electrode (surf) + C2H4OCOOLi + Li+Li2CO3+ C2H4 3 Li2CO3(surf) + C2H4OCOOLi + Li+Li2CO3+ C2H4 4 Li2CO3+ EC-Li+C2H4OCOOLi 5 C2H4OCOOLi + C2H4OCOOLi Li2EDC + C2H4 6 Li2EDC + Li2EDC (Li2EDC)2 7 (Li2EDC)2+ Li2EDC SEI cluster 8 Li2EDC + SEI cluster SEI cluster 9 (Li2EDC)2+ SEI cluster SEI cluster 10 (Li2EDC)2+ (Li2EDC)2SEI cluster 11 SEI cluster + SEI cluster SEI cluster In this form, surface reaction equation is prepared to be coupled to Eq. (5) in a form of a source term in accordance with Eq. (7). 2.1.4. Nucleation of porous SEI Nucleation of porous SEI was modeled based on the probability of nucleation driven by nucleation activation energy and concentration of reactants (mobile precursors). Probability for the nucleation was modeled by the equation 𝑃𝑛𝑢𝑐 =𝜅 𝑐 .(12) Here, 𝑃𝑛𝑢𝑐 represents the spatially dependent probability for the nucleation of porous SEI, while 𝜅signifies the rate of nucleation. The rate constant 𝜅is derived directly from the nucleation energy, as governed by the Arrhenius Equation [35], and 𝑐denotes the local concentration of SEI precursors. 2.2. Isolation of the limiting reactions and the predominant species In the original kinetic Monte Carlo (KMC) model [19], eleven major reaction steps governing SEI formation were identified. These reactions are listed in Table 1. EC represents ethylene carbonate molecule, Li+is lithium ion, Li2EDC denotes dilithium ethylene decarbonate and SEI cluster is a cluster of several (more than two) dilithium ethylene decarbonate molecules [19]. From the comprehensive list of reactions and reacting species in Table 1, we discerned the most significant (limiting) reactions and species. This determination relied on the sensitivity analysis of simulation results obtained from the original KMC model, which provided insights into the temporal evolution of the system as a function of reaction rates. To maximize the trade-off between preserving the mechanistic background of lower-scale models and achieving adequate computational complexity in the derived modeling approach, we considered the most significant species (C2H4OCOOLi — later referred to as SEI precursor, Li2CO3, Li2EDC, and SEI cluster), and reactions, which are implemented at the continuum level and are summarized in Table 2. Isolating significant reactions must be approached with caution, as while it can significantly improve the model’s computational efficiency, excessive lumping may compromise both the accuracy and physicochemical consistency of the model section. Systematic analysis of all reaction rates for the reaction from Table 2was performed ba original KMC model. The impact of each parameter in this sensitivity analysis was evaluated based on the final SEI observables — SEI thickness, porosity, and volume fraction (described later in 3.5), which are detailed in Section 3.5 of the manuscript. The most influential parameters were identified as those whose perturbation resulted in the largest average deviation across these three SEI observables in the KMC simulation results. During the isolation of significant reactions, we considered each single reaction step in the continuum model that corresponds to the Journal of Power Sources 627 (2025) 235814 4 K. Zelič et al. Table 2 List of the species and reactions used in the presented continuum model. Entry Reactants Products 1 Electrode (surf) + EC-Li+SEI precursors 2 Electrode (surf) + SEI precursors Li2CO3+ C2H4 3 Li2CO3(surf) + EC-Li+SEI precursors 4 SEI precursors + SEI precursors Li2EDC + C2H4 5 Li2EDC + Li2EDC SEI cluster formation of the same species depicted through multiple possible events in the original KMC model. The only exception was for the formation of SEI precursors, which primarily occur as intermediates on two distinct surfaces — graphite and Li2CO3. Therefore, reactions 1 and 4 from Table 1in the manuscript correspond to reactions 1 and 3 in Table 2, respectively. In Table 1, reactions 2 and 3, which detail the formation of the compact SEI component Li2CO3, are consolidated into a single reaction (reaction 2) in Table 2. Reaction 5 in Table 1, representing the main component of porous SEI, is aligned with reaction 4 in Table 2. Additionally, reactions 6–11 in Table 1, which signify the aggregation and clustering processes for forming porous SEI, are collectively represented by reaction 5 in Table 2. The color coding in Tables 1and 2highlights the correspondence between reactions across both tables. This consolidation reduces the total number of reactions from eleven to five, streamlining the presentation based on the chemistry of the reaction processes. To establish the final set of equations describing the system, Table 2 must be linked to the approaches detailed in Section 2.1 All mobile species in the system (EC molecules, SEI precursors, and Li2EDC) are modeled using the equation derived in Section 2.1.1. The diffusive concentration fields determine the concentration-dependent reaction rates for all species, as outlined in Table 2. Reactions 1, 3, and 4 from Table 2are modeled using the methods described in Section 2.1.2. Reaction 2 represents the adsorption of newly formed Li2CO3onto the graphite surface and is therefore modeled using the approach outlined in Section 2.1.3. Reaction 5 involves the nucleation of porous SEI and subsequent adsorption to the nucleated clusters, and is consequently modeled by the methods in Sections 2.1.3 and 2.1.4. All reaction rates that correspond to the reactions listed in Table 2, were determined by the calibrating the model to the KMC simulation results as later described in Section 3.3. 2.3. Dimensionality reduction The original KMC model was implemented on a two-dimensional lattice, a common approach for describing diffusive media using a gas lattice model. To further reduce computational complexity, the derived continuum model was designed to deliver high fidelity results in one dimension. This dimensional reduction did not significantly compromise accuracy due to the system’s distinct geometry. The system namely exhibits rotational symmetry around an axis perpendicular to the graphite surface, and most dynamic processes occur along this axis. 3. Integrated continuum model 3.1. Bulk equations Integration of Eqs. (4)–(11) provides the final set of differential equations, describing time evolution of SEI growth on the graphite surface in the surrounding electrolyte, which represents bulk of the calculation domain: 𝜕(𝜀1+𝜀2)𝑐1 𝜕 𝑡=𝐷1 𝜀𝛽∇(∇𝑐1)+𝐷1 𝜀𝛽∇(1 1 −𝑐1−𝑐2 ∇𝑐1) −𝑘1𝑃𝑡 𝑐1 𝑐2 −𝑘4∫(1 −𝜀2)𝑐1𝑑 𝑉 ,(13) 𝜕(𝜀1+𝜀2)𝑐2 𝜕 𝑡=𝐷2 𝜀𝛽∇(∇𝑐2)+𝐷2 𝜀𝛽∇(1 1 −𝑐1−𝑐2 ∇𝑐2)+𝑘1𝑃𝑡 𝑐1 𝑐2 −𝑘2∫(1 −𝜀1)𝑐2𝑑 𝑉−𝑘3∫(1 −𝜀2)𝑐2𝑑 𝑉 ,(14) 𝜕 𝜀1 𝜕 𝑡= −𝑘2∫(1 −𝜀1)𝑐2𝑑 𝑉 ,(15) 𝜕 𝜀2 𝜕 𝑡= −𝑘3∫(1 −𝜀2)𝑐2𝑑 𝑉−𝑘4∫(1 −𝜀2)𝑐1𝑑 𝑉 .(16) Effective diffusivities were prescribed to the porous media with the Bruggeman coefficient 𝛽[36]. Here, 𝜀1and 𝜀2denote the porosities of the porous and compact SEI layers, respectively. A value of one indicates the absence of the SEI layer at that coordinate. Given that porosity is defined as the fraction of volume unoccupied, it is effectively treated as a non-dimensional concentration. Additionally, 𝑐1and 𝑐2represent the concentrations of EC molecules and SEI precursors, while 𝐷1and 𝐷2stand for the diffusivity coefficients of EC molecules and SEI precursors within the surrounding environment. The reaction coefficients are denoted by 𝑘, where 𝑘1,𝑘2,𝑘3, and 𝑘4correspond to specific reactions as outlined in Table 2, with 𝑘4representing a joint reaction rate for reactions 4 and 5. Consistent with the assumptions of the original KMC model, it is postulated here that the migration coefficient for Li+ions is significantly higher compared to that of EC molecules and SEI precursors. Consequently, the diffusion of Li+ions is not explicitly modeled, as their rapid migration enables omitting their detailed tracking at the envisaged timescales and accuracy of the simulation. 3.2. Boundary conditions The exact same boundary conditions as in KMC [19] were used in continuum model: ∇𝑐1|𝑧=0 = 0,∇𝑐2|𝑧=0 = 0,∇𝜀1|𝑧=0 = 0,∇𝜀2|𝑧=0 = 0, ∇𝑐1|𝑧=0 =𝑐0 1,∇𝑐2|𝑧=∞ = 0,∇𝜀1|𝑧=∞ = 0,∇𝜀2|𝑧=∞ = 0, where 𝑐0 1represents initial (bulk) concentration of lithium ions in electrolyte. 3.3. Model calibration The first step in scale bridging methodology presented in Section 2.2 Isolation of the limiting reactions and the predominant species, ensured low computational complexity of the continuum model, whereas activation energies of the KMC model need to be adapted to enable their use as input parameter for reaction rates and diffusion coefficients in the continuum model. The data driven approach was used to determine those parameters from the results of the KMC model. Differential evolution algorithm [37] was applied to calibrate parameters of the continuum model. For the calibration with the use of differential evolution seventy different KMC generated time traces were used. For each of seventy used time traces, calibration was done for six distinctive time steps including the last one, representing the equilibrium state at the end of simulation. The cumulative square root of a differences between continuum model result and KMC model result for each point in space and time was calculated and use as a cost function value in optimization. Journal of Power Sources 627 (2025) 235814 5 K. Zelič et al. Fig. 3. Schematic representation illustrating the model domain and its discretization. 3.4. Solving modeling equations The system of equations was discretized using the finite volume method in one dimension, specifically in the direction normal to the graphite surface. The system was solved employing the implicit Powell’s hybrid iteration technique. Discretization of the domain was designed to account for all boundary conditions established at an infinite distance from the surface. This domain discretization is visually depicted in Fig. 3. Computational domain is comprised of two distinct parts: the KMC domain and the auxiliary domain. The KMC domain is of primary interest to us as it is the region where SEI formation occurs. This part of the domain is referred to as the KMC domain because it applies to the same region as the original KMC model, maintaining the same resolution of 1 nm. The auxiliary domain, on the other hand, is introduced solely to satisfy the three boundary conditions defined at an infinite distance from the surface. Since the dynamics of the system become less significant when moving away from the surface, we have employed a resolution reduction scheme within the auxiliary domain to maintain computational efficiency of the continuum model. This resolution reduction is defined using an exponential function, expressed as 𝛥𝑧 =𝛥𝑧0exp (𝑧∕𝜁), where 𝜁is a numerical parameter, that was tuned to ensure the maximal boundary condition accuracy and computational speed. 3.5. SEI observables In line with the original KMC model from Ref. [19], three SEI observables were calculated based on concentration trends obtained from continuum model simulations. These SEI observables — thickness, porosity, and SEI fraction were computed in the same manner as in Ref. [19] to enable a direct comparison of results between the two models. The equations defining these observables are as follows: 𝑉𝑓=1 𝐿∫𝐿 0[(1 −𝜀1)+(1 −𝜀2)]𝑑 𝑧, (17) 𝜀 =𝜀1+𝜀2 𝑉𝑓 ,(18) 𝛿=∫𝐿 0[𝑑(𝜀1+𝜀2) 𝑑 𝑧𝑧]𝑑 𝑧, (19) where 𝑉𝑓denotes the volume fraction of the SEI, 𝜀 is the average porosity of the SEI, and 𝛿is the estimated SEI thickness. A detailed explanation of these equations can be found in the supplementary information of Ref. [19]. 4. Results To showcase the accuracy of the continuum model, we replicate the results obtained from the KMC model, as originally presented in the Ref. [19]. This serves as a validation case. In Fig. 4, a comparison is made between the results of the original KMC model and the presented continuum model at six different time points. Fig. 4demonstrates excellent agreement between the results of the KMC simulation and the continuum model simulation. A minor discrepancy between the simulation results is evident in the case of SEI precursors. This was anticipated, as reactions 2 and 3 from Table 1 were consolidated into a single reaction 2 in Table 2. Nevertheless, the obtained approximate agreement in this trend supports and justifies the feasibility of combining these two reactions into one. The observed agreement justifies the scale bridging steps detailed in Section 2, which were implemented on the original KMC model to reduce computational complexity. In the presented case, the computational time of the continuum model simulation was 30 times shorter compared to the computation time of the presented KMC simulation. Fig. 4affirms that developed continuum model aptly captures the spatial dependencies of all predominant species in the system. To assess temporal dependencies, Fig. 5compares the evolution of species in both the KMC and continuum models. The plotted values in Fig. 5represent the temporal dependence of the overall quantity of each species within the computational domain. These values were obtained by integrating spatially resolved concentrations across the entire domain for each simulation time step. In Fig. 5, the green line in the continuum model result plot represents the sum of the orange, green, and blue lines from the KMC result. This aggregation is necessary as the continuum model describes all SEI precursors as a single species (Table 2). The alignment between trends in the KMC and continuum model results indicates that characteristic times for reactions and species transport in the system are consistent between both models. To assess the predictive capabilities of the developed model, we conducted a parametric study, systematically varying the model parameters. Subsequently, we compared the obtained results with those derived from the original Kinetic Monte Carlo (KMC). The noteworthy consistency observed between the outcomes of our model and the original KMC, particularly in cases not included in the optimization process (validation cases), underscores the elevated predictive accuracy of our continuum model. The parametric study was conducted through an examination of parameter sensitivity in the presented model. Our objective was to evaluate the impact of individual reaction rates and diffusion coefficients on the final SEI thickness — a crucial parameter in the system, determining the quality of the SEI layer that culminates in either a ‘good’ or ‘bad’ SEI. This analysis provides valuable insights for designing an optimal electrode surface and electrolyte to achieve the desired performance. Sensitivity analysis of the SEI thickness on model parameters was performed using Fisher information matrix (FIM). Due to the stochastic nature of porous SEI nucleation event in the continuum model, simulations were run 3072 times to obtain the average values of SEI thickness and its derivatives, needed for FIM calculation. Fig. 6illustrates the obtained FIM. Analysis of the matrix reveals that the most influential parameter on the final SEI thickness is the reaction rate between EC molecules and the graphite surface (denoted by 𝑘1in Eqs. (13) and (14)), forming SEI precursors. The diffusion coefficient of the SEI precursor species emerges as the second most influential parameter. These two parameters exhibit a slight correlation, which is expected given their association with the same species. All other parameters do not significantly influence final SEI thickness and they do also not show any significant correlation between each other, which once again justifies adequacy of the isolation of the limiting reactions and the predominant species in the scale bridging methodology. Utilizing information derived from the FIM matrix, we conducted a systematic analysis to investigate the impact of the reaction rate between EC molecules and the graphite surface on the model outcomes, recognizing this as the most influential parameter. The results revealed a significant influence of this parameter on the quality of the SEI formed during the cycling process. The classification proposed by [19] divides the SEI layer into three categories: ‘good SEI’, ‘bad SEI’, and ‘inorganic Journal of Power Sources 627 (2025) 235814 6 K. Zelič et al. Fig. 4. Comparison between simulation results of the KMC and continuum models. Solid lines depict the results of the continuum model simulations, while dots represent the results of the KMC model simulations. Fig. 5. Comparison of time dependencies of cumulative species amounts in the domain between (a) KMC (Reproduced with permission from [19]. ©2023 Wiley-VCH GmbH. Published by Wiley-VCH.) and (b) continuum model simulation results. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Journal of Power Sources 627 (2025) 235814 7 K. Zelič et al. Fig. 6. Fischer information matrix corresponding to the case presented in Fig. 4. In this matrix, the values of the elements reflect the sensitivity of the model to specific parameters (diagonal elements) and the correlation between different model parameters (off-diagonal elements). These values are obtained by observing the model’s response to slight perturbations in the parameters [38]. Notably, the highest magnitude of matrix elements occurs for the parameter 𝑘1(yellow color). This suggests that the studied system is most sensitive to perturbations in the kinetics of reaction 1, as outlined in Table 2. Following this, the second most influential parameter is identified as 𝐷2, the diffusion coefficient of SEI precursors. Moreover, the absence of any remarkably high correlations between parameters justifies the scale-bridging steps described in Section 2.2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) SEI’. A ‘good SEI’ is characterized by high thickness and low porosity, whereas a ’bad SEI’ exhibits the opposite characteristics. The term ‘inorganic SEI’ refers to conditions where only compact inorganic SEI is formed. Based on the results of the KMC simulations, Ref. [19] shows the phase diagram that proves the existence of three regions of parameter combinations in which ‘‘good SEI’’, ‘‘bad SEI’’ and ‘‘compact SEI’’ grow. We show here that parameter 𝑘1is the most influential parameter of the model, as it makes possible shifting the system across all there regions of such a phase diagram. It can, therefore, serve as a single parameter influencing formation of good and bad SEI. Fig. 7illustrates the final equilibrium states of the simulation for various values of 𝑘1, with all other parameters held constant as in the simulation presented in Fig. 4. Analysis of the 𝑘1parameter, as illustrated in Fig. 7shows that a systematic reduction in the reaction rate (𝑘1) led to the formation of a thinner Solid SEI, eventually transitioning from a ‘good SEI’ to a ‘bad SEI’. These observed trends were then compared to the ‘bad SEI’ example obtained from the KMC model in Ref. [19]. Solely by manipulating the reaction rate 𝑘1, we not only facilitated the transition between regions in the bad/good SEI phase diagram but also maintained a high degree of alignment between the results of the continuum model and the SEI model. The comparison of both models for the case of ‘bad SEI’ is depicted in Fig. 8. The capability of our model to identify the primary influence of the reaction rate of EC molecules on the graphite surface and the diffusion of SEI precursors in the electrolyte underscores its efficacy as a powerful digital platform for researching and designing improved electrode/electrolyte interfaces in respect to SEI observables of interest. The parametric study, presented in Figs. 7and 8, was also utilized for model validation. Given that experimental procedures to obtain insights into key observables (such as SEI porosity, thickness, and composition) essential for validation are extremely demanding, we rely on previously published experimental data. For this validation, data from Ref. [39] were selected, as it provides valuable insights into SEI evolution during calendar aging of Li-ion cells with a graphite anode. This study was chosen due to the similarity between the conditions in presented simulations and experimental setup used in Ref. [39], specifically the SEI formation during calendar aging on a non-lithiated graphite surface, without cycling. The data from [39], acquired using a combination of SEM-FIB and XPS techniques, were used for comparison. The final SEI thicknesses reported for electrodes with low SOC in [39] range from 8 to 9 nm (Table 2 in [39]), which aligns with the thickness range predicted by presented model for the same conditions (Fig. 8) predicting value of 9.81. Similar values for final SEI thickness are also reported in other experimental studies that estimated the SEI growth during cell cycling [5,40,41] in the range of 8 to 25 nm. Further validation of the model is demonstrated through the parametric study in Fig. 7. Regardless of the magnitude of the most influential parameter (𝑘1), presented simulations produced SEI thicknesses between 3.99 nm and 22.40 nm, which closely match the experimental values from Ref. [39]. The average thickness from all simulations in presented parametric study was 9 nm, exactly matching the reported value in [39]. Additionally, we compared the final SEI compositions from presented simulations with the experimental results presented in Figures 4 and 5 of [39]. The experimental data clearly indicate that SEI formed during calendar aging consists of a thin, compact layer near the graphite surface, with a thicker, porous layer extending into the electrolyte — this same pattern is observed in presented simulation results. Furthermore, the model presented offers potential for further development. Currently tailored to estimate SEI evolution during calendar aging, the next step involves implementing variable electrode potentials. This enhancement will significantly broaden research and design possibilities at the electrode/electrolyte interface. Additionally, it will provide further opportunities to validate the entire modeling framework and the corresponding upscaling methodology. 5. Conclusions This paper introduces a continuum model for SEI layer growth, which is derived through consistent upscaling from the KMC model and which features three significant contributions to the field. The first major contribution lies in the development of a scalebridging methodology, employed to derive the proposed model from the existing KMC model. This methodology involves three pivotal steps: (i) identification and isolation of limiting reactions and predominant species, (ii) derivation of dynamical equations for all system species, and (iii) implementation of a dimensionality reduction methodology. Demonstrating its efficacy, this scale-bridging approach significantly reduces computational times (by up to 30 times in the presented case), while concurrently preserving accuracy and physico-chemical relevance. Importantly, this methodology is versatile and can be applied to various molecular dynamics models, facilitating the creation of upscaled continuum models capable of describing phenomena with comparable accuracy and reduced computational complexity. The second key contribution is the development of an SEI growth model, which serves as a valuable tool for fundamental research and design exploration at the electrode/electrolyte interface. Notably, this is the first continuum model to differentiate between ‘‘bad’’, ‘‘good’’, and ‘‘inorganic’’ SEI growth scenarios, as defined in [19]. The model reliably predicts critical properties of the final SEI layer, such as thickness, porosity, and volume fraction, which were key observables used to demonstrate the strong agreement between the simulation results and published experimental data. In this paper, the model’s capabilities are demonstrated using a simple combination of an EC-based electrolyte and a pure graphite electrode. Due to its flexibility, the model can be easily parameterized for different chemistries, including alternative electrode materials or electrolytes, and can accommodate simulations involving mixtures of multiple solvents. The third and final contribution underscores the suitability of the presented model for integration into higher-scale models. With its low computational complexity, continuum nature, scalability potential, and high accuracy, the model is tailored for coupling with electrochemical Journal of Power Sources 627 (2025) 235814 8 K. Zelič et al. Fig. 7. Parameter Sensitivity Analysis: Dependency of the equilibrium state, calculated by the continuum model, on the reaction rate 𝑘1with respect to the quality of the SEI. Journal of Power Sources 627 (2025) 235814 9