Multi-Species, Multi-Reaction Model for Porous Intercalation Electrodes: Part I. Model Formulation and a Perturbation Solution for Low-Scan-Rate, Linear-Sweep Voltammetry of a Spinel Lithium Manganese Oxide Electrode

We formulate and implement a porous electrode model that includes multiple lithium-insertion species and associated electrochemical and homogenous reactions. This multi-species, multi-reaction (MSMR) model can account for the multiple current peaks observed in voltammetry data at slow scan rates and is also consistent with X-ray diffraction analysis of several different electrode materials. The MSMR model is used to simulate linear-sweep voltammetry data of a porous electrode of spinel lithium manganese oxide at different scan rates, and an analytic solution to these equations at very slow scan rates is given, based on a perturbation analysis. The greatest sensitivity at these slow scan rates is to the thermodynamic parameters, because the electrode is close to a state of dynamic equilibrium. The poor sensitivity to transport parameters means that it is not possible to distinguish between transport effects and a small amount of hysteresis that could be present in the slow-scan voltage measurements. These facts are illustrated using both the analytic solution and its comparison to numerical calculations. In Part 2, the analytic solution will be combined with numerical simulations at higher scan rates to complete the fitting process of transport, kinetic and thermodynamic parameters, along with a small amount of hysteresis. © The Author(s) 2018. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0771816jes]

In recent work, 1,2 it has been noted that current peaks during slowscan voltammetry of lithium host materials can be interpreted as evidence of multiple reactions involving different species of lithium in the host. To interpret the experimental data, the authors developed a model incorporating multiple lithium species and associated reactions (the MSMR model). The MSMR model provided quantitative agreement with experimental data for the open circuit voltage of lithium host materials of Si, 1 graphite, iron phosphate, and spinel manganese oxide. 2 The different species of lithium correspond to different galleries in the host materials and were shown to be consistent with X-ray analyses. The MSMR model describes the relationship between the Open Circuit Voltage (OCV) and the mole fraction of each species of intercalated lithium. The model also includes multiple Butler-Volmer equations to describe the kinetics of each lithium insertion reaction, as summarized below. A corresponding description of transport in the solid phase is potentially rather complicated and could also involve exchange reactions of lithium between the galleries. 2 In this work, we will make the simplifying assumption that all of these exchange reactions are in equilibrium, an assumption that allows one to derive a single lithium diffusion equation that is then coupled to multi-species kinetics equations.
We shall implement our approach on lithiated spinel manganese oxide (Li x Mn 2 O 4 , with 0 ≤ x ≤ 1, which we refer to as LMO) for the higher voltage plateaus that provide stable cycling. The material is used in our Gen 1 and Gen 2 Volts and is of strong interest to General Motors for advanced 12 V batteries capable of high discharge rates and regenerative braking, leading to substantial cell charging and energy capture over a driving profile. LMO is a well-studied electrode material. [4][5][6][7][8][9][10][11][12][13][14][15][16] Helpful reviews of various positive electrodes for lithiumbased batteries can be found in Refs. 18-20. A recent review of LMO is provided in Ref. 21. Our first goal in this work is to show how the resultant kinetics and solid-phase transport equations associated with the aforementioned MSMR model can be put into a porous electrode model, which has not been done previously, and to demonstrate the utility of this model by using it to describe linear-sweep voltammetry (LSV) data of an LMO electrode. A least squares optimization routine is then used to fit physical parameters for this model to the LSV data. A total number of twenty parameters were simultaneously fit to the highest and lowest scan rates, and the size of this problem poses a serious computational challenge, because each evaluation of the model for a specific set of parameters involves the numerical solution of a system of partial differential equations containing two independent spatial variables r and z, and time t.
The LSV data were taken over a range of scan rates from 0.05 mV/s to 0.5 mV/s. In addition to numerically solving the equations of the porous electrode model described above, a perturbation analysis of these equations was done in the limit of very slow scan rates, and an analytic solution was derived that predicts the current as a function of voltage for scan rates that yield a small perturbation from dynamic equilibrium, where dynamic equilibrium denotes the response when the system moves reversibly through a set of equilibrium states uniquely determined by the OCV at each state. The analytic solution can be used to approximate to good accuracy the current response at the lowest scan rate of 0.05 mV/s. (See Figures 2 and 3 and comments in the Discussion section.) By incorporating it into the solution strategy for parameter fitting, significant savings in computational time can be achieved. This is Part 1 of a two-part series on this subject. In Part 1, we describe the MSMR approach used to determine both the OCV and resistance associated with irreversible processes (e.g., charge-transfer kinetics and solid-phase diffusion resistances). The full set of equations for the resulting porous electrode model is given in dimensionless form in Table III. A perturbation analysis of these equations is then given to provide an analytic approximation to the current density as a function of voltage at low scan rates. This will then be compared to numerical simulations based on the complete equation system to illustrate its accuracy. The analytic solution also provides some insight into the sensitivity at low scan rates to different resistances. Part 2 of this series will describe the details of the (simultaneous) parameter fitting and the subsequent sensitivity analysis to the fit parameters.
In the Discussion Section, we summarize and consider some of the consequences of this analysis. Because the system is so close to equilibrium at the lowest scan rate, it is not very sensitive to the values chosen for transport and kinetics parameters, and the primary sensitivity is to the thermodynamic parameters that determine the OCV. One of the issues that arose was the inability to distinguish between a small amount of hysteresis in the OCV and the small impacts of A3953 transport and kinetics resistances at the lowest scan rate. These two effects differ from each other only in the fact that the impact of transport resistance grows with increasing scan rate whereas (OCV) hysteresis is rate independent. In contrast, at the highest scan rate, a small amount of hysteresis in the OCV has only a small impact, because the biggest difference between the forward and reverse scan rates is due to transport and kinetics resistances. The lowest scan rate is thus best suited for fitting thermodynamic parameters and the highest scan rate is best suited for fitting transport parameters; in this sense, the two complement each other and can be used together to distinguish between transport resistance and hysteresis. A more complete treatment of hysteresis can be found in Ref. 24, along with references to related work.

Model Description
Thermodynamics.-The thermodynamics and kinetics of the lithiation reaction are given using the methods developed in Refs. 1,2. References to related work can be found in Refs. 1,2, and Birkl et al. 3 use a similar methodology for the thermodynamic model in their recent publication. The mole fraction of sites occupied by lithium in the LMO is given as where x 1 and x 2 are the mole fractions of sites occupied in each of the galleries. Each of these species is involved in an insertion reaction at the interface between the host particles and the electrolyte of the form Thus, for each species j, a vacant host site H j can accommodate one lithium leading to a filled host site (Li − H) j . The OCV for this reaction is then written as X j represents the total fraction of available host sites which can be occupied by species j, U 0 j can be viewed as a standard electrode potential, and the parameter ω j is an additional parameter which can be fit to experimental data. Equation 3 represents the Reaction 2 in a quasi-Nernstian form, where the Nernst equation with unit activities is recovered when ω j = 1. When ω j << 1, the OCV U j (x j ) develops a plateau approaching two-phase behavior. The form of Equation 3 was chosen because it allows one to approximate OCV behavior of many lithium host materials, as shown in Ref. 1. At a particle interface with the electrolyte, local equilibrium requires that [4] where the common potential U (x) is the OCV of the electrode with respect to a lithium reference. For each gallery j, Equations 3 and 4 can be inverted to yield Combining Equations 1 and 5, we obtain an explicit closed form expression for the inverse of the OCV The parameters X j , U 0 j and ω j must be chosen to fit experimental measurements of OCV data. Similar procedures have been shown capable of reproducing OCV data for a variety of different lithium host materials. 1,2 The driving forces for diffusive transport in the solid phase are the gradients of the chemical potentials of each species μ j and the corresponding chemical potentials for vacancies μ H, j , and we seek an explicit form for these chemical potentials. An exchange reaction between species 1 and 2 in the particle interior can be written as [7] In this work, we assume that Reaction 7 is fast enough so that it is in local equilibrium everywhere within the host particles. The local equilibrium assumption implies that μ H,1 + μ 2 = μ H,2 + μ 1 [8] In addition, Reaction 2 implies that [9] where μ re f Li is the chemical potential of a lithium reference electrode, also assumed to be in the same electrolyte as the host material. Equations 8 and 9 imply that Equation 4 holds, even in the interior of the host particles where there is no contact with electrolyte.
The assumption of a locally equilibrated Reaction 7 allows one to define a total chemical potential as In a similar fashion, one defines [11] From Equations 9 and 4 one has [12] Since the volumetric expansion of LMO under lithium insertion is small, it will be neglected, so that one can ignore the impact of stress on the chemical potential. The differential of the Gibbs free energy at constant temperature and pressure then becomes The usual derivation of the Gibbs-Duhem relationship at constant temperature and pressure then yields Equations 12 and 14 then lead to Kinetics.-The kinetics of the insertion reaction are given as where i s,k is the current density due to reaction k, j is the total interfacial current density, β k are the symmetry factors for each reaction, 1 is the potential in the solid phase and 2 is the potential in the electrolyte phase with respect to a lithium reference electrode. R f is a film resistance due to a surface layer on the electrode, c 2 is the salt concentration in the electrolyte, and c 0 is the bulk concentration of c 2 .
Solid-phase transport.-Substitutional diffusion in a lithium host material consisting of two species, the occupied and unoccupied sites, can be described by Eq. 7 of Ref. 22 [17] where N is the flux of lithiated sites in the substitutional material, N H is the flux of unlithiated sites, D is a diffusion coefficient, which we will assume to be constant, c T is the total concentration of both lithiated and unlithiated sites (assumed constant), and RT is the product of the universal gas constant with temperature in Kelvin. Because we are ignoring volumetric expansion during lithiation, the total flux of sites in the host material can be assumed to vanish 22 Equations 15, 17 and 18 then yield It then follows that material balance of lithium in the solid phase can be given as It is assumed that the LMO particles are spherical so that, for radially symmetric diffusion, Equation 20 takes the form Boundary conditions for Equation 21 are where a is the particle radius. The overpotential η s in Equations 16 requires evaluation of the function U (x) and Equations 20-22 require the derivative dU/dx, but unfortunately these two functions cannot be explicitly evaluated, and the evaluations must be done by numerically inverting Equations 1 and 5. This problem can be avoided by replacing the dependent variable x in these equations with a new dependent variable U , subject to the transformation Equation 21 then takes the form  where Equations 22 then become Equations 1-26 can be incorporated into a porous-electrode half-cell battery model as shown in Figure 1. The interface between the counterelectrode and the separator is assumed to be at z = 0, the interface between the separator and the working electrode is at z = L 2 and the current collector of the working electrode is at z = L + L 2 . Table I gives the field equations and the boundary conditions for voltammetry of a porous electrode. This table assumes that a lithium reference electrode is placed at the interface between the counter-electrode and the separator, at z = 0, and that the liquid-phase potential 2 is measured with respect to this reference. Superficial current density in the solid phase of the porous electrode is given as i 1 and in the electrolyte phase as i 2 . The parameter ε 1,F represents the fraction of the solid-phase volume that is available for the lithium insertion reaction. R c is a contact resistance. The meaning of the remaining parameters can be found in the Nomenclature Index. The LSV potential excitation function is given by the function For this work, V min = 3.5 Vand V max = 4.25 V.
The equations in Table I are first order in i 1 , i 2 , 1 , and 2 , but for the purposes of numerical solution they were converted into a system of second order equations in the variables 1 and 2 . In order to do ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 207.241.232.120 Downloaded on 2019-11-14 to IP Table II. Second order system of dimensioned equations for the porous electrode (upper table) and the spherical particle (lower table). Sample parameter values are given in Table IV. Equations 85-88 give parameter relations that depend on concentrations or other variables. The current density at the particle surface j is given in Equation 16.
The corresponding system of second order equations is given in Table II. Table II in dimensionless form we introduce the following resistances [29] This leads to the following dimensionless quantities

Nondimensionalization, scaling, and dependent variable transformation.-Please refer to the Notation index for definitions of symbols. To put the equations from
The scaling parameter i s should be chosen to appropriately scale the interfacial current j. We have chosen the value for i s = i re f 0,1 for this purpose. The resulting dimensionless equations are listed in Table III. The dimensionless form of Equation 27 is The dimensionless forms of Equations 24-26 are Equations 16 becomē Table III. Dimensionless equations as a second order system for the porous electrode (upper table) and the spherical particle (lower table). The dimensionless current density at the particle surfacej is given in Equation 35. In contrast to Table II, the dependent variable for the particle isŪ instead of x.  Table III, were numerically solved using a linear voltage scan from 3.5 V to 4.25 V and back at different scan rates. The initial conditions at t = 0 were

Perturbation Solution of the Equations at Low Scan Rates
From a thermodynamic standpoint, if the voltage is scanned slowly enough, the electrode should move reversibly through the equilibrium states determined at each voltage, and this process is referred to as dynamic equilibrium. 23 The current density at dynamic equilibrium is given as where Q is the total charge capacity per geometric area of the electrode and υ is the voltage scan rate. At least two factors prohibit direct measurements of dynamic equilibrium: (1) side reactions in the electrolyte start to impact voltage measurements at very slow scan rates and (2) the dynamic equilibrium current densities i DE become too small to measure accurately. One must then try to correct measurements at the slowest sweep rates that are not influenced by these issues for irreversible effects in order to infer what the function x(U ) should be. The current densities i DE (t) can be simulated from the model equations in Table III in the limit as the time t 0 to sweep from the minimum to maximum voltage tends to infinity, or equivalently as the voltage scan rate υ tends to zero. This corresponds to the limit as the dimensionless parameters γ andγ tend to zero (see Equation 30), and it suggests that one do a perturbation analysis in the limit of small values of γ andγ in order to try to calculate correction terms to Equation 37 that account for the impact of irreversible processes. This is the analysis that we now describe. We start by restating the dimensionless equations in Table III in a form that facilitates the perturbation analysis and allows us to reference equations by equation numbers. In the separator 0 <z < 1: In the electrode 1 <z < 2:ī Equation 41 assumes β 1 = β 2 .
In the spherical particles 0 <r < 1: Boundary condition at the particle surfacer = 1: Boundary conditions at the separator interface with the counter elec- At the interfacez = 1 between separator and porous electrodē [47] At the current collector of the working electrodez = 2: Integration formulas.-If one integrates Equation 43 over the spherical particle, one obtains [49] If one then integrates Equation 49 over the porous electrode and uses Equation 40, one obtains Equation 50 will be the main tool that will be used to derive perturbation approximations for the current under slow-sweep voltammetry.
Perturbation assumptions.-Our interest is the limiting case when t 0 tends to infinity or equivalently when the scan rate υ tends to zero. In this limit, the dimensionless parameters γ,γ also tend to zero, and a perturbation analysis is done for γ,γ << 1. However, as t 0 → ∞ the ratioγ/γ remains fixed, and this is the assumption we make in the following analysis. We note that, generally speaking,γ/γ is a very small number, so that terms involvingγ will also be small if γ is small. In addition, the dimensionless ratios of the resistances R e , R s , R e,di f f , R k , R d , R sep (and indeed the resistances themselves) appearing in Equations 40-48 also remain constant as γ tends to zero.
Perturbation analysis.-We adopt the same notation for each dependent variable above to designate order of correction with an index number, for example, the mole fraction in gallery j will be written as [51] where . . . represents higher order terms in the expansion. For the total current density or total mole fraction, this notation can lead ) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 207.241.232.120 Downloaded on 2019-11-14 to IP A3957 to confusion; hence, a comma will be used for clarity before the subscripts of the perturbation approximation: [52] (If a comma did not precede the subscript, we would see terms such as i 1 andī 2 in the expansion for the total current, which could be confused with the solid and solution phase current densities, respectively.) As γ tends to zero, the system approaches a state of dynamic equilibrium, where U = V (t). It follows that Under conditions of dynamic equilibrium, Equation 50 yields Also under conditions of dynamic equilibrium, We turn now to calculating the next higher order terms inc 2 and i . Equations 55 and 42 and the boundary condition 48 imply that Equation 41 also implies that With these calculations, one can now use Equations 40-48 to calculate the next higher order perturbation terms: Boundary condition at the particle surfacer = 1: The boundary condition at the separator interface with the counter electrodez = 0 is given by¯ The boundary conditions at the separator electrode interfacez = 1 correspond to At the current collector of the working electrodez = 2, the boundary conditions are¯ 1,1 =ī 1,1

Rc
Rs From Equation 62, it follows thatj 1 is constant in the porous electrode and not a function ofz, and this makes it easy to integrate the first of Equation 59: but from Equations 59 and 62 ) unless CC License in place (see abstract [73] From Equations 71-73 it follows that From Equation 53 We compactify Equation 75 with the following definitions: and Equation 76 becomes Lastly, we note that From Equations 54 and 81, it follows that [82] To express this approximation in dimensioned form, note that Then, in dimensioned form, Equation 82 becomes 0.01 S/cm Measured * Corresponds to 118 mAh/g measured capacity for the active material 2 and an Mn 2 O 4 molecular weight and density of 173.9 g/mol and 4.54 g/cm 3 , respectively. parison to i DE . Each of these resistances in the correction term is associated to the corresponding type of transport loss. One last comment needs to be made about the validity of the above analysis. Equation 69 implies that the expression x 1 (r ) has a fully-developed concentration profile, i.e., the same profile seen at long times when lithium is diffusing in from the boundary at a constant rate. This profile takes a short time to develop (on the order of γ in dimensionless time) when starting from a state of static, as opposed to dynamic, equilibrium. During scan reversal, this profile will again be interrupted, but develops again after a dimensionless time period on the order of γ. It is important to understand that in these transitional periods after the start of LSV and after scan reversal, the above analysis does not hold. More discussion of this is given in the appendix of Ref. 25. See also the discussion of Figures 3c and 3d below, where the period after scan reversal is clearly visible in the curves done by numerical simulation. The following relationships were used in solving the model equations. Christensen 26 provides a short review of the literature relative to the quantities provided in Equations 85 for the solution-phase conductivity and lithium ion transference number: Last, as is commonly done, we employ the following relation,

Accuracy of the Perturbation Analysis and Sensitivity to Parameters
which is valid for spherical particles. Additional parameter values used are also given in Table IV. When simulating LSV, the model must calculate current as a function of voltage. In this regard, it is helpful to note that the total current density is given as The dimensionless equations in Table III were solved using a Fortran code based on finite differences. 32 Figure 2 shows a comparison of the numerical solution a for current density given in Equation 89 and the approximate expression given in Equation 84 at a scan rate of 0.05 mV/s. Also shown is the curve for dynamic equilibrium. One sees that, although the current is very close to dynamic equilibrium, a small amount of transport losses is also present. Figure 3 shows comparisons a The accuracy of the numerical solutions will be addressed in Part 2 of this work, where both lower and higher scan rates are considered. between numerical solution and Equation 84 for the correction term i − i DE , which is the difference between the total current and the current at dynamic equilibrium, given by Equations 37 and 84, at scan rates of 0.05 and 0.1 mV/s. Fig. 3a is the forward scan direction at 0.05 mV/s and Fig. 3b is the same plot at 0.1 mV/s. One sees clearly that the error in the perturbation formula increases with scan rate. Fig. 3c shows the reverse scan direction at 0.05 mV/s and Fig. 3d is the same plot at 0.1 mV/s. At 0.05 mV/s, one sees a good agreement between the numerics and the analytic approximation, except in a small voltage region directly after scan reversal in Fig. 3c. As noted above, there is a brief period of time after scan reversal before the lithium concentration profiles in the particles reach a fully-developed state, as assumed in the analytic solution. Scan reversal occurs at 4.25 V, and Fig. 3c indicates that the fully developed concentration profiles are present at voltages lower than approximately 4.22 V. Figs. 3e and 3f, show the same plots as Figs. 3a and 3c, except that the solid phase diffusion coefficient was decreased by a factor of 5 from its value in Table IV, thereby increasing the value of R d also by a factor of 5. Again one sees that the errors in the perturbation solution grow, not only with scan rate, but also with the sizes of the various resistances that appear in formula 84; thus the error in Figures 2e and 2f are larger than in Figures 3a and 3c.
In part 2 of this work, we will compare the analytical solution 84 to measured voltammetry data, and one would like to know how sensitive the current in Equation 84 is to the values of the individual resistances that appear in the formula. In fact, it is not possible to discern between the effects of R e , R s , R sep , R c , and R f in the formula by comparing to measured voltammetry, because they all impact the current density in the same way. This is because the coef- ficients A 2,s , A 2,e , A 3 , A 2,c and A 4  independent of the mole fraction x 0 of lithium in the particles. One can thus only determine the linear combination by comparing Equation 84 to measured voltammetry data. In Equation 90, the resistances R c and R f drop out, because they were set to zero in Table IV. The remaining terms in Equation 90 are shown in Table V. Clearly the solid phase resistance R s is also negligible. It must be emphasized, however, that the inability to distinguish between these resistances only holds at the very low scan rates considered here. A detailed analysis of the sensitivity to these resistances at higher scan rates goes beyond the scope of this work, but it seems likely that current data at higher scan rates exhibit different sensitivities to each of the resistances in R total . The correction term i − i DE in Equation 84 can be rewritten as a linear sum of terms, each of which is proportional to one of the resistancesR total , R k , and R d . Figure 4 shows each of the summands in the current density associated with R total , R k , and R d in Equation 84 as a function of voltage. One sees that the impact of each of these resistances on current density depends on the particular voltage value, but the kinetic resistance appears to have much less impact than that of the solid-phase diffusion or R total . It is also important to note that each of the terms plotted in Figure 4 is independent of the direction in which voltage is scanned. This is not true of the current density i DE at dynamic equilibrium, which changes sign when the scan rate reverses.

Discussion of Results
Because we have assumed that the inter-species bulk reaction given in Equation 7 is always in a state of local equilibrium, the diffusion Equation 20 is no different from what it would be in a standard treatment without multiple species (see for example Ref. 22). The assumption of local equilibrium also means that the surface overpotentials for the reactions in Equations 16 are equal to each other, but if the symmetry factors for each reaction are different, then Equation 16 are not equivalent to a single Butler-Volmer equation for a single species. Thus the major difference between the treatment here and a single-species analysis is that one has four parameters to fit to kinetics (β 1 , β 2 , i re f 0,1 , i re f 0,2 ) instead of only a single symmetry factor and exchange current density. Although it is not assumed in the current model, it is possible that, for example, at high current densities the inter-species Reaction 7 no longer satisfies local equilibrium. In this case, a single diffusion equation no longer suffices to describe transport and a number of other possibilities must be considered. In addition, the kinetics expressions would then be given by two independent interfacial reactions, one for each species.  Figure 5, and the right plot provides the ratio of the lithium filled site fraction of particles at the separator-electrode interface,z = 1, to its average value atz = 1. (The profiles in the right plots can be compared to Equations 69 and 70 of the perturbation analysis, which exhibit a quadratic dependence onr .) We see that the ratio is largest at the end of the forward potential scan, when the fraction of filled sites at the center of the particles is about 17 percent higher than that of the average atz = 1, and the fraction of filled sites at the particles' surface is about 10 percent lower than the same average. Despite this 27 percent difference in the fraction of filled sites throughout the particles, the perturbation analysis does provide good accuracy for the calculated current density at the end of the forward scan (cf. Figures 2 and 3a, V = 4.25 V). The highest current density during the forward scan is near τ = 0.9 and V = 4.15 V; as seen in Figures 5, there is about a 20 percent difference in the filled site fractions between the particle center and surface for particles at the separator-electrode interface at τ = 0.9, and, again, the perturbation analysis maintains reasonable accuracy (cf. Figures  2 and 3a; the peak current is near V = 4.15 V).   Figure 5).
Three dimensional surface-contour plots of the filled site fraction x (left) and the potential U(x) (see Equation 6), the dependent variable employed in the numerical calculations (see Table III) are provided in Figure 6 for the 0.05 mV/s scan rate and τ = 0.9, near the peak current density of the forward scan, as discussed immediately above. The smooth surfaces are indicative of why accurate numerical calculations are possible. We see much less variation in x and U(x) in the zdirection than in the r-direction, which is consistent with R d being the largest resistance in Table V. Figure 7 provides surface-contour plots of the filled site fraction x at the end of the forward scan (left) and near the start of the backward scan for the numerically calculated voltammogram of Figure 2. The perturbation solution works well for the entire 0.05 mV/s forward scan, as the currents begin with near zero values at 3.5 V vs. Li, and they reach their fully developed profile [quadratic, see Equation 69, for the correction term x 1 (r )] for the remainder of the forward scan (left plot of Figure 7). However, upon scan reversal (right plot of Figure 7), the x(r) profiles are not quadratic in r and are not fully developed, and, as discussed in the context of Figure 3c, the perturbation solution is inaccurate for the initial portions of the backward scan. Figure 2 shows clearly that there is only a small sensitivity to transport resistance at the very low scan rate of 0.05 mV/s, the system is close to dynamic equilibrium, and the perturbation solution is close to that of the numerical solution to the full system of equations. At the same time, this reduced sensitivity to transport makes it possible Table VI. Thermodynamic values with hysteresis, shown in Figure 8. to confuse the effects of transport with a small hysteresis effect in the OCV. Figure 8 shows current densities for two different sets of parameter values, based on Equation 84, both of which produce close to identical current densities. One set assumes the transport and thermodynamic values in Table IV, which exhibit no hysteresis. The second set of values assumes that all transport resistances R total , R k , and R d are zero, but allows for hysteresis in the thermodynamic values, as shown in Table VI. The hysteresis shown in Figure 5 is represented by choosing a different set of thermodynamic parameters for the forward and backward scan rates, so that the function x(U ) now depends on the scanning direction. Figure 7 underscores the fact that at very slow scan rates, one cannot distinguish between the effects of transport and kinetics resistances and a small amount of hysteresis. A detailed model of hysteresis is complicated, for example, by the fact that the hysteretic   Table IV, and assuming dynamic equilibrium (no transport losses) but assuming hysteresis with thermodynamic values taken from Table VI. The vertical dashed lines show the difference in the voltages of the peak current densities in the forward and reverse scans and correspond the U 0 i values indicated and in Table VI. OCV depends on when the voltage scan is reversed, and whether it is reversed during LSV or scanning at constant current. 24 The model we have used here simply fits different thermodynamic parameters to the forward and backward scan rates and assumes that these values determine a direction-dependent OCV-function x(U ). The simulation in Figure 8 determined hysteresis by assuming that there was no transport present, i.e. the scan at 0.05 mV/s was in dynamic equilibrium. On the other hand, it is also possible to assume other nonzero values for transport parameters and to fit a different hysteretic function x(U ) to the same current data. Thus there does not appear to be any way to distinguish between the effects of transport resistance and hysteresis based solely on fitting data to the single scan rate of 0.05 mV/s. In Part 2 of this work, we will show that LSV data at higher scan rates can resolve this ambiguity in the transport parameters, because the higher scan rates have much greater sensitivity to transport and less sensitivity to hysteresis than at the low scan rates considered here. By fitting to voltammetry at both high and low scan rates, it will be possible to determine more accurate values for both the transport parameters and for a function x(U ) that exhibits a small amount of hysteresis.

Conclusions
This is Part 1 of a two-part series describing how to fit model parameters to linear sweep voltammetry (LSV) data taken with a lithium manganese oxide electrode (LMO) electrode. The previously developed multi-species, multi-reaction (MSMR) model has been incorporated into a porous electrode model capable of describing the behavior of many materials of current commercial interest. In Part 2 of this series, this model will be used to simulate LSV data at sweep rates between 0.05 and 0.5 mV/s for an LMO electrode. Solution of the resulting system of partial differential equations (PDEs) must be done numerically, and the resulting calculations require substantial amounts of computation time during the process of optimizing parameter values to measured data. However, the perturbation analysis in Part 1 gives an analytic formula for current density as a function of voltage during LSV in the limit of very slow scanning rates. This formula offers insights into parameter sensitivity at slow scanning rates and it can significantly reduce the amounts of computation time needed to fit parameters to data.
Using estimated parameter values for LMO electrodes, taken from the literature, (see Table IV) it appears that the analytic formula provides sufficient accuracy at the lowest scan rate of 0.05 mV/s to be useful in the parameter fitting process. However, numerical simulation of the full set of PDEs is needed at higher sweep rates. Comparison of the analytic formula with numerical solutions illustrates how the accuracy of the analytic formula is impacted by increasing the sweep rate or by increased resistances (associated with irreversible processes) at a given sweep rate.
At very slow sweep rates, the electrode is close to a state of dynamic equilibrium, so that the sensitivity to transport and kinetic parameters is very poor. In addition, it is not possible to distinguish between the effects of a small amount of apparent hysteresis in the opencircuit voltage and the effects of transport and kinetic resistance at these slow rates. In contrast, at higher sweep rates, the sensitivity to transport increases and the impact of a small amount of hysteresis decreases; thus the fitting process can be enhanced by simultaneously considering both slow and fast sweep rates.
In Part 2 of this work, the experimental details of the voltammetry will be given, as well as the details of how to fit parameters using a least squares optimization of data based on numerical solutions of the fastest scan rate (0.5 mV/s) as well as the analytic solution at the lowest scan rate (0.05 mV/s).

Acknowledgment
The authors wish to recognize Wentian Gu, former General Motors researcher within our China Science Lab, for providing a link to our CSL electrode fabrication facility, taking data on LMO electrodes, and related discussions.