# A microstructure-resolved model of sodium-ion battery hard carbon electrodes

Imelda Cardenas-Sierra, Francisco Fernandez, Martin Petit, Alejandro Franco

► To cite this version:

Imelda Cardenas-Sierra, Francisco Fernandez, Martin Petit, Alejandro Franco. A microstructure-resolved model of sodium-ion battery hard carbon electrodes. Journal of Power Sources, 2025, 655, pp.237878. 10.1016/j.jpowsour.2025.237878. hal-05214382

HAL Id: hal-05214382

https://u-picardie.hal.science/hal-05214382v1

Submitted on 27 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LIBs. The main reason for choosing HC was the discovery that this carbonaceous material could reach a specific capacity of 300 mAh g⁻¹ at an average potential of 0.3 V vs. Na⁺/Na [6,7]. This value contrasts with the one provided by graphite, which could only offer a fraction of this capacity when it was first tested as a negative electrode for SIBs [8]. Later, it was revealed that the low capacity was caused by a thermodynamically unfavorable formation of graphite-intercalation compounds [9–11]. Another reason for the interest in HC is its inherent environmental friendliness, as it can be readily synthesized from the pyrolysis of different organic precursors, such as glucose [8,12] and sucrose [13,14], to mention some of the most studied ones. More recently, biomass coming from different sources such as coconut sprouts [15], or tannins from the mimosa plant [16], has also been studied. However, synthesis of HC is not constrained to the aforementioned precursors, since materials like resins, and synthetic polymers, can also be used [17].

HC is a non-graphitizable carbon, and its structural characteristics depend on the precursor choice and synthesis conditions. But even if HC does not present a long-range structural order, it has small nanosized domains with local crystallinity, caused by the stacking of curved graphene sheets in a disordered fashion [18,19]. These microstructural characteristics ultimately affect its electrochemical properties, and depend heavily on the quality of the precursor, especially in the case where HC is derived from biomass. These variations in microstructure have led to some discrepancies in the electrochemical data analysis. Depending on the characterization technique, the probing time scale, and the bulk and surface sensitivities, different sodiation mechanisms can be identified [20].

When HC electrodes are subjected to a galvanostatic discharge, two main zones can be identified. First, a zone where the voltage reduces rapidly upon sodiation (going from around 1 V to 0.1 V vs. Na⁺/Na), referred to as the “sloping region”. Second, a zone where the voltage remains almost constant at very low potentials (from 0.1 V to 0 V vs. Na⁺/Na), referred to as the “plateau region” [21]. Despite the fact that these two zones are easily distinguishable, the scientific community has struggled to identify a universal mechanism that explains which process

produces each characteristic voltage behavior. Nevertheless, they have managed to identify the occurrence of the following storage processes in HC [18]:

i. Storage of Na ions on the particle's surface in contact with the electrolyte (which can be analogous to a capacitive process where an electrical double layer is formed);

ii. Storage of Na ions on the internal surface area of HC (this can be related to point defects, heteroatoms, or functional groups, and is commonly called “adsorption process”);

iii. Storage of Na between curved graphene layers in the turbostratic domains, which may also be affected by the presence of defects or heteroatoms (commonly known as “intercalation” process);

iv. Storage of Na clusters, located inside the particle nanopores that remain inaccessible to the electrolyte (commonly known as a “nanopore filling” process).

A schematic of all these processes is shown in Fig. 1a. However, even if these main storage processes have been identified, there is currently no agreement on attributing these processes to a specific region of the galvanostatic charge/discharge profile [20]. This debate is well reported in the literature, since many different research groups have ascribed the sloping and plateau regions to different sodiation processes. So far, the main sodiation mechanisms proposed [22] can be categorized in two types: (i) the “two-step” mechanisms [6,23–27], where one process is ascribed to each voltage region; and (ii) the “multi-stage” mechanism [13,28], where at least three sodiation processes are ascribed to the electrochemical footprint of HC, and one of the zones can be associated to multiple sodiation processes. As an example, the first proposed mechanism was published in the 2000s by Stevens and Dahn [6] which explained that the sodiation in HC occurred through an “Intercalation-Filling” process. In the cited work, the sloping region of the electrochemical discharge curve was associated to the intercalation of sodium in graphene layers, while the plateau region was attributed to the nanopore filling, as shown in Fig. 1b.

Given that the sodiation mechanism of the HC electrode is still not

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** a) Sodium storage phenomena in Hard Carbon (HC) electrodes, b) Electrochemical signature of HC electrodes upon the first discharge and charge at constant current, and c) Proposed “Intercalation-Filling” model for HC with 3D-resolved microstructures at the mesoscale.</figcaption></figure>
well established, it is crucial to employ alternative techniques to study this system, which will lead to key findings that improve its electrochemical performance. Computational simulations are one method that researchers have employed to gain insights into electrode performance. These simulations can be complementarily used with experimental techniques for the development of next-generation batteries, as they can help predict and understand different phenomena that cannot be easily examined by experiments [29]. A vast collection of physics-based models has been used for studying LIB electrodes behavior at the continuum scale [30]. In previous works of our research group, and in the context of the ARTISTIC Digitalization Initiative [29], mesoscopic 3D-resolved models were used to correlate the electrochemical performance behavior with the manufacturing parameters of the electrodes and cells [31–33,53]. These microstructure-resolved electrochemical modeling studies of LIB electrodes have pointed out the importance of microstructural heterogeneities (such as the interfaces between Active Material (AM), Carbon Binder Domain (CBD), and pores) in determining the practical electrode performance [34–38]. This approach proves to be particularly useful to study different features of an electrode’s microstructure, while providing insights into how manufacturing parameters affect the overall electrochemical behavior [39].

When it comes to SIBs, electrochemical modeling approaches are still in the early stages, with only a few studies having been published. For example, Chayambuka et al. proposed using the pseudo-two-dimensional (P2D) model in combination with experimental parameter determination for SIBs modeling [40,41]. This led to an experimentally validated full-cell electrochemical model (for Na₃V₂(PO₄)₂F₃ vs. HC) that was used to study galvanostatic discharge profiles at different C-rates. In addition, Garapati et al. followed the previous work by developing a reduced order model by considering a Single Particle Model (SPM) version of the same SIB full cell configuration [42]. They compared their results with experimental data from literature and simulation results obtained from the previous P2D model established by Chayambuka et al. [40,41]. Their findings showed that the SPM approach was less computationally expensive than the P2D one, albeit with lower accuracy at high C-rates. Following that, Neumann et al. developed a stochastic model for constructing virtual electrode microstructures of a Na₃V₂(PO₄)₃/carbon composite [43]. This was achieved using a Pluri-Gaussian technique to create a parametric stochastic model that used imaging data inputs to predict the composite electrode microstructures with varying contents of AM, carbon additive, and pore phases. However, this work focused on developing the positive electrode microstructures, and neither electrode microstructures nor electrochemical simulations were reported for the SIB negative electrodes.

It is worth underlining that a 3D-resolved electrode microstructure has not been considered in any of the earlier HC electrode representations or electrochemical models. Additionally, the complex sodiation mechanism was not investigated, as it was assumed that the sodium insertion in HC occurred by intercalation in spherical particles. Thus, a research gap still exists in developing an electrochemical model that enables the study of the heterogeneities in the electrode’s microstructure. In this article, we address these aspects for the first time by considering a 3D-resolved electrochemical model for HC, using Finite Element Method (FEM) computations (as schematized in Fig. 1c).

In this work, a half-cell model of a HC electrode is developed by using FEM. The 3D-resolved electrode microstructures are stochastically generated to represent the formulation and porosity found in an experimental study. Then, the structures are used inside the electrochemical model, one that is parameterized with experimental data from the same HC sample. The model is used to elucidate how changing the formulation and porosity of the electrode affects the discharge performance of the HC material at the electrode level. In addition, a particular interest is placed on exploring the sodiation processes of HC, in this case we model the Intercalation-Filling mechanism by taking inspiration from models describing blended electrodes. As for our knowledge, this is the first time that a 3D-resolved HC electrochemical model solved using

FEM has been proposed.

The next sections of the article are organized as follows. In Section 2, we present the experimental setup that inspired our model development. Then, in Section 3, we describe our modeling workflow, from the generation of stochastic electrode microstructures, to the potential calculation (adapted to study the two-step sodiation mechanism), and finally, the formulation of the electrochemical model. In Section 4, we discuss the parameterization and validation of the model by comparing simulation results to the experimental data, in addition, we show the effect of changing the electrode’s manufacturing parameters (formulation and porosity), on the electrochemical performance. Finally, we conclude and indicate perspectives for our work in Section 5.

## 2. Simulation setup

### 2.1. Cell-material specifications

A HC electrode previously synthesized in our laboratory from a mimosa tannin extract, was used for the parameterization and validation of our model. Specifically, the HC1600 sample (synthesized at 1600 °C) from the referred study was chosen due to its electrochemical performance [16]. This HC1600 electrode sample was prepared from the AM powder by creating an aqueous slurry including 92 wt% HC, and 8 wt% of additives (in this case including binder, conductive carbon and surfactant). The HC electrodes were not calendered, and their porosity was calculated to be around 60 %. The AM powder consisted of platelet-like particles of different sizes, ranging from 5 to 13 µm in length, and from 2.5 to 5 µm in thickness (as reported in the Supplementary Information of the referred work).

For electrochemical testing, the authors assembled CR2032 coin cells, using HC and sodium metal electrodes with a glass fiber as separator and an electrolyte composed of sodium hexafluorophosphate salt (NaPF₆, 1 mol L⁻¹) inside a mixture of ethylene carbonate and dimethyl carbonate (EC/DMC, DodoChem, 50:50 by mass). These cells were cycled galvanostatically vs. Na⁺/Na at C/20 for 10 cycles, and at C/10 for 100 cycles, (18.6 mA g HC⁻¹ for C/20). The applied current was calculated using the theoretical capacity of NaC₆ (372 mAh g HC⁻¹).

## 3. Model description

To ensure that the electrochemical simulations reflected the real electrode characteristics and experimental conditions, the model was parameterized using data from the earlier experimental work [16]. The data consisted of manufacturing parameters of the HC1600 electrode sample (including formulation, solid wt% of AM and additives in the slurry, and electrode porosity), as well as characterization data of the AM (like particle shape and size, density, and textural properties), up to the electrochemical data (cycling protocol, experimental capacity, and discharge profiles). For the electrochemical data, information of HC vs. sodium metal was used to parameterize and validate the model. Since this model neglected the formation of Solid Electrolyte Interphase (SEI) on HC during the first cycle, the electrochemical data was taken from the second cycle. We have also focused on modeling the two C-rates (C/20 and C/10) reported in the experimental work [16], as to have a direct comparison. The general workflow for this work is outlined in Fig. 2.

### 3.1. Stochastic generation of electrode microstructures

To consider the effect of HC electrode heterogeneities on the electrochemical model, 3D electrode microstructures were generated by using the GrainGeo module of the commercial software GeoDict®, version 2024 (Math2Market) and a standard desktop computer [44]. The generated electrode microstructures consisted of the open macro-pore domain, as well as two solid domains corresponding to the AM and the CBD, the latter one being a phase that aggregates the conductive carbon additive and the binder, and in this case, also the surfactant. This
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Schematic of our modeling workflow. From the injection of the experimental data (manufacturing, characterization, and electrochemical), to the generation of the 3D-resolved electrode microstructure and electrochemical input parameters, until the FEM electrochemical simulations.</figcaption></figure>

consideration of an aggregated phase has been previously discussed by our group in earlier works for microstructure-resolved electrode manufacturing models of LIBs [34–38]. Here, we are demonstrating that this concept can also be transferred to SIBs.

The generation of the HC electrodes was done by stochastically filling a computational box of 30 × 30 × 60 µm³ with particles of each solid phase to target a formulation of 92 % AM and 8 % CBD, while considering an overall electrode porosity of 60 %, as in the experimental work [16]. The stopping criterion for the simulated electrodes was a target value of Solid Volume Percentage (SVP), which was calculated from experimental weight percentage inputs and density values for each one of the solid phases. The target values of SVP were calculated considering the internal porosity of the solids, and were found to be 35.15 % and 4.85 % for the AM and CBD phases, respectively. The HC internal nanoporosity was calculated to be 4.85 %, according to experimental gas adsorption data, and assumed to be 50 % for the CBD [34]. Since these values had to be included in the overall electrode porosity, when calculating the SVP from the wt% reported in the experimental work, the density of the materials was modified to include the effect of the internal porosity of each phase, as shown in previous studies [45]. For our particular case, considering that HC has a density of 1.5 g cm⁻³

[16] and CBD one of 1.8 g cm⁻³ [36], incorporation of the aforementioned nanoporosities led to modified values of 1.42 g cm⁻³ and 0.9 g cm⁻³, respectively. The detailed calculations of the SVP, and HC nanoporosity are explained with detail in equations S1-S8, and S9-S12 of the Supplementary Information.

For the generation of the electrode microstructures, the AM particles were the first ones to be stochastically placed inside the simulation box, using the “Create Grains” function of the GrainGeo module of GeoDict®, with different random seeds for each microstructure. The particles were chosen to be non-uniform and non-spherical, with a planar polyhedron shape, according to scanning electron microscope measurements [16]. The parameters used for the stochastic generation of the AM particles are shown in Table S1 of the Supporting Information. After having the AM particles in the simulation box, the CBD particles were included by selecting the “Add Binder” function, also from the GrainGeo module of GeoDict®. By using this function, the CBD phase gets placed in the shape of a concave meniscus where the software detects agglomerated particles of AM. For this step, the default value of 0° was chosen for the contact angle of the CBD. A 3D visualization of the final domains included in the stochastic 3D structures can be seen on Fig. S1 of the Supplementary Information.
After the electrode microstructures were generated in GeoDict®, they were further exported as 2D image stacks into INNOV (our in-house meshing tool developed using Matlab) with the purpose of generating a meshed version of the HC microstructure, using a 0.25 µm voxel resolution [46]. The benefit of using this tool is that, in addition to meshing, a homogeneous separator and Current Collector (CC) can also be included. The thicknesses of the separator and CC were chosen to be 10 µm in height.

After the meshing, the electrode microstructures were imported into the electrochemical model, for performing the FEM computations. The number of elements, resulting from the meshing step, is shown in Table S2 of the Supporting Information for all electrode microstructures generated during the development of this work.

## 3.2. Potential calculation of Intercalation-Filling mechanism

In order to obtain insights into the sodiation mechanism of HC and decouple the contributions of the two sodiation processes, we took previous electrochemical models for blended electrodes as an inspiration [38,47]. For this type of model, electrochemical inputs for the different AMs that make up the blended electrode, are necessary.

Following the hypothesis that the HC electrode behaved as described by the sodiation mechanism proposed by Stevens and Dahn [6], the HC electrode was treated as if it was composed of two different AMs, one that stores sodium through an intercalation process and one that does so by undergoing the nanopore filling process. This means that our model assumes that during the HC discharge versus Na metal, sodium is stored through intercalation (at the sloping potential region), followed by a subsequent sodium storage in the HC nanoporosity (at the plateau potential region).

To assess if the potential of our HC sample follows the Intercalation-Filling model potential, a decoupling of the electrochemical signature measured experimentally into two different input equilibrium potentials was performed. To evaluate if these decoupled potentials agreed with the overall HC behavior, we used blended electrode formulas as to evaluate the global potential resulting from the mixture of the two AMs. The following calculation was used to assess the necessary concentration of Na to reach the global potential, in this case, the potential value measured experimentally. This calculation was done inside a working potential window comprised between the minimum potential ($U_{min}$) and the maximum one ($U_{max}$), which in this case, were calculated from the individual sodiation phenomena as shown in equations (1) and (2),

$$U_{min} = max_i \left( \min \left( U_{ref}^i \right) \right) \qquad (1)$$

$$U_{max} = min_i \left( \max \left( U_{ref}^i \right) \right) \qquad (2)$$

where $i$ represented the sodiation phenomenon, either intercalation or nanopore filling. To ensure a good range of potential we set the condition for $U_{ref}^i(0) > 2$ and $U_{ref}^i(1) \simeq 0$. In this model, $N_{global}$ values of the resulting global potential of the blended electrode $U_{global}^k$ are taken regularly in the interval inside the range [$U_{max} \cdot U_{min}$] as the potential decreases with the sodiation extent. From the initial concentration values of a phenomenon $i$ at $U_{max}$ or $U_{global}^0$, we described this by the expression of equation (3),

$$c_s^i = U_{ref}^i \cdot U_{global}^0 \cdot c_{s,i}^{max} \qquad (3)$$

Having established this, the Na concentration needed for each phenomenon $i$ to reach $U_{global}^k$ could then be obtained using equation (4),

$$c_s^k = U_{ref}^i \cdot U_{global}^k \cdot c_s^{i,max} - c_s^0 \qquad (4)$$

This previous consideration was necessary, since it allowed to start evaluating potential values which do not correspond to a concentration

of zero in all of the sodiation phenomena. The concentration in the material could then be calculated by equation (5),

$$c_s^{to} = \sum_i c_s^i \qquad (5)$$

At the end, the total maximum concentration that can be loaded into the material, is the maximum value of $c_s^{to}$. The result of this calculation was the full Intercalation-Filling potential as a function of the total concentration, thus reflecting the global behavior of the HC electrode, which considered the two sodiation phenomena. This calculated potential could then be compared to the experimentally measured one, as to ensure that the two decoupled potentials reflected the behavior of the HC material.

## 3.3. Electrochemical model

The electrochemical model was built using the experimental half-cell configuration, considering the following domains in the electrode's 3D microstructure: AM, CBD, electrolyte, separator, CC and Na metal, the latter one modeled implicitly through the use of a boundary condition, as shown in the lower part of Fig. 2. In our particular case, resolving nanometric subdomains to represent each sodiation processes, while keeping the electrode microstructure in the order of microns, would be computationally prohibitive. For this reason, our model does not explicitly resolve a domain for each sodiation phenomenon and instead, the AM is treated as an homogenized volume. In this way, the two phenomena are described mathematically within the same AM domain, by using a corresponding solid volume fraction, according to the hypothesized amount of each phenomenon in our particular HC. As a reminder, this model also considered the internal nanoporosity of HC implicitly, by modifying the parameters of the microstructure generation, previously described in Section 3.1. In addition, the CBD, and separator phases were also considered as homogeneous domains with assumed electrolyte volume fractions of 0.5, and 0.9, respectively.

The electrochemical model was built using the version 6.2 of COMSOL Multiphysics® software, all calculations were run on the MatriCS platform (Université de Picardie Jules Verne) using one node, 500 GB of RAM, and the Intel® Xeon® CPU E5-2680 processor, with 28 cores. A multifunctional massively parallel sparse direct solver was chosen to solve the Partial Differential Equations (PDEs) using a backward differentiation formula method for time stepping. The simulation time varied with the applied discharge rate, and on average, the calculations took 5, 4 and 3 h for simulating C/20, C/10 and C/5 rates, respectively.

The equations used for the model are based on the “Doyle-Fuller-Newman” model equations [48–50], extended to describe 3D geometries, like shown in previous works [34,37,38] and using different parameters for each sodiation process. Table S3 of the Supporting Information presents a full description of these equations.

The model accounts for the electronic conduction in the CBD and AM domains, using Ohm's law and conductivity values of 100 S m⁻¹ and 0.48 S m⁻¹, respectively. The ionic charge transport is considered in the bulk electrolyte, using concentrated solution theory, and an electrolyte conductivity value of 1.02 S m⁻¹. Likewise, the ionic charge transport is also considered to happen in the implicit pores of the separator and CBD phases, and corrected for by using the Bruggeman relation, alongside the solid volume fraction of the phases, and a Bruggeman coefficient of 1.5. In an analogous way, the Na⁺ transport is also considered in the bulk electrolyte, separator, and CBD phases, because these last two phases are permeable to Na⁺ ions, due to their porosity. The sodium transport happening from the surface of the HC to the center of the homogenized HC domain is considered for both sodiation phenomena, through the use of Fick's law, and using a solid diffusion coefficient for sodium in HC of 1 x 10⁻¹² m² s⁻¹, this value was obtained through optimization approaches. This value was used as is for modeling the intercalation phenomenon, but modified for modeling the transport of sodium related to
the nanopore filling process by using the Bruggeman relation, the internal nanoporosity of the HC domain, and a Bruggeman coefficient of 1.5; resulting in a value of $1 \times 10^{-14} \text{ m}^2 \text{ s}^{-1}$.

In this model, the intercalation and nanopore filling reactions are considered as two distinct charge transfer reactions which happen at the AM/electrolyte boundary. At the surface of HC domain, mass transport of sodium is coupled to the flux of each reacting species, and in this way two concentration variables can be tracked along the model (the solid concentration of the intercalation, and filling reactions, $c_{s,int}$ and $c_{s,fill}$, respectively). The kinetics for the two reactions are modeled through the use of Butler-Volmer expressions, also one for each considered sodiation phenomenon. Individual exchange current densities for the intercalation ($I_{0,int}$) and nanopore filling ($I_{0,fill}$) reactions are established, expressed in terms of each specific concentration variable.

Due to the porous nature of the CBD and separator domains, the intercalation and nanopore filling reactions were also considered to have an effect at the AM/CBD and AM/separator boundaries, for this reason, they are considered to be permeable to sodium, as it has been done in previous works for LIBs [34]. In this case, the kinetics of these two sodiation reactions were adjusted by using their electrolyte volume fractions to correct each reaction's exchange current density, as denoted in Table S3. For the AM/electrolyte boundary, in addition to the charge transfer reactions, we considered the formation of an electrical double-layer capacitance, while for the AM/CBD and AM/separator boundaries, only the intercalation and nanopore filling reactions were considered.

At the external CC and Na foil boundaries, boundary conditions were considered to control the applied discharge current ($I_{app}$), and set the Na metal potential to 0 V, respectively. In all other external boundaries of the simulation box, an insulation and a no mass flux boundary conditions were applied.

In the end, this electrochemical model solves for five dependent variables: electrolyte potential ($\psi_l$), electric potential in the solid ($\psi_s$), salt concentration in the electrolyte ($c_l$), and the two solid Na concentrations in the AM, corresponding to intercalation and nanopore filling.

# 4. Results and discussions

## 4.1. Decoupling intercalation and nanopore filling potential contribution

As described in Section 3.3, our model deals with two different sodiation reactions, so it requires two distinct equilibrium potential inputs for both intercalation and filling processes, ($E_{eq,int}$ and $E_{eq,fill}$, respectively). Since we only had a single electrochemical response when cycling HC vs. Na, we propose a method to assign two individual functions to the experimental data as to obtain two different equilibrium potential expressions (in the form of E versus Na fraction). These expressions were derived by using the experimental discharge profile at C/20. These functions were obtained using an in-house Python code and a piecewise cubic interpolation approach, by the CubicSpline interpolator from the SciPy library.

Since the hypothesis that we wanted to test in this first model accounting for a two-step sodiation mechanism was the Intercalation-Filling theory [6], the code minimizes two cubic spline functions until reaching the best fit to the experimental data in an specified subinterval. In this case, the slope and the plateau sections of the discharge profile. The way to specify the end of the first subinterval and the start of the second one was by defining an internal parameter, referred to as $x_{break}$. Different values for $x_{break}$ where explored, ranging from 0.10 to 0.40. In the end, the best fit was obtained by using $x_{break} = 0.30$. This obtained value of this parameter was related to the percentage of each sodiation phenomena inside this particular HC material, and utilized to select the volume fraction of each phenomenon in the electrochemical model (in this case 0.3 and 0.7 for intercalation and nanopore filling, respectively). If the methodology described in this paper were to be applied to a different HC material, this hyperparameter would need to be optimized

to the particular electrochemical behavior of that AM.

The obtained potentials as a function of Na fraction for each sodiation process are shown in Fig. 3a, along with the experimental data to which they were fitted to. To be sure that these two potentials were relevant to the electrochemistry of the HC electrode, the global Intercalation-Filling potential of the HC electrode was calculated by adding the contribution of these two inputs with the approach described in Section 3.2. The obtained Intercalation-Filling potential model was compared to the experimental data and shown in Fig. 3b, where a good agreement is observed.

Due to the fact that our decoupling approach resulted in two equilibrium potentials with very different starting potentials (in this case 1.3 V and 0.05 V), inputting these into the electrochemical model failed due to a lack of numerical convergence. This failure was expected, and associated with the extreme current redistribution that occurred during the first steps of the model, something hard to handle for a numerical solver. To solve this problem, the obtained $E_{eq,int}$ and $E_{eq,fill}$ potentials were modified to be in the same range on both x and y axes at the initial Na concentration. To account for this, two starting points were introduced to the equilibrium potential expressions obtained from the Python code. These points corresponded to the potential at zero and at the initial state of charge, the first being located at 2 V and the latter corresponding to the potential at which the experimental data was obtained (either 1.3 or 1.0 V, for the C/20 and C/10 discharge rate, respectively).

It should also be noted that the expression for $E_{eq,int}$ was rapidly decreasing to zero after the Na state of charge reached 0.3 ($x_{break}$), not fully capturing the effect of the slope at the start of the discharge, this is why the $E_{eq,int}$ potential was normalized by the last value before the plateau. The final equilibrium potential inputs used in our electrochemical model built on COMSOL Multiphysics®, are shown in Fig. S2 of the Supporting Information.

## 4.2. Model parameterization and validation

A full description of the optimized parameters of our electrochemical simulations is provided in Table S4 of the Supplementary Information. A sensitivity analysis was performed to either optimize or take some parameters from the literature, in the few cases where they could not be obtained directly or calculated from the referred experimental work [16]. For our model, we identified the kinetic rate constants ($k_{int}$ and $k_{fill}$), and the solid-state diffusion coefficient ($D_{HC}$) as the most sensitive parameters, so we adopted a sequential approach to optimize them. First, we explored a wide range of value combinations for $k_{int}$ and $k_{fill}$, ranging from $1 \times 10^{-15}$ and $1 \times 10^{-10} \text{ m s}^{-1}$. We evaluated numerous combinations of values in three scenarios: first, when $k_{int}$ was of the same order of magnitude as $k_{fill}$, second, when $k_{int} < k_{fill}$, and lastly, when $k_{int} > k_{fill}$. The best results were obtained for the last scenario, with the optimal values found to be $9 \times 10^{-11}$ and $8 \times 10^{-12} \text{ m s}^{-1}$ for the intercalation and filling kinetic constants, respectively. After fixing the kinetic parameters, we optimized the $D_{HC}$ by testing several values between $1 \times 10^{-16}$ and $1 \times 10^{-10} \text{ m}^2 \text{ s}^{-1}$. For this parameter, we determined that $1 \times 10^{-12} \text{ m}^2 \text{ s}^{-1}$ was the optimal value. To exemplify the optimization process, a schematic figure showing the optimization of the kinetic parameters can be found in the Supplementary Information (Fig. S3). To decrease the computational cost, the sensitivity analysis was performed on smaller versions of the electrode microstructures ($15 \times 15 \times 30 \text{ µm}^3$) and then transferred to the larger ones used along the work.

The parameterization was performed with the discharge profile obtained at C/20. Then, this parameterization was validated by directly comparing the experimental discharge profile obtained at C/10. Both parameterization and validation were performed in the same potential window as in the experimental study. All the results presented in this work were obtained using this same set of parameter values for all discharge rates at a constant temperature of 298.15 K. For statistical purposes, three distinct stochastic electrode microstructures (with
<figure><img src="image_4.png" /><figcaption>Fig. 3. a) Decoupled intercalation (blue) and filling (red) potentials obtained through fitting of the experimental discharge profile data at C/20 (black). b) Experimental C/20 data (black) shown in relation to the calculated Intercalation-Filling potential (purple) obtained through the consideration of the individual potentials, as described in Section 3.2 (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article).</figcaption></figure>

identical target SVP criteria, but different initial random seeds) were generated and evaluated with the electrochemical model at both C-rates. The results of the three configurations were averaged out, presenting the average result along with the standard deviation. The three different stochastic configurations are shown for all simulated scenarios in Fig. S4 of the Supplementary Information. The applied discharge current was calculated from the mass of the stochastic HC microstructures, and using a theoretical capacity of 372 mAh g⁻¹, as reported in the experimental study [16]. Table S5 of the Supporting Information depicts the specific discharge currents applied to each structure.

Fig. 4 shows a comparison of the experimental (scatter points) and mean simulated (solid lines) discharge profiles, as a function of the obtained specific capacity for the electrodes subjected to a C/20 (Fig. 4a) and C/10 (Fig. 4d) rate. The standard deviation of the three simulated microstructures is shown in the plots as a darker color around the solid lines. The thin aspect of this standard deviation line indicates that there is a minor variation in the results of the three generated stochastic

configurations.

For the C/20 case, the mean simulated discharge falls between 1.28 V and ends at 0.32 mV. It can be observed that the model is able to capture the main features and trends of the experimental data. Nevertheless, the simulated discharge at C/20 slightly overestimates the potential at the intercalation zone (from the start of discharge and up to 0.7 V). On the other hand, the potential values at the plateau zone are better matched to the experimental dataset after 100 mAh g⁻¹, however, at the transition zone (around xbreak), the model does not fully capture the curvature, indicating the switch from one sodiation process to another.

For the C/10 rate, the mean simulated discharge starts at 0.97 V and ends at 0.67 mV, providing an overall better match to the experimental dataset. Nevertheless, some overestimation (smaller than in the C/20 case) can still be seen for the intercalation zone, observed from the start of the discharge until 0.2 V.

In terms of the specific capacity obtained at the electrode level, the

<figure><img src="image_5.png" /><figcaption>Fig. 4. (a) Comparison of experimental versus simulated discharge profiles for C/20 rate. (b) Evolution of the concentration variables in relation to the simulated capacity for C/20 rate. (c) The spatial distribution of sodium concentration variables in one of the 3D electrode microstructures for C/20 rate. (d) Comparison of experimental versus simulated discharge profiles for C/10 rate. (e) Evolution of the concentration variables in relation to the simulated capacity for C/10 rate. (f) The spatial distribution of sodium concentration variables in one of the 3D electrode microstructures for C/10 rate.</figcaption></figure>
simulations at C/20 and C/10 reached 254.82 mAh g⁻¹ and 226.92 mAh g⁻¹, respectively. When compared to the experimental reference values (291 mAh g⁻¹ and 235 mAh g⁻¹, for C/20 and C/10 respectively), both simulations underestimate them. The discrepancy between the simulated and experimental capacities may be partly due to numerical issues arising from the solver, particularly in the plateau region where values approach zero, since this can interfere with the stability and accuracy of solving PDEs using numerical methods like FEM and commercial software like COMSOL Multiphysics®. Additionally, the discrepancies may suggest the presence of a supplementary sodiation mechanism in the experimental sample, that is not able to be captured by the current model, such as the adsorption of Na⁺ on defect sites or heteroatoms.

Fig. 4b and e show the evolution of the concentration variables in relation to the simulated capacity for C/20 and C/10 rates, respectively. These concentration variables represent the sodium stored in AM domain in the form of either intercalation (c<sub>s,int</sub>) or nanopore filling (c<sub>s,fill</sub>) processes. Moreover, the global sodiation contribution can be tracked by summing both concentrations and plotting the total concentration variable (c<sub>s,tot</sub>). For the C/20 and C/10 cases, the highest values of c<sub>s,tot</sub> were 14400 and 12842 mol m⁻³, respectively. In contrast,

the maximum concentration calculated from the experimental capacity of the material was 16286.41 mol m⁻³, and not reached by either of the simulations. For the C/20 simulation, the c<sub>s,int</sub> and c<sub>s,fill</sub> values reached 4281.13 and 10119 mol m⁻³, respectively. On the other hand, for C/10 the c<sub>s,int</sub> and c<sub>s,fill</sub> values were 4152.80 and 8689.53 mol m⁻³, respectively. We can also comment that the nanopore filling contribution at the End of Discharge (EoD) was, on average, 70 % of the total concentration in the C/20 example, compared to 68 % in the C/10 case.

When comparing both Fig. 4b and e, for both of the C-rates studied, the filling phenomenon is stationary at the beginning of the discharge, as indicated by the c<sub>s,fill</sub> plateau. Nonetheless, nanopore filling begins to play a substantial role and increases rapidly after surpassing 55 mAh g⁻¹. For the intercalation behavior, we can observe a consistent growth since the start of the discharge. However, after reaching 55 mAh g⁻¹, the c<sub>s,int</sub> variable enters a stage where it reaches a small plateau and then starts to slowly increase. For the two rates studied, there is a point in the evolution of concentration where the c<sub>s,int</sub> and c<sub>s,fill</sub> curves cross each other (around 100 and 110 mAh g⁻¹ for the C/20 and C/10 cases, respectively).

Fig. 4c and f depict the spatial distribution of the total, intercalation,

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** On top, the spatial distribution of the solid potential (φ<sub>s</sub>) in the AM, CBD and CC regions for a) C/20 and b) C/10 rate. In the bottom, the spatial distribution of the liquid potential (φ<sub>l</sub>) in the electrolyte and separator domains for c) C/20 and d) C/10 rate.</figcaption></figure>
and nanopore filling sodium concentrations inside the homogenized HC domain for the C/20 and C/10 rates, respectively. At C/20 (Fig. 4c) the concentration is distributed evenly across the microstructure, for all three cases. However, when subjected to a rate of C/10 (Fig. 4f) we can notice a decreased $c_{s,tot}$ concentration, indicated by the pink areas in the structure showing the $c_{s,tot}$ distribution. More importantly, we observe the presence of increased concentration gradients for the $c_{s,fill}$ variable, which contributes with most of the total concentration, and which influence is reflected on the concentration gradients observed in the $c_{s,tot}$ distribution.

Fig. S5 shows the concentration of Na⁺ in the electrolyte for both C-rates. This concentration remained evenly distributed around the starting values of 1000 mol m⁻³, with some slight concentration decrease (around 997 mol m⁻³) in the current collector side, for the C/10 case. This behavior can be attributed to the slow nature of both studied C-rates, where the ionic transport is not a limiting process.

Fig. 5 shows the results of the obtained potential variables. More specifically, in Fig. 5a and b, we can observe the φs distribution at EoD in the solid parts of the electrode microstructure (AM, CBD and CC) for the C/20 and C/10 case, respectively. We identify that the solid potential is higher on the AM close to the separator, and lower on the AM close to the CC. On Fig. 5c and d we can observe the distribution of φl in the domains with liquid (bulk electrolyte, and separator and CBD), for C/20 and C/10, respectively. The same trend of higher potentials close to the

separator side, and lower potentials on the CC side is maintained. Overall, the gradients of both potential variables augment when a higher C-rate is applied, which is an expected behavior.

### 4.3. Influence of manufacturing parameters on electrode level performance

In addition to generating 3D electrode microstructures that matched the experimental data, consisting of a 92:8 AM:CBD formulation with 60 % porosity, these two manufacturing parameters were changed in order to generate two new electrode microstructures and later assess their effect on the electrochemical performance at the electrode level. The choice of the new manufacturing parameters was based on previous experimental works reported in the literature. First, a different formulation of 95:5 AM:CBD [51] was used, while maintaining the original value of porosity at 60 %. Second, the formulation of the original electrode (92:8 AM:CBD) was kept, but the value of the overall electrode porosity was lowered to a value of 41 % [52]. For readability sake, these two new electrodes are referred to as by 95.60p and 92.41p, respectively, while the original electrode is denoted as 92.60p. In all these cases, the first number indicated the AM content and the second one the overall electrode porosity. These additional electrode microstructures were generated using the same stochastic generation protocol described in Section 3.1. It is important to clarify that in the case of the 92.41p

Experimental
Formulation + Porosity

92:8 AM:CBD
60% porosity

a)

<figure><img src="image_7.png" /></figure>

b)

<figure><img src="image_8.png" /></figure>

Formulation
change

95:5 AM:CBD
60% porosity

c)

<figure><img src="image_9.png" /></figure>

Porosity
change

92:8 AM:CBD
41% porosity

<caption>**Fig. 6.** Electrode microstructures showing the difference in phase distribution for the three different cases simulated: a) Experimental electrode sample “92_60p”, b) Additional electrode sample “95.60p”, and c) Additional electrode sample “92.41p”. CBD phase is shown in colors purple, orange and teal, for the 92_60p, 95_60p and 92.41p samples, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)</caption>
electrode, the redistribution of the CBD phase resulting from electrode calendering is not captured, since the use of our physics-based manufacturing simulation (the ARTISTIC framework, previously demonstrated for LIBs) was not used in this work [29]. Instead, the CBD was stochastically added on top of the AM particles, as described in Section 3.1. Fig. 6 compares the three electrode cases analyzed: 92_60p, 95_60p and 92_41p.

It is necessary to say that even though no experimental raw data was available to directly compare these new electrode simulations, the intent of this section was to take advantage of our model (previously optimized for a specific type of HC), and show that it can be used to model different HC electrode configurations made with the same AM but different formulation parameters, thereby optimizing the electrode design stage.

The new electrode microstructures (95_60p and 92_41p) were imported into our electrochemical model and galvanostatic discharges at

three different C-rates were performed: C/20, C/10 and finally, C/5. Fig. 7 displays the electrochemical performance findings of the two new electrode microstructures evaluated at discharge rates of C/20 and C/5, as the most contrasting cases. With the purpose of comparing the overall capacity performance of the three simulated cases at three different discharge rates, the original electrode microstructure was also subjected to a C/5 discharge rate (although this discharge rate was not tried in the experimental study). In Fig. S6 of the Supplementary Information, we include the performance of the simulated electrode cases when subjected to a C/10 discharge rate.

In terms of discharge profiles (Fig. 7a, c, e and g), the new electrodes behave quite similarly to the model shown in section 4.2. However, a slight decrease in specific capacity is obtained for the new manufacturing conditions tested.

In terms of the spatial distribution of Na in the HC phase, the trend

<figure><img src="image_10.png" /><figcaption>Fig. 7. a) Simulated discharge profiles and evolution of the concentration variables in relation to the simulated capacity of the 95_60p electrode microstructure, at a C/20 rate. b) Spatial distribution of sodium concentration variables in 95_60p electrode microstructure at a C/20 rate. c) Simulated discharge profiles and evolution of the concentration variables in relation to the simulated capacity of the 95_60p electrode microstructure, at a C/5 rate. d) Spatial distribution of sodium concentration variables in 95_60p electrode microstructure at a C/5 rate. e) Simulated discharge profiles and evolution of the concentration variables in relation to the simulated capacity of the 92_41p electrode microstructure, at a C/20 rate. f) Spatial distribution of sodium concentration variables in 92_41p electrode microstructure at a C/20 rate. g) Simulated discharge profiles and evolution of the concentration variables in relation to the simulated capacity of the 92_41p electrode microstructure, at a C/5 rate. h) Spatial distribution of sodium concentration variables in 92_41p electrode microstructure at a C/5 rate.</figcaption></figure>
observed for the 92_60p electrode (Fig. 4) was maintained for the new simulated electrodes, this meant that greater concentration gradients are caused by faster discharge rates, particularly impacting the c<sub>s,fill</sub> distribution. Another observation is that even if the 95_60p and 92_41p simulations used the same AM parameters as the experimental 92_60p electrode, the electrochemical performance was reduced for the two additional electrode cases, simply because of the microstructural changes in the electrode. The lowered performance from the 95_60p electrode in this case was associated with the decreased CBD quantity in the formulation (only 5 %), which resulted in poor electronic conduction channels. For the specific microstructure shown in Fig. 7b and d, we can observe that the CBD does not connect well to all the AM particles, as there are zones left completely unsodiated, due to this poor connection. The lowered performance of the 92_41p electrode on the other hand, was associated to the 19 % decrease in porosity, which hinders the sodium transport through the bulk electrolyte phase and making the overall electrode to suffer from higher concentration gradients, as can be seen in Fig. 7h. For illustrative purposes, animations showing the spatial distribution of the c<sub>s,tot</sub> variable are shown for each simulated electrode case, at three different C-rates (see Videos S1-S3 of the Supplementary Information).

### 4.4. Electrode microstructure properties and performance comparison

With the purpose of analyzing the electrochemical simulation results in conjunction to the electrode microstructure metrics, we first compared the final capacity obtained for each simulated case. Fig. 8a shows the mean specific capacity obtained at EoD for the three electrode microstructures tested across different C-rates and electrode configurations. For all electrodes, the slowest discharge rate (C/20) resulted in the highest capacity, while the fastest one (C/5), yielded the worst performance. In addition, error bars are higher for the C/5 rate, indicating that with increasing rates, the model can experience a greater deal of variation. At C/20, the performance of the three simulated scenarios remained similar. However, differences in performance become noticeable when the electrodes are subjected to C/10 and C/5 rates, where we can start to see the impact of changing formulation and porosity values. In this comparison, we also observe a decrease in capacity going from 92_60p to 95_60p to 92_41p, indicating that a decreased porosity has a more detrimental impact in performance than a reduced CBD content.

Another way to analyze the obtained performance is through careful analysis of the geometrical features of the electrode microstructures, mainly tortuosity factor and electronic conductivity. These properties were calculated for all electrode cases by using GeoDict<sup>®</sup> in a standard desktop computer, and in a similar way as some of our previous works have detailed [53]. The “Simulate Diffusion Experiment” function was selected from the DiffuDict module in order to obtain the effective diffusivity for each electrode microstructure. This was done by using stationary state diffusion calculations and solving the Fick’s laws throughout the macropore network. This value of effective diffusivity was then used to calculate the electrode’s tortuosity factor (τ), which is determined as the ratio between porosity and the effective diffusivity (τ = ε/D<sub>eff</sub>). For the electronic conductivity, the “Compute Electrical Conductivity” function of the ConductoDict module was used. Here, we considered the individual conductivity values of our AM and CBD phases (0.48 S m<sup>-1</sup> and 100 S m<sup>-1</sup>, respectively), and Ohm’s law as well as the Poisson equation, to simulate the conductivity experiment. For both metrics, the calculations considered the use of periodical boundary conditions in the tangential directions (x and y), while considering a Dirichlet condition in the z direction.

Fig. 8b shows a comparison of the mean value for the tortuosity factor and electronic conductivity at the electrode level for the three simulated scenarios. These two metrics were calculated from the three stochastic electrode microstructures utilized for each simulated case. We can see that in comparison to the electrode microstructure based on experiments (92_60p), the case with increased AM content (95_60p) showed reduced values for both analyzed electrode metrics, however, for the case where the overall porosity is reduced (92_41p), both the tortuosity factor and the electronic conductivity increased.

Thanks to this electrode level analysis, we confirmed that the 95_60p and 92_41p cases exhibited a lowered electrochemical performance due to two different reasons. First, even if a higher AM content was expected to improve the capacity of the electrode, the lower CBD content (for the 95_60p case), decreased the overall electrode’s electronic conductivity, due to the loss of electronic percolation pathways. And second, even if a higher electrode electronic conductivity was achieved by compacting the electrode, the lower macroporosity generates more tortuous paths for Na<sup>+</sup> transport in the electrolyte located in the macro-pores, which results in ionic transport limitations (for the 92_41p case). These findings highlight the need to carefully optimize the formulation and porosity through a careful electrode design.

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Comparison of the mean and standard deviation of a) electrode level capacity at EoD and b) tortuosity factor and electronic conductivity, obtained for the three different simulated electrode cases: 92_60p, 95_60p and 92_41p.</figcaption></figure>
## 5. Conclusions

In this work, we have presented a new physics-based electrochemical model for HC electrodes applied to SIBs, in a half cell configuration. This model targets the study of the Intercalation-Filling sodiation mechanism by tuning both the kinetics and transport aspects of each sodiation phenomenon for a particular HC sample. The electrochemical and transport processes occurring in the HC are modeled with 3D resolution at the mesoscale, using a homogenized approach for the two sodiation processes. The model used inputs of 3D electrode microstructures, which were stochastically generated. These microstructures explicitly considered three phases: AM, CBD and pores. All these aspects make this model the first of its kind for application in SIB technology.

First, the 3D electrode microstructures were stochastically generated through the use of the commercial software GeoDict®. The electrodes were then meshed using our in-house Matlab software, INNOV. The model was parameterized using experimental AM data like particle shape and size, density, and textural properties. In addition, the discharge profiles at C/20 and C/10 were used to both parameterize and validate the model, respectively. The experimental HC electrode configuration with a formulation of 92:8 AM:CBD, a porosity of 60 %, and an electrolyte of NaPF₆ salt in a mix of 50:50 EC/DMC was successfully modeled. The electrochemical performance was evaluated through FEM computations performed using COMSOL Multiphysics®. Three discharge profiles were obtained from the different stochastic microstructures for each simulation, and the average of these profiles was calculated and plotted against the experimental data. For the simulations performed at C/10 rate, the simulated profile closely matched the experimental data, while at C/20, the model overestimated the voltage values at the intercalation region and underestimated the final capacity. The deviation between the experimental and simulated discharge profiles can be attributed in part to the fact that not all of the known sodiation processes were considered in this proof-of-concept model. In addition, the discrepancy between the simulated capacity at EoD might also be influenced by numerical challenges, like solving PDEs involving values close to zero, which may cause localized differences when modeling the potentials in the plateau region.

The 3D-resolved electrode character of the model has been useful to demonstrate the transferability of our ARTISTIC mesoscopic modeling approach to study the electrode level behavior of the HC material, with electrode microstructures arising from different manufacturing parameters. We tested first an electrode with a different formulation (95:5 AM: CBD), but the same experimental porosity (60 %), as well as an electrode with different porosity (41 %), while keeping the same experimental formulation. In this case, these complementary simulations have shown that for the 95_60p case, increasing the HC content while decreasing the CBD content, does not increase the electrode level performance. This can be attributed to the loss of electronic percolation, that eventually leads to a loss of efficiency, highlighting the important role of the CBD phase and its distribution. And in the 92_41p case, decreasing the porosity while keeping the same formulation, increases the electronic conductivity of the electrode, but at the same time results in a decreased performance due to a higher tortuosity factor, hindering the performance because of increased ionic transport limitations. The use of our 3D model to exhaustively study the influence of the manufacturing parameters of SIBs (not only limited to porosity and formulation), is envisioned as one potential use of our model, and represents a promising direction for future studies emerging from this approach.

As future perspectives, several directions can be suggested for an improved version of the model. The first is the adoption of more representative HC electrode architectures, which can be derived from physics-based manufacturing simulation approaches, like the one pursued in our ARTISTIC digitalization initiative [29], or experimental imaging techniques, such as computer tomography or focused ion beam scanning electron microscopy. A second approach could be to improve the parameterization of the kinetic and transport properties,

particularly, using data from galvanostatic intermittent titration technique, for a more tailored approach to the specific HC modeled, introducing parameters as functions of the electrode's state of charge. The assumption of the porosity value of the porous phases (like the CBD), could also be improved as to be determined with new methodologies that combine both experimental and simulation techniques, like shown by us in a recent work [45]. Third, the inclusion of additional sodiation phenomena (like the adsorption of Na⁺ on defect zones, or hetero-atoms) could also enhance the prediction of both potential and specific capacity, as this consideration may allow the model to describe a more complex sodiation mechanism composed of multiple stages. Lastly, it would also be beneficial to expand the model to investigate a broader range of C-rates, particularly for studying SEI formation (at low C-rates), or testing the rate performance of electrodes (at higher C-rates). The inclusion of aging phenomena and its effect on the performance indicators is an area of opportunity and targeted as a next step for our work.

Still, we believe that the novel model proposed in this work paves the way towards a more comprehensive understanding of the interplay between electrode microstructure, transport, reaction processes and the resulting electrochemical performance of the HC. Our modeling approach is intended to support the efficient engineering and design of this key electrode material for SIB technology. Inspired by the tremendous impact of Professor Bruno Scrosati's work [54] whose scientific contributions over several decades were instrumental in shaping the development of modern battery technology, we hope that this work humbly honors his legacy.

## Materials availability

This study did not generate new unique materials.

## CRediT authorship contribution statement

**Imelda Cardenas-Sierra:** Investigation, Writing – original draft, Visualization, Validation, Software, Methodology, Formal analysis, Data curation, Conceptualization. **Francisco Fernandez:** Writing – review & editing, Visualization, Software, Methodology, Data curation, Conceptualization. **Martin Petit:** Writing – review & editing, Supervision. **Alejandro A. Franco:** Supervision, Resources, Project administration, Funding acquisition, Conceptualization, Writing – review & editing, Methodology.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

The authors are deeply grateful to Raphaël Janot, Da Huo, and Hélène Tonnoir, all from LRCS, for the very useful scientific discussions and insightful perspectives, while also providing the raw experimental data used in this work. A.A.F., M. P., and I.C.S. acknowledge funding from the French National Research Agency (ANR) on the context of the "France 2030 program" and the Priority Research Programs and Equipment (PEPR) "Batteries" (Project BATMAN, Reference 22-PEBA-0002). A.A.F and F.F. acknowledge financial support from the EU Horizon Europe project DigiCell (Grant Agreement no. 101135486). The authors also acknowledge Jean-Baptiste Hoock and Laurent Renault for managing the MatriCS HPC platform from Université de Picardie Jules Verne on which the simulations were performed. A.A.F acknowledges Institut Universitaire de France for the support.
## Nomenclature

$E_{eq, int}$ - Equilibrium potential for intercalation phenomenon in HC electrode

$E_{eq, fill}$ - Equilibrium potential for nanopore filling phenomenon in HC electrode

$E_{eq, Na}$ - Equilibrium potential for sodium metal electrode
$c_{s,max}$ - Maximum concentration of Na in HC

$c_{s,max, int}$ - Maximum concentration of Na in HC due to the intercalation phenomenon

$c_{s,max, fill}$ - Maximum concentration of Na in HC due to the nanopore filling phenomenon

$c_{s,0, int}$ - Minimum concentration of Na in HC due to the intercalation phenomenon

$c_{s,0, fill}$ - Minimum concentration of Na in HC due to the nanopore filling phenomenon

$\text{soc}_{0,HC}$ - Minimum electrode state of charge (Stoichiometry at 0 % discharge)

$\text{soc}_{\text{max},HC}$ - Maximum electrode state of charge (Stoichiometry at 100 % discharge)

$c_{I,0}$ - Electrolyte salt concentration (Initial value)

$c_{I,ref}$ - Reference concentration of salt in electrolyte

$t_+$ - Transport number of the salt in the electrolyte

$\varepsilon_{I, CBD}$ - Electrolyte volume fraction in Separator

$\varepsilon_{I, sep}$ - Electrolyte volume fraction in Separator

$\varepsilon_{HC,fill}$ - Intercalation surface fraction in HC domain

$\varepsilon_{HC, int}$ - Nanopore filling surface fraction in HC domain

$\rho_{HC}$ - Density of HC material

$\text{vol}_{HC}$ - Volume of the generated HC structures

$m_{HC}$ - Mass of HC in the generated structures

$I_{app}$ - Applied discharge current at the current collector boundary

$Q_{th}$ - Theoretical capacity of NaC₆ compound

$D_{HC}$ - Solid diffusion coefficient for sodium in the HC

$\sigma_I$ - Conductivity of electrolyte

$\sigma_{s,HC}$ - Conductivity of HC

$\sigma_{s,CBD}$ - Conductivity of CBD

$D_I$ - Liquid diffusion coefficient for Na⁺ in the electrolyte

b - Bruggeman exponent

$k_{\text{int}}$ - Rate constant for HC electrode (intercalation mechanism)

$k_{\text{fill}}$ - Rate constant for HC electrode (nanopore filling mechanism)

$k_{\text{Na}}$ - Rate constant for sodium metal electrode

$C_{dI}$ - Electrical Double Layer Capacitance

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jpowsour.2025.237878.

## Data availability

Data will be shared upon reasonable request.

## References

[1] D. Schäfer, K. Hankins, M. Allion, U. Krewer, F. Karcher, L. Derr, R. Schuster, J. Maibach, S. Mück, D. Kramer, R. Mönig, F. Jeschull, S. Daboss, T. Philipp, G. Neusser, J. Romer, K. Palanisamy, C. Kranz, F. Buchner, R.J. Behm, A. Ahmadian, C. Kübel, I. Mohammad, A. Samoson, R. Witter, B. Smarsly, M. Rohnke, Multiscale investigation of sodium-ion battery anodes: analytical techniques and applications, Adv. Energy Mater. (2024), https://doi.org/10.1002/aemn.202302830.

[2] T.C. Wanger, The Lithium future—resources, recycling, and the environment, Conserv. Lett. 4 (2011) 202–206, https://doi.org/10.1111/j.1755-263X.2011.00166.x.

[3] M.D. Slater, D. Kim, E. Lee, C.S. Johnson, Sodium-Ion batteries, Adv. Funct. Mater. 23 (2013) 947–958, https://doi.org/10.1002/ADFM.201200691.

[4] N. Yabuuchi, K. Kubota, M. Dahbi, S. Komaba, Research development on sodium-ion batteries, Chem. Rev. 114 (2014) 11636–11682, https://doi.org/10.1021/cr500192f. JPEG.

[5] C. Vaalma, D. Buchholz, M. Weil, S. Passerini, A cost and resource analysis of sodium-ion batteries, Nat. Rev. Mater. 3 (4 3) (2018) 1–11, https://doi.org/10.1038/natrevmats.2018.13.

[6] D.A. Stevens, J.R. Dahn, An in situ small-angle X-Ray scattering study of sodium insertion into a nanoporous carbon anode material within an operating electrochemical cell, J. Electrochem. Soc. 147 (2000) 4428, https://doi.org/10.1149/1.1394081.XML.

[7] P. Thomas, D. Billaud, Electrochemical insertion of sodium into hard carbons, Electrochim. Acta 47 (2002) 3303–3307, https://doi.org/10.1016/S0013-4686(02)00250-5.

[8] D.A. Stevens, J.R. Dahn, The mechanisms of lithium and sodium insertion in carbon materials, J. Electrochem. Soc. 148 (2001) A803, https://doi.org/10.1149/1.1379565.

[9] Z. Wang, S.M. Selbach, T. Grande, Van der Waals density functional study of the energetics of alkali metal intercalation in graphite, RSC Adv. 4 (2013) 4069–4079, https://doi.org/10.1039/C3RA47187J.

[10] H. Moriwake, A. Kuwabara, C.A.J. Fisher, Y. Ikuhara, Why is sodium-intercalated graphite unstable? RSC Adv. 7 (2017) 36550–36554, https://doi.org/10.1039/C7RA06777A.

[11] O. Lenchuk, P. Adelhelm, D. Mollenhauer, New insights into the origin of unstable sodium graphite intercalation compounds, Phys. Chem. Chem. Phys. 21 (2019) 19378–19390, https://doi.org/10.1039/C9CP03453F.

[12] M.M. Titirici, H. Alptekin, H. Au, A.C.S. Jensen, E. Olsson, M. Gokts, T.F. Headen, P. Adelhelm, Q. Cai, A.J. Drew, Sodium storage mechanism investigations through structural changes in hard carbons, ACS Appl. Energy Mater. 3 (2020) 9918–9927, https://doi.org/10.1021/acsaem.0c01614.

[13] C. Bommier, T.W. Surta, M. Dolgos, X. Ji, New mechanistic insights on Na-Ion storage in nongraphitizable carbon, Nano Lett. 15 (2015) 5888–5892, https://doi.org/10.1021/acs.nanolett.5b01969.

[14] P. Bai, Y. He, X. Zou, X. Zhao, P. Xiong, Y. Xu, Elucidation of the sodium-storage mechanism in hard carbons, Adv. Energy Mater. 8 (2018), https://doi.org/10.1002/AENM.201703217.

[15] V. Surendran, R.K. Hema, M.S.O. Hassan, V. Vijayan, M.M. Shajumon, Open or closed? Elucidating the correlation between micropore nature and sodium storage mechanisms in hard carbon, Batter Supercaps 5 (2022) e202200316, https://doi.org/10.1002/BATT.202200316.

[16] H. Tonnoir, D. Huo, R.L.S. Canevesi, V. Fierro, A. Celzard, R. Janot, Tannin-based hard carbons as high-performance anode materials for sodium-ion batteries, Mater. Today Chem. 23 (2022), https://doi.org/10.1016/j.mtcchem.2021.100614.

[17] C.S. Del Mar Rios, A. Beda, L. Simonin, C.M. Ghembeu, Hard carbon for Na-ion batteries: from synthesis to performance and storage mechanism, Na-Ion Batteries (2021) 101–146, https://doi.org/10.1002/9781119818069.CH3.

[18] D. Saurel, B. Orayech, B. Xiao, D. Carriazo, X. Li, T. Rojo, From charge storage mechanism to performance: a roadmap toward high specific energy sodium-ion batteries through carbon Anode optimization, Adv. Energy Mater. 8 (2018), https://doi.org/10.1002/aenm.201703268.

[19] B. Xiao, T. Rojo, X. Li, Hard carbon as sodium-ion battery anodes: progress and challenges, ChemSusChem 12 (2019) 133–144, https://doi.org/10.1002/cscs.201801879.

[20] C. Matei Ghimbeu, A. Beda, B. Réty, H. El Marouazi, A. Vizintin, B. Tratnik, L. Simonin, J. Michel, J. Abou-Rjeily, R. Dominko, Review: insights on hard carbon materials for sodium-ion batteries (SIBs): synthesis – properties – performance relationships, Adv. Energy Mater. 14 (2024), https://doi.org/10.1002/aenm.202303833.

[21] D.A. Stevens, J.R. Dahn, High capacity anode materials for rechargeable sodium-ion batteries, J. Electrochem. Soc. 147 (2000) 1271, https://doi.org/10.1149/1.1393348.XML.

[22] N. Sun, J. Qiu, B. Xu, Understanding of sodium storage mechanism in hard carbons: ongoing development under debate, Adv. Energy Mater. 12 (2022), https://doi.org/10.1002/aenm.202200715.

[23] Y. Cao, L. Xiao, M.L. Sushko, W. Wang, B. Schwenzer, J. Xiao, Z. Nie, L.V. Saraf, Z. Yang, J. Liu, Sodium ion insertion in hollow carbon nanowires for battery applications, Nano Lett. 12 (2012) 3783–3787, https://doi.org/10.1021/NL3016957.

[24] S. Qiu, L. Xiao, M.L. Sushko, K.S. Han, Y. Shao, M. Yan, X. Liang, L. Mai, J. Feng, Y. Cao, X. Ai, H. Yang, J. Liu, Manipulating adsorption-insertion mechanisms in nanostructured carbon materials for high-efficiency sodium ion storage, Adv. Energy Mater. 7 (2017) 1700403, https://doi.org/10.1002/AENM.201700403.

[25] B. Zhang, C.M. Ghimbeu, C. Laberty, C. Vix-Guterl, J.M. Tarascon, Correlation between microstructure and Na storage behavior in hard carbon, Adv. Energy Mater. 6 (2016), https://doi.org/10.1002/aenm.201501588.

[26] Y. Li, Y.-S. Hu, M.-M. Titirici, L. Chen, X.Y. Huang M Li, Y. Hu, L.Q. Chen, X. J. Huang, Y.M. Li, M. Titirici, Hard carbon microtubes made from renewable cotton as high-performance anode material for sodium-ion batteries, Adv. Energy Mater. 6 (2016) 1600659, https://doi.org/10.1002/AENM.201600659.

[27] J.M. Stratford, P.K. Allan, O. Pecher, P.A. Chater, C.P. Grey, Mechanistic insights into sodium storage in hard carbon anodes using local structure probes, Chem. Commun. 52 (2016) 12430–12433, https://doi.org/10.1039/C6CC06990H.

[28] H. Au, H. Alptekin, A.C.S. Jensen, E. Olsson, C.A. O'Keefe, T. Smith, M. Crespo-Ribadeneira, T.F. Headen, C.P. Grey, Q. Cai, A.J. Drew, M.M. Titirici, A revised mechanistic model for sodium insertion in hard carbons, Energy Environ. Sci. 13 (2020) 3469–3479, https://doi.org/10.1039/d0ee01363c.

[29] J.F. Troncoso, F.M. Zanotto, D.E. Galvez-Aranda, D. Zapata Dominguez, L. Denisart, A.A. Franco, The ARTISTIC battery manufacturing digitalization initiative: from fundamental research to industrialization, Batter Supercaps 8 (2025) e202400385, https://doi.org/10.1002/BATT.202400385.
[30] F. Brosa Planella, W. Ai, A.M. Boyce, A. Ghosh, I. Korotkin, S. Sahu, V. Sulzer, R. Timms, T.G. Tranter, M. Zyskin, S.J. Cooper, J.S. Edge, J.M. Foster, M. Marinescu, B. Wu, G. Richardson, A continuum of physics-based lithium-ion battery models reviewed, Prog. Energy 4 (2022) 042003, https://doi.org/10.1088/2516-1083/AC7D31.

[31] A.C. Ngandjong, A. Rucci, M. Maiza, G. Shukla, J. Vazquez-Arenas, A.A. Franco, Multiscale simulation platform linking lithium ion battery electrode fabrication process with performance at the cell level, J. Phys. Chem. Lett. 8 (2017) 5966-5972, https://doi.org/10.1021/acs.jpclett.7b02647.

[32] T. Lombardo, A.C. Ngandjong, A. Belhcen, A.A. Franco, Carbon-Binder migration: a three-dimensional drying model for lithium-ion battery electrodes, Energy Storage Mater. 43 (2021) 337-347, https://doi.org/10.1016/J.ENSM.2021.09.015.

[33] J. Xu, A.C. Ngandjong, C. Liu, F.M. Zanotto, O. Arcelus, A. Demortière, A. A. Franco, Lithium ion battery electrode manufacturing model accounting for 3D realistic shapes of active material particles, J. Power Sources 554 (2023) 232294, https://doi.org/10.1016/J.JPOWSOUR.2022.232294.

[34] M. Chouchane, A. Rucci, T. Lombardo, A.C. Ngandjong, A.A. Franco, Lithium ion battery electrodes predicted from manufacturing simulations: assessing the impact of the carbon-binder spatial location on the electrochemical performance, J. Power Sources 444 (2019), https://doi.org/10.1016/j.jpowsour.2019.227285.

[35] M. Chouchane, E.N. Primo, A.A. Franco, Mesoscale effects in the extraction of the solid-state lithium diffusion coefficient values of battery active materials: physical insights from 3D modeling, J. Phys. Chem. Lett. 11 (2020) 2775-2780, https://doi.org/10.1021/acs.jpclett.0c00517.

[36] A.C. Ngandjong, T. Lombardo, E.N. Primo, M. Chouchane, A. Shodiev, O. Arcelus, A.A. Franco, Investigating electrode calendering and its impact on electrochemical performance by means of a new discrete element method model: towards a digital twin of Li-Ion battery manufacturing, J. Power Sources 485 (2021) 229320, https://doi.org/10.1016/J.JPOWSOUR.2020.229320.

[37] C. Liu, T. Lombardo, J. Xu, A.C. Ngandjong, A.A. Franco, An experimentally-validated 3D electrochemical model revealing electrode manufacturing parameters' effects on battery performance, Energy Storage Mater. 54 (2023) 156-163, https://doi.org/10.1016/j.ensm.2022.10.035.

[38] C. Liu, O. Arcelus, T. Lombardo, H. Oularbi, A.A. Franco, Towards a 3D-resolved model of Si/Graphite composite electrodes from manufacturing simulations, J. Power Sources 512 (2021), https://doi.org/10.1016/j.jpowsour.2021.230486.

[39] A. Marcato, J.E. Santos, C. Liu, G. Boccardo, D. Marchisio, A.A. Franco, Modeling the 4D discharge of lithium-ion batteries with a multiscale time-dependent deep learning framework, Energy Storage Mater. 63 (2023), https://doi.org/10.1016/j.ensm.2023.102927.

[40] K. Chayambuka, M. Jiang, G. Mulder, D.L. Danilov, P.H.L. Notten, Physics-based modeling of sodium-ion batteries part I: experimental parameter determination, Electrochim. Acta 404 (2022), https://doi.org/10.1016/j.electacta.2021.139726.

[41] K. Chayambuka, G. Mulder, D.L. Danilov, P.H.L. Notten, Physics-based modeling of sodium-ion batteries part II. Model and validation, Electrochim. Acta 404 (2022), https://doi.org/10.1016/j.electacta.2021.139764.

[42] V.K. Garapati, N.N. Dingari, M. Mynam, B. Rai, Physics-based reduced order model for sodium-ion batteries, J. Electrochem. Soc. 170 (2023) 010517, https://doi.org/10.1149/1945-7111/ACB01B.

[43] M. Neumann, T. Philipp, M. Häringer, G. Neusser, J.R. Binder, C. Kranz, Stochastic 3D modeling of nanostructured NVP/C active material particles for sodium-ion batteries. https://doi.org/10.1002/batt.202300409, 2023.

[44] Math2Market GmbH, Geodict Simulation Software Release 2024, (n.d.). doi.org/10.30423/release.geodict2024 (accessed May 3, 2024).

[45] S. Pinilla, F.M. Zanotto, D.Z. Dominguez, T. García, A.A. Franco, Carbon-binder-domain porosity extraction through lithium-ion battery electrode impedance data, Energy Storage Mater. 74 (2025) 103818, https://doi.org/10.1016/J.ENSM.2024.103818.

[46] M. Chouchane, A. Rucci, A.A. Franco, A versatile and efficient voxelization-based meshing algorithm of multiple phases, ACS Omega 4 (2019) 11141-11144, https://doi.org/10.1021/acsomega.9b01279.

[47] P. Albertus, J. Christensen, J. Newman, Experiments on and modeling of positive electrodes with multiple active materials for lithium-ion batteries, J. Electrochem. Soc. 156 (2009) A606, https://doi.org/10.1149/1.3129656/XML.

[48] J.S. Newman, C.W. Tobias, Theoretical analysis of current distribution in porous electrodes, J. Electrochem. Soc. 109 (1962) 1183, https://doi.org/10.1149/1.2425269.

[49] M. Doyle, T.F. Fuller, J. Newman, Modeling of galvanostatic charge and discharge of the Lithium/polymer/insertion cell, J. Electrochem. Soc. 140 (1993) 1526-1533, https://doi.org/10.1149/1.2221597/XML.

[50] T.F. Fuller, M. Doyle, J. Newman, Simulation and optimization of the dual lithium ion insertion cell, J. Electrochem. Soc. 141 (1994) 1-10, https://doi.org/10.1149/1.2054684/XML.

[51] R. Gond, H.D. Asfaw, O. Hosseinaei, K. Edström, R. Younesi, A.J. Naylor, A lignosulfonate binder for hard carbon anodes in sodium-ion batteries: a comparative study, ACS Sustain. Chem. Eng. 9 (2021) 12708-12717. https://pubs.acs.org/doi/10.1021/acssuschemeng.1c05263.

[52] A. Beda, F. Rabuel, M. Morcrette, S. Knopf, P.L. Taberna, P. Simon, C. Matei Ghimbeu, Hard carbon key properties allow for the achievement of high Coulombic efficiency and high volumetric capacity in Na-ion batteries, J. Mater. Chem. A Mater. 9 (2021) 1743-1758, https://doi.org/10.1039/d0ta07687b.

[53] D.E. Galvez-Aranda, F. Fernandez, A.A. Franco, Physics-Assisted Machine Learning for the Simulation of the Slurry Drying in the Manufacturing Process of Battery Electrodes: A Hybrid Time-Dependent VGG16-DEM Model, ACS Appl Mater Interfaces 17 (2025) 32150-32162, https://doi.org/10.1021/acsami.4c23103.

[54] Y. Aihara, G.B. Appetecchi, M. Armand, S. Passerini, C.A. Vincent, Bruno Scrosati (1937-2024), Nat. Nanotechnol. 20 (2025) 178, https://doi.org/10.1038/s41565-025-01861-z.