Physics-based modeling of sodium-ion batteries part II. Model and validation

Sodium-ion batteries (SIBs) have recently been proclaimed as the frontrunner ’post lithium’ energy storage technology. This is because SIBs share similar performance metrics with lithium-ion batteries, and sodium is 10 0 0 times more abundant than lithium. In order to understand the electrochemical characteristics of SIBs and improve present-day designs, physics-based models are necessary. Herein, a physicsbased, pseudo-two-dimensional (P2D) model is introduced for SIBs for the first time. The P2D SIB model is based on N a 3 V 2 ( P O 4 ) 2 F 3 (NVPF) and hard carbon (HC) as positive and negative electrodes, respectively. Charge transfer in the NVPF and HC electrodes is described by concentration-dependent diffusion coefficients and kinetic rate constants. Parametrization of the model is based on experimental data and genetic algorithm optimization. It is shown that the model is highly accurate in predicting the discharge profiles of full cell HC//NVPF SIBs. In addition, internal battery states, such as the individual electrode potentials and concentrations, can be obtained from the model at applied currents. Several key challenges in both electrodes and the electrolyte are herein unraveled, and useful design considerations to improve the performance of SIBs are highlighted. © 2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Several research groups and startups have recently shown great interest in developing sodium-ion batteries (SIBs) [1] . SIBs are, at present, regarded as the most promising complementary technology to the ubiquitous lithium-ion batteries (LIBs) [2] . This is because sodium is 10 0 0 times more abundant than lithium on the earth's crust [ 3 , 4 ]. In addition, both battery types share similarities in performance, physical structure, and manufacturing infrastructure. As a result, the scientific knowledge and know-how accumulated in developing LIBs have been transferred in the past decade to accelerate the commercialization effort s of SIBs [5] . The most encouraging outcome of this recent drive is the availability of a large repository of cathode material choices for SIB applications, all of which are based on earth-abundant elements [3] . Therefore, monolithic slab with straight pores and high conductivity. Newman et al . [ 12 , 13 ] later improved this simplified model structure by treating electrode particles as a distinct phase in intimate contact with the electrolyte. Using the principles of homogenization, the particles were treated as a macro-homogeneous phase. Newman's model, therefore, provided the basic framework for the rigorous treatment of charge transport in discrete and conceptually spherical electrode particles. This multi-phase, multi-scale coupling is often referred to as the pseudo-two-dimensional (P2D) model structure because of the 1D representation of the electrode thickness and an additional, pseudo-dimension, representing the spherical radius of active particles at different electrode positions.
Various P2D models have been applied to different battery chemistries, such as lead-acid and nickel-metal hydride [ 14 , 15 ]. Although P2D models are widely accepted and demonstrate unparalleled accuracy and reliability, there remain practical challenges to parametrize new chemistries and integrate the models in BMS microcontrollers [16] . This is because the models are based on systems of coupled partial differential equations (PDEs), which are computationally expensive and potentially non-convergent during execution [17] . This fundamental challenge has propelled a growing trend of using reduced-order models such as single-particle models [18][19][20] , equivalent circuit models [ 21 , 22 ], and data-driven semi-empirical models [ 23 , 24 ].
Details of the reduced-order battery models have been published in thematic reviews for the interested reader [ 25 , 26 ]. Nevertheless, as models become increasingly simplified, the danger is obtaining parameters that are detached from electrochemistry and physics. In the end, the simplified models cannot be reliably used as predicting tools to improve cell design because the underlying parameters lack physical meaning and are not valid outside the conditions of model parametrization. Therefore, the development of physics-based models remains an important undertaking to understand internal battery dynamics and provide a link with experimentally derived parameters. This development should carefully consider all relevant electrochemical processes involved in the given battery chemistry to achieve an accurate physical model.
Herein, a physics based, P2D model of a SIB full cell is presented for the first time to understand and improve the design of this emerging battery chemistry. The experiments used to derive parameters of a SIB based on hard carbon (HC) as the anode/negative electrode and N a 3 V 2 ( P O 4 ) 2 F 3 (NVPF) as cathode/positive electrode are described in a separate, preceding paper [27] . The electrolyte is composed of 1 kmol m −3 NaP F 6 salt dissolved in equal weight mixtures of ethylene carbonate (EC) and propylene carbonate (PC), E C 0 . 5 : P C 0 . 5 (w/w) solvent. Based on the experimental evidence, the SIB electrode and electrolyte parameters are concentration-dependent. As a result, concentrationdependent diffusion coefficients, kinetic rate constants, and conductivity are introduced to the P2D model.
The full cell SIB model is scripted in MATLAB, which disposes of a global optimization toolbox to determine parameters that could not be experimentally deduced. Using the genetic algorithm for the optimization procedure, the SIB model is validated by comparing the simulation results with experimental voltage data for the positive and negative electrode potentials. Analyses of the model results reveal mass transport limitations in the 1 kmol m −3 NaP F 6 E C 0 . 5 : P C 0 . 5 (w/w) electrolyte and in the HC and NVPF active particles. This can further guide the design of SIB systems, which are expected to operate at high-power demanding applications. Fig. 1 shows a detailed layout of a three-electrode SIB setup composed of an HC negative electrode, a separator, an NVPF positive electrode, and a metallic sodium reference electrode (Na-RE).

Description of the system components
These are assembled in a PAT-Cell (EL-Cell GmbH). HC and NVPF electrodes, in this case, act as two working electrodes, and the Na-RE acts as reference electrode of the first kind. The Na-RE is carefully positioned between the two working electrodes for accurate potential measurements while being electronically isolated from either electrode by the separator.
The thickness of the negative electrode, the separator, and the positive electrode are δ n = 64 , δ s = 25 and δ p = 68 [μm], respectively. Based on experimental scanning electron micrographs, the average radii of the HC and NVPF electrode particles are R n = 3 . 48 and R p = 0 . 59 [μm], respectively [27] . R n is, therefore, 6 times larger than R p on average. Both electrodes additionally contain graphitic conductive additives in order to enhance the electrical conductivity of the composite electrodes. These are represented in Fig. 1 by the black spheres. More conductive additives are needed in the NVPF positive electrode because of the low electronic conductivity of the NVPF material [ 28 , 29 ]. Fig. 1 also shows the cable connections from the potentiostat to the three electrodes for automated cycling and cell voltage measurements. A current I app [A] is specified in the potentiostat program to either charge or discharge the SIB. The Na-RE is connected to a high impedance lead of the potentiostat, which ensures a very low current passes through the RE. The three-electrode setup, therefore, allows for the accurate determination of individual electrode potentials v s. Na-RE. As a result, positive electrode potential ( V p ), and the negative electrode potential ( V n ), can be deconvoluted from the full cell voltage ( V bat ) at different values of I app .
From a modeling perspective, this knowledge of the individual electrode potentials is important for two reasons. First, the parameters for both electrodes can be independently optimized instead of relying only on V bat , which is a combination of the two electrode potentials. In this way, parameters for the individual electrodes can be independently optimized. Second, the number of simultaneously optimized parameters is reduced, which increases optimization speed and model fidelity. Care, however, must be exercised on the position of the reference electrode to minimize overpotentials and voltage crosstalk between the anode and the cathode.

Model description
In the isothermal P2D model described herein, the active particles are considered spherical. Another assumption is that the particle sizes are homogeneous and represented by the average particle radius. The model variables include the Na-concentration ( c θ ,m ), the potential ( ϕ θ ,m ), and current ( i θ ,m ). The subscript θ symbolizes the phase of the variable, which can either be the solid phase ( θ = 1 ) or the liquid/electrolyte phase ( θ = 2 ), subscript m symbolizes the domain inside the battery stack, which can either be the negative electrode ( m = n ), the positive electrode ( m = p) or the separator ( m = s ).

Mass transport in electrode particles
Fick's second law expresses the time-dependent radial transport of intercalated N a + inside the electrode active particles where j m is the interfacial flux of species [ mol m −2 s −1 ]. The boundary conditions in Eqs. (2) and 3 express the surface reaction flux and spherical symmetry, respectively. An initial condition is further required for the particle phase concentrations, which is defined as where c 0 1 ,m is the initial N a + concentration inside electrode particles. c 0 1 ,m depends on the initial state-of-charge (SOC) of the active materials.
In the case of N a + intercalation with constant D 1 ,m and R m , Eqs. (1) to 4 can be solved by fast analytical methods [ 30 , 31 ]. R m changes in intercalation active materials are generally very low. In the case of NVPF and HC electrode materials, unit cell volume changes of approximately 2% have been reported [32] . However, experimental and modeling galvanostatic intermittent titra-tion technique (GITT) results showed that D 1 ,m is strongly dependent on c 1 ,m [ 27 , 33 , 34 ]. In other studies, D 1 ,p has also been shown to vary by two orders of magnitude [33] . For this reason, a concentration-dependent D 1 ,m is used for the NVPF//HC SIB. As a result, the numerical method of the hybrid backward Euler control volume method (HBECV) is applied instead of the analytical methods. The HBECV obtains fast and accurate solutions and has been reported before elsewhere [35] .

Electrode kinetics model
At the particle surface, the electrode kinetics can be described by the Butler-Volmer expression [36] where j 0 ,m is the exchange flux density of Na-ions across the elec- The charge transfer overpotential, η ct m can then be expressed as where ϕ 1 ,m , ϕ 2 ,m and U m are the electrode, electrolyte, and equilibrium potentials (EMF), respectively [V]. The EMF potentials for both the HC and NVPF electrodes were experimentally determined as described in an accompanying publication [27] .

Current distribution
Throughout the separator and porous electrode regions, the current is distributed between the electronic current density ( i 1 ,m ) and the ionic current density in the solid and electrolyte phases ( i 2 ,m ). Both i 1 ,m and i 2 ,m are related to the total applied current density i app = I app / A cc as At the current collector boundaries of the porous electrodes The boundary condition in Eq. (10) specifies that only an electronic current is present at the current collectors. Applying Eq. (8) , therefore, results in In the porous electrodes where m = { n, p } ), i 1 , m can be modeled by Ohm's law where σ e f f m is the effective electronic conductivity in the porous electrode [ −1 m −1 ]. In addition, the derivative of i 2 ,m is proportional to j m in Eq. (5) . This relation is expressed as where a m is the electrode active surface area, volume ratio [ m −1 ]. a m is calculated from the particle radius and the electrode porosity as where ε el m and ε f il l er m are the electrolyte and additive filler volume fractions, respectively. ε f il l er m includes the binder and conductive filler additives.

Electrolyte potential and mass distribution
The dilute solution theory governs the electrolyte potential distribution in the liquid electrolytes. This theory is based on the Nernst-Planck equation, a classical description of the transport of charged ionic species in electrolyte media [37] . The dilute solution theory essentially considers binary interactions between ionic species and the solvent and neglects ion-ion pairing effects [38] .
Because the NaP F 6 E C 0 . 5 : P C 0 . 5 (w/w) electrolyte does not show extensive ion-pairing effects in the concentration range of 0 to 2 mol k g −1 as earlier determined [39] , electrolyte potential can thus be modeled by the expression [40] .
where κ eff m is effective ionic conductivity in the electrolyte in cell domain m [ −1 m −1 ], and t + is the cationic transference number [-], defined as the fraction of i 2 ,s due to cationic migration in the absence of diffusion and convection forces. The definition of ϕ 2 , m results in the multiplier ( 1 − 2 t + ) , which is different from the multiplier 2( 1 − t + ) used in other works [41] .
A Dirichlet boundary condition is thus defined for the electrolyte potential at the anode Eq. (16) sets ϕ 2 ,n as the reference potential for all potential difference measurements. Another option is to set ϕ 1 ,n | x =0 = 0 , which is equivalent to grounding the anode. In either case, the position and choice of the reference potential do not affect the overpotentials and the overall cell voltage [12] . However, the convenience of the boundary condition in Eq. (16) is that electrode potentials ϕ 1 ,n and ϕ 1 ,p have values similar to those obtained experimentally using a reference electrode of the first kind. This property makes the model validation based on individual electrode potentials using the Na-RE straightforward.
The electrolyte concentration mass balance in the porous electrode region is expressed as and in the separator region as where  [40] .
Two symmetrical Neumann boundary conditions are needed to resolve the concentration profile. These are expressed at the negative electrode/current collector boundary ( x = 0 ) and at the positive electrode/current collector boundary ( x = L ) as Boundary conditions Eqs. (19) and 20 state that there is no flux of ionic species at the current collector/electrode interface. These are identical in the case of a full cell battery with two porous electrodes.
In addition, an initial condition is needed for c 2 , which is defined as where c 0 2 is the initial concentration at rest/equilibrium, equal to 1 kmol m −3 . Based on experimental conductivity studies, the conductivity of NaP F 6 in E C 0 . 5 : P C 0 . 5 (w/w) electrolyte is also highest around this concentration [39] .

Relation between bulk transport properties and porous electrode properties
Bulk properties such as the diffusion coefficients and the electrode conductivity need to be related to the volume fraction of the bulk material in the porous electrode. The Bruggeman correlation [42] is herein used to define the effective electrolyte transport properties of conductivity and diffusion coefficient where κ and D 2 represent the bulk electrolyte conductivity and diffusion coefficient, separately determined by conductivity experiments and the advanced electrolyte model (AEM) version 2.19.1 [39] . The AEM is a statistical-mechanics-based simulation model for determining electrolyte properties [43][44][45] . A brief discussion on the Bruggeman exponent value of 1 . 5 in Eqs. (22) and 23 is necessary. This exponent includes a hidden factor of 1, which correlates the bulk electrolyte concentration to the porous media concentration where c 2 ,m is the electrolyte concentration, which includes the volume of the solid phases in porous media [ mol m −3 ]. It is equally valid to use the Bruggeman exponent of 0.5, in which case, the correlated c 2 ,m is used in the model equations. The advantage, however, of using the uncorrelated concentration c 2 is that electrolyte concentrations in different cell domains can be treated as continuous functions. The Bruggeman correlation was nevertheless not used for σ e f f m since there was no prior knowledge of the bulk electrode conductivity and the electrode porosity was not fixed. Therefore, the optimization of the electrode conductivity was performed for the effective property without including a correlation.

Modeling interfaces
The first-order and second-order PDEs in the electrolyte were discretized by a finite difference method (FDM). Because of porosity differences in the battery domains, the transport properties can change across the interface. Thus the FDM mass balance may be inaccurate across the interface. To overcome this, some authors have proposed an effective interfacial diffusion coefficient [12] . In this work, however, a control volume method (CVM) was used to discretize the interfacial boundary of the separator and porous electrode. Detailed descriptions of the CVM have been published [46] . For the interested reader, Botte et al . [ 47 , 48 ] compared the FDM and the CVM in battery modeling applications. In general, the CVM results in perfect mass conservation across the interface. To obtain the concentration at the interface node, an imaginary control volume is defined across the interface node, and fluxes across the faces of the imaginary volume are calculated, assuming a linear concentration profile. Fig. 2 illustrates the elements of the CVM herein used to derive the interfacial concentration.
The interfacial concentration at the negative electrode/separator interface can therefore be determined as where x n and x s are the node spacings in the negative electrode and separator [m], respectively, while b, e, and w are representing the interface, the outermost east boundary, and the west outer boundary of the imaginary control volume, respectively. A similar mass balance in a control volume expression can be derived at the positive electrode/separator interface.

Battery voltage
Finally, a Dirichlet boundary condition for the electrode potential on each electrode must be considered [49] . Suppose that the electrode potentials at the left boundary of each porous electrode ( ϕ 1 ,n | x =0 and ϕ 1 ,p | x = δn + δs ) are set to some arbitrary values. These values influence the ionic current densities at the opposing electrode end at i 1 ,n | x = δn and i 2 ,p | x = L , respectively. To ascertain if the imposed boundary conditions are correct, Eq. (9) and the boundary conditions of Eqs. (10) and 11 must be satisfied. The functions to be solved at the negative electrode and the positive electrode can therefore be expressed by and where I 1 ,n represents electronic current at the interface between the anode and separator, while I 2 ,p represent the ionic current at the interface with the cathode current collector. These currents are considered as generated by coupled PDEs for the negative and positive electrodes, respectively.
Eqs. (26) and 27 can be solved by various iterative root-finding methods, in which an approximate value of the electrode potential is supplied as an initial guess. In the literature, such a method is also called 'shooting' [ 50 , 51 ]. In this model, a combination of the dichotomy method and the secant method is applied to optimize the solution convergence [52] . Because the system of equations in porous battery electrodes is nonlinear, it is possible to obtain intermediate solutions of extreme magnitude for minor deviations in the initial guess. A robust root finding method is therefore needed in the first iterations. To our knowledge, the dichotomy method is the only method capable of finding the root in cases where solutions on the right-hand side of the equation can have infinite values. The dichotomy method uses two guess values of the electrode potential, an overestimate and an underestimate, and then reduces the estimated range until the solution is found within the error tolerance. In this case, the tolerance is set at 0.01% of i app . However, the dichotomy method requires many iterations and is therefore slow to converge to the root. For faster root finding, the solution method switches to the secant method after 5 iterations of the dichotomy method, when the right-hand side value is finite.
Once the solution is found within the error tolerance, the full cell battery voltage can then be determined by (28) where R contact,n and R contact,p [ m 2 ] are the negative and positive electrode current collector/porous electrode contact specific resistances, respectively. The importance of the contact resistance is revealed at high currents [ 13 , 53 ]. In this work, the contact resistances of individual electrodes were calculated from the Na-RE measurements [27] .
It is worth highlighting that the validation of the model vs. experimental data is done based on the individual electrode potentials and not the full cell voltage, V bat . This scheme is applied because, from Eq. (28) , it is possible to have wrong values of both ϕ 1 ,p and ϕ 1 ,n and yet still manage to have a correct V bat . The Na-RE electrode deconvolutes the individual electrode potentials from V bat and thereby allows model validation on two separate electrodes.
The electrode potentials used for validation of the model vs. experimental data can thus be expressed as where V p and V n [V] are the positive and negative electrode potentials, respectively.

Parameter identification and optimization
For the developed model to provide physically meaningful results, the model parameters should be inferred from an extensive experimental data set. Because of the minimal assumptions in P2D models, experimentally derived parameters should ideally result in a fitting model. However, due to the disparity in the definitions of key parameters, such as the diffusion coefficients and transference number between experimental and modeling techniques, this is seldom the case. Another challenge is that experimentally derived parameters are technique-dependent. There is, therefore, a great need to bridge the gap between model and experimental parameters. Nevertheless, experimental parameters are an ideal starting point and provide insight into the order of magnitude of the model parameters.
In this work, the geometric parameters of the SIBs were determined from experiments and used without modification [27] . These include the thicknesses of the electrodes and separator as well as the particle radii of the positive and negative electrodes. Fig. 3 shows the cell dimensions of the HC//HVPF SIB, concluded from the experiments.
The negative electrode is therefore defined between 0 ≤ x ≤ 64 μm, the separator in the 64 ≤ x ≤ 89 μm region, while the negative electrode is defined in the 89 ≤ x ≤ 157 μm region. The HC and NVPF electrode particles are separately modeled in a homogenous P2D domain, in which the r m -and x -axis represent the particle radii and particle positions in the porous electrode, respectively. Having fixed cell dimensions allows the mesh of the cell components to be defined in the model. Changes in cell dimensions during optimization require a new mesh to be defined, which may inadvertently affect parameters with length scale, such as the conductivity and diffusion coefficients. This problem can also lead to numerical instabilities and inconsistent model results. It is therefore advised to maintain constant cell dimensions unless experimental evidence proves otherwise.
The EMF of the NVPF and HC electrodes were also determined from experimental data based on the slow discharge rate experiments [27] . The particle-phase diffusion coefficients, D 1 ,m and kinetic rate constants, k m were determined from a combination of experimental GITT and P2D modeling approach (P2D GITT model) [34] . A half-cell P2D GITT model was used in combination with half-cell GITT experimental data to determine the HC and NVPF transport parameters at different electrode SOC. As a result, concentration-dependent D 1 ,n ( c 1 ,n ) , D 1 ,p ( c 1 ,p ) , k n ( c 1 ,n ) and k p ( c 1 ,p ) were thus obtained. On the other hand, electrolyte properties of D 2 ( c 2 ) , κ ( c 2 ) and t + were determined from the AEM modeling and experiments [39] .
Appendix A lists the parameters used in the model. Constant value parameters are listed in Table A1 , while concentrationdependent parameters for the HC negative electrode, NVPF positive electrode and NaP F 6 E C 0 . 5 : P C 0 . 5 (w/w) electrolyte are shown in Figs. A1 , A2 , and A3 , respectively. These figures also show a comparison between the experimental and model optimized parameters.
For optimizing the unknown model parameters, the root-meansquare error between the model and the experimental results was defined as the objective function. The MATLAB genetic algorithm (GA) was used to obtain the error minimum for multiple discharge curves at different rates [54] . The GA is necessary due to the nonlinearity of the P2D model equations and parameter identification complexity. Because of the availability of two experimentally determined electrode potentials for each discharge curve, a two-step optimization procedure was therefore followed, in which parameters for the positive electrode were optimized in the first step, separate from the parameters of the negative electrode, which were then optimized in the second step. The concentration-dependent parameters were optimized by means of a scaling factor. This strategy resulted in improved optimization results at high rates. Fig. 4 shows the measured (symbols) and simulated (solid lines) SIB discharge voltage profiles as a function of the transferred charge during discharge. The model results are obtained using a single set of optimized parameters for all discharge rates. Discharge current densities of 1, 5, 10, and 12 A m 2 were applied, corresponding to 0.1, 0.6, 1.2, and 1.4 C-rate, respectively. Fig. 4 a shows the measured and simulated voltage profiles of an HC//NVPF, full cell SIB. The fully-charged cell voltage starts at 4.2 V and terminates at the cutoff voltage of 2 V. Note that the practically recommended cutoff voltage is 2.5 V [55] . As the current increases, the battery voltage and the maximum transferred charge decrease. This is because of an increase in mass transport and charge transfer overpotentials at higher currents. An accurate physics-based model is, therefore, the only way to account for the various kinetic and mass transport effects at different rates. Fig. 4 b and c shows the measured and simulated voltage profiles of the NVPF positive electrode ( V p v s. Na-RE) and the HC negative electrode ( V n v s. Na-RE), respectively. At the different discharge rates, V p varies between 4 . 3 − 3 . 4 V v s. Na-RE while V n varies between 0.1-1.5 V v s. Na-RE. Therefore, based on the potential range and the current dependence of the voltage profiles, both electrodes contribute significantly to the overpotential losses in the full cell and consequently to the capacity losses at high currents. In both cases, however, the P2D model is in good agreement with the experimental voltage profiles of V p and V n v s. Na-RE. Table 1 shows the percentage error in V bat and the mean absolute errors in V bat , V p and V n at different rates. The largest percentage error in V bat is 1.47%, corresponding to 48.1 mV in absolute error terms. Therefore, the model results match the experimental full cell voltage and the individual electrode potentials at different discharge rates. The accurate P2D model herein presented goes beyond the current density and terminal voltage data by providing additional information on internal battery states. In the subsequent figures, the Na + electrolyte concentration, the Na + concentration in the active particles, and the ionic current distribution are compared for the applied current of 1 and 12 A m −2 , to investigate how different discharge rates influence the battery performance. c 2 increases as a function of time due to N a + deintercalation in the HC negative electrode, while in the positive electrode region ( 89 ≤ x ≤ 157 μm), c 2 decreases due to N a + intercalation in NVPF. This ionic transport in the electrolyte is driven by migration and diffusion mechanisms, as expressed in Eq. (15) . Note that due to electroneutrality condition, at any given time, the average of c 2 remains constant and equal to the equilibrium and initial concentration of 1 kmol m −3 . Although the results in Fig. 5 a do not appreciably deviate from1 kmol m −3 , these slow-discharge rate profiles, however, show dynamic waves compared to the fast-discharge rate profiles shown in Fig. 5 b (see Fig. S1 in the Supplementary Information at higher  magnification). On the other hand, the results in Fig. 5  . This situation can be mitigated by optimizing the electrolyte conductivity, electrode porosity, and coating thickness [56][57][58] . The optimization objective here is to reduce the electrolyte mass transport limitations, which induce high overpotentials at high current densities.  ( Fig. 6 a and c), while the positive electrode concentrations, c 1 ,p and c s 1 ,p increase ( Fig. 6 b and d). Fig. 6 a and b show the simulated results of the intercalated concentration of N a + ions during a 1 A m −2 discharge, in the negative and positive electrode active particles, respectively. For the slow-discharge rate operation, c 1 ,m and c s 1 ,m profiles are shown to evolve uniformly along x . In addition, the profiles remain very close at all times. This behavior indicates that slight concentration gradients develop in the electrode active particles during the slow discharge rate. From a modeling perspective, such concentration profiles can be simulated quite accurately using computationally efficient analytical methods [30] . Fig. 6 c and d show the simulated results during 12 A m −2 discharge in the negative and positive electrodes, respectively. In contrast to the results shown in Fig. 6 a and b, the fast-discharge exhibits non-uniform c 1 ,m and c s 1 ,m profiles along x . This is most apparent in the positive electrode ( Fig. 6 d), where the active particles close to the separator receive 58% higher average concentration at compared to particles at the current collector (compare c 1 ,p at x = 89 and x = 157 μm in the dashed curves of Fig. 6 d). These profiles can only be obtained accurately using numerical methods because of the concentration dependence of D 1 ,n and D 1 ,p .

Results and discussion
An analysis of the difference c s 1 ,m −c 1 ,m is shown in Supplementary Information Fig. S3 as a function of discharge time.  −c 1 ,p ). A small difference indicates negligible concentration gradients in the electrode active particles. Such behavior is desirable because high concentration gradients induce mechanical strain in the particles [59] . At slow discharge rates ( Fig. S3 a and b), c s 1 ,m −c 1 ,m is not large. The maximum differences are about 0.8 and 0.25 kmol m −3 in the negative and positive electrodes, respectively. At fast discharge rates, however ( Fig. S3 b and c), c s 1 ,m −c 1 ,m is significantly large. The maximum differences are about 3.5 and 2 kmol m −3 in the negative and positive electrodes, respectively. In addition, compared to the negative electrode profile ( Fig. S3 c), the positive electrode profile monotonically increases with time ( Fig. S3 d). This increase indicates diffusion transport limitations in the NVPF particles as a result of the comparatively low D 1 ,p of the material. Such behavior is also consistent with diffusion length calculations, wherein the diffusion times for HC and NVPF were estimated to be 1 and 2 h, respectively [27] . Supplementary Information Fig. S4 shows profiles of the interfacial flux, j m at discharge rates of 1 ( Fig. S4 a and b) and 12 A m −2 ( Fig. S4 c and d). In order to analyze the distribution of the flux, the profiles are averaged over quarterly intervals of t max . It can be observed that the initial stages (red lines) are characterized by high j m rates at the electrode/separator boundary. Another observation is that profiles in the negative electrode at 1 and 12 A m −2 , Fig. S4 a and c respectively are nearly identical and scaled versions of each other. This similarity indicates that the distribution of j m in the negative electrode is not altered by the increase in the discharge rate. In contrast, profiles in the positive electrode ( Fig. S4 b and d) show considerable differences, especially in the intermediate periods 0 . 25 t max < t ≤ 0 . 5 t max (blue line) and 0 . 5 t max < t ≤ 0 . 75 t max (green line). This indicates that high discharge rates are influencing the distribution of j m in the positive electrode. Fig. 7 shows 2D simulation results of the intercalated N a + concentration in the negative ( Fig. 7 a-d) and positive ( Fig. 7 e-h) electrodes during a 1 A m −2 discharge rate. For an illustration of the relationship between 1D and 2D coordinates, refer to Fig. 3 . The 2D results in Fig. 7 show in more detail the 1D profiles shown in Fig. 6 a and b. The 2D concentrations are, however, expressed as SOC, which is defined for the negative and positive electrodes as The SOC scale is convenient for a side-by-side comparison of the 2D concentration profiles in two battery electrodes because the SOC is scaled between 0 and 1 or 0 and 100%.
At low discharge rates, the SOC is uniformly distributed within the active particles (along the r m -axis) and for particles located at different positions in the electrodes (along the x -axis). In addition, toward the end of discharge ( Fig. 7 d and h, t = 80% t max ), the SOC is low and uniformly distributed in both electrodes. This behavior signifies that the intercalated N a + is optimally utilized and that the maximum extractable capacity is attained. The electrode thickness can also be safely increased without harming the discharge performance of the cell at this discharge rate. Fig. 8 shows the 2D SOC profiles in particles of the HC ( Fig. 8 ad) and NVPF ( Fig. 8 e-h) electrodes at a discharge rate of 12 A m −2 . The results in Fig. 8 further elaborate the 1D profiles shown in Fig. 6 c and d. In contrast to the uniform concentration profiles observed in Fig. 7 , the fast discharge rate reveals non-uniform SOC distribution in both the x -and r m -axis. Fig. 8 a-d show the evolution of SOC profiles in the HC negative electrode at the various indicated times. Compared to slow-discharge profiles ( Fig. 7 a-d), the SOC is non-uniformly distributed during the fast-discharge rate. It can also be observed that, for the negative electrode, the differences in SOC mainly develop inside the particles (along the r n -axis) compared to the electrode thickness (along the x -axis). For example, toward the end of discharge in Fig. 8 d, although the SOC at the surface ( r n = 3 . 48 μm) is low, approximately 0.2, it remains high at the center of the HC particles ( r n = 0 ), approximately 0.8. In fact, the changes in SOC at the center of HC particles are insignificant at all times shown. This implies that the HC particles are too large for efficient charge transfer at fast and continuous discharge rates. Nevertheless, because of the high D 1 ,n [27] , this SOC at the cen-ter of the particles can still be recovered by setting the SIB in relaxation. Fig. 8 e-h shows the evolution of SOC profiles in the NVPF positive electrode along both the x -and r p -axis at various indicated times. During the fast-discharge rate, and similar to the SOC profiles shown in the negative electrode ( Fig. 8 a-d), the SOC profiles are also non-uniformly distributed. However, in contrast to the negative electrode profiles in Fig. 8 a-d, differences in SOC develop both inside the particles (along the r p -axis) and along the electrode thickness (along the x -axis). In addition, there are significant changes in the SOC at the center of the particles ( r p = 0 ) during discharge. The small NVPF particle radius (0.59 μm) compared to The SOC profiles along the x -axis in Fig. 8 e-h also mirror the electrolyte concentration profiles shown in Fig. 5 b. The low electrolyte concentration close to the positive electrode current collec-tor ( x = 157 μm) means the NVPF particles in this region are underutilized compared to the active material close to the separator. Such a variation of SOC along the electrode thickness is detrimental to the battery's performance because intercalated N a + cannot diffuse between adjacent particles. This issue is analogous to a cell  balancing problem in a battery module, although we are talking here of imbalances occurring along the electrode thickness. Consequently, increasing the electrode coating thickness without improving the electrolyte mass transport will result in greater imbalances and a huge penalty in terms of capacity loss for the HC//NVPF SIB. Fig. 9 shows the simulated results of ionic current density, i 2 in 2D color plots and the EMF of the positive electrode, U p in 1D plots as a function of the discharge time t (expressed as a percentage of t max ). The results of i 2 are also shown as a function of position x , and they are normalized with respect to i app , as indicated in the color code at the right-hand side of Fig. 9 a and c. Here, i 2 represents the flux of N a + due to migration and diffusion in the electrolyte phase (see Eq. (15) ). The results show that i 2 = 0 at the current collector ( x = 0 and x = 157 μm) and i 2 = i app at the separator ( 64 ≤ x ≤ 89 μm). This is in accord with the boundary conditions and thus validates the solution method in Eqs. (26) and 27. Fig. 9 a and b show i 2 profiles and the corresponding U p during 1 A m −2 discharge rate, respectively, while Fig. 9 c and d show i 2 and U p during 12 A m −2 discharge rate, respectively. It can be observed that i 2 profiles are linear in the negative electrode and nonlinear in the positive electrode, irrespective of the discharge rate. In addition, based on the side-by-side comparison of i 2 and U p , it can be observed that the nonlinear i 2 profiles in the positive elec-trode align with the step changes in the corresponding EMF of the NVPF. Therefore, the 'staircase' NVPF EMF results in nonlinear i 2 in the positive electrode. Results in Fig. 9 demonstrate that even in the cases of a slow discharge rate, the profiles of i 2 can be quite dynamic, which can pose a challenge for reduced-order models to be accurate in the case of SIBs. Fig. 10 shows a comparison of the experimental and simulation results of the HC//NVPF SIB in a Ragone plot. This figure compares the energy and power characteristics of the SIB. The simulation results are very close to the experimental results up to the 1-hour discharging rate. At higher rates, deviations appear, which can be explained by phase changes in the NVPF active material, which are not included in the solid-solution model. To improve accuracy at higher rates, a multi-phase diffusion mechanism is therefore needed to model phase transformations in the NVPF active material [ 60 , 61 ]. Nevertheless, the improvements brought by concentration-dependent D 1 ,m and k m result in a close match between the experiment and model predictions while maintaining a single set of parameters for all rates.
Often, in battery design, increasing the battery's energy density results in decreased power density. As a result, optimizing battery performance is nontrivial. However, the model herein presented can be used to determine design parameters, such as electrode thickness and porosity, based on the accuracy shown in the Ragone plot. At the same time, battery manufacturing costs should not be neglected, which can also be part of a multi-objective optimization procedure [62] . For example, increasing the coating thickness reduces the cost but simultaneously reduces the battery's power [ 63 , 64 ]. These factors can be combined and investigated using this SIB P2D model as a strategy and tool to avoid the often-expensive experimental trial and error methods.

Conclusions
A pseudo-two-dimensional (P2D) model is herein shown to model the voltage of a sodium-ion battery (SIB) composed of N a 3 V 2 ( P O 4 ) 2 F 3 (NVPF) and hard carbon (HC) as positive and negative electrodes, respectively. The HC//NVPF SIB model uses a coupled set of partial differential equations (PDEs) for the current and concentration profiles. An iterative root finding method is applied to determine the solution of the coupled system of PDEs. The model parameters are obtained by a genetic algorithm, a gradientfree optimization method. The negative and positive electrode parameters are optimized separately to reduce the number of simultaneously optimized parameters and improve accuracy. It is shown that the model is least accurate by 1.47% at a 1.4 C-rate using a constant set of parameters.
The developed P2D model can be rapidly parametrized using experimentally derived data. The voltage profiles for individual electrodes were obtained at different C-rates using a reference electrode and were used to determine parameters for each electrode. Using the validated P2D SIB model, more information concerning internal cell dynamics was obtained, which allowed an analysis of the limiting factors. It is shown that the high C-rate performance of the HC//NVPF SIB is limited by the poor mass transport in the HC and NVPF electrodes and in the electrolyte. Mass transport in HC electrodes can be improved by reducing the particle sizes. In contrast, the NVPF particles suffer from a low diffusion coefficient.
The model herein shown can be used as a design tool to improve the performances of SIBs, starting with the limiting fac-tors already identified. Future works will thus focus on multiobjective optimization of the cell design, including electrode thickness and material costs as additional design considerations. In addition, model accuracy can be improved by including temperature effects in large format cells and multi-phase intercalation dynamics in the NVPF electrode material. Another important aspect to consider is online parameter identification, estimation, and life prediction using the developed P2D model.

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.

Appendix A. Model parameters
This section lists the parameters used in the model and compares them to the experimental values where available.
Figs. A1 and A2 compare the concentration-dependent parameters, D 1 ,m and k m obtained from the P2D GITT model (symbols) [34] to the parameters used in the optimized full cell P2D model (solid line). In all cases, the obtained full cell model parameters are higher than the P2D GITT model parameters, although the same qualitative trend is maintained. These differences could be the result of model uncertainties and/or temperature effects at high discharge rates since the P2D GITT model parameters are obtained at comparatively very low currents (approximately C/30). were the same in both cases. Additional experimental measurements have validated the AEM electrolyte conductivity results [39] .    However, the D 2 values used in the optimized P2D model are qualitatively similar but quantitatively lower than the ones obtained from the AEM. Because two models are used to determine these parameters, further experimental investigations are necessary to investigate the origin of the differences.