Modeling the Mesoscale Transport of Lithium-Magnetite Electrodes Using Insight from Discharge and Voltage Recovery Experiments

A multi-scale mathematical model, which accounts for mass transport on the crystal and agglomerate length-scales, is used to investigatetheelectrochemicalperformanceoflithium-magnetiteelectrochemicalcells.Experimentaldischargeandvoltagerecoverydataarecomparedtothreesetsofsimulations,whichincorporatecrystal-only,agglomerate-only,ormulti-scaletransporteffects.Masstransportdiffusioncoefﬁcientsaredeterminedbyﬁttingthesimulatedvoltagerecoverytimestoexperimentaldata.Inaddition,afurtherextensionofthemulti-scalemodelisproposedwhichaccountsfortheimpactofagglomeratesizedistributionsonelectrochemicalperformance.Theresultsofthestudyindicatethat,dependingonthecrystalsize,thelowutilizationoftheactivematerialiscausedbytransportlimitationsontheagglomerateand/orcrystallength-scales.Forelectrodescomposedofsmallcrystals(6and8nmdiameters),itisconcludedthatthetransportlimitationsintheagglomerateareprimarilyresponsibleforthelongvoltage recovery times and low utilization of the active mass. In the electrodes composed of large crystals (32 nm diameter), the slow voltage recovery is attributed to transport limitations on both the agglomerate and crystal length-scales. In the present a more in depth study is conducted using the multi-scale model. The performance-limiting processes within the agglomerate and crystal length-scales of a Li/Fe 3 O 4 are simulated and compared to experimental discharge and voltage recovery data for electrodes composed of 6, 8, and 32 nm nanocrystals. Results provide information on how

Large increases in the use of distributed and intermittent energy sources (i.e., wind and solar) have increased the need for cost effective, reliable, and efficient energy storage technologies. 1[7][8][9][10][11][12][13][14] Despite these advantages, one of the major challenges limiting the advancement of magnetite electrodes is a considerable difference between the maximum, theoretical capacity and the observed, experimental capacity of the active material.This difference increases the anticipated cost of magnetite batteries because it requires the electrodes to be overdesigned with excess amounts of active material.
The difference between the theoretical and experimental capacity is related to the close-packed inverse spinel structure of Fe 3 O 4 , which restricts the transport of lithium in the material.0][11][12][13][14] The smaller path length increases the utilization of the active material by making it possible for ions to penetrate to the center of the crystals, especially at high rates of discharge.Electrodes fabricated with nano-crystalline magnetite have shown significant improvement in capacity; however, the theoretical capacity has still proven difficult to obtain. 11Further improvements in capacity may require a more detailed understanding of the ancillary effects associated with fabricating an electrode from nano-crystalline active materials.For instance, due to the materials synthesis and electrode fabrication processes, Fe 3 O 4 nanocrystals tend to form micron-sized agglomerates. 15These agglomerates could decrease the utilization of the active material by hindering ion transport toward the crystals at the center of the agglomerate.
At present, it has been difficult to directly quantify the impact of agglomerates on electrochemical performance due to the complex structure of the battery electrodes.With the addition of agglomerates, there are three length-scales within the electrode that can impact the overall battery performance: the bulk electrode (macro-scale), the agglomerates (mesoscale), and the crystals (nanoscale).An understanding of the processes, especially ion transport, occurring on all three length-scales is needed to further optimize the nanocomposite magnetite electrodes.
One way to help clarify which physical processes on which lengthscales influence the battery performance is through the development of a mathematical model with predictive capabilities.A variety of modeling efforts exist in the literature for a variety of lithium ion batteries and electrodes; [16][17][18][19][20][21][22][23][24][25][26][27][28][29] however, to the best of the authors' knowledge, there have been no attempts to simulate the performance of a Fe 3 O 4 electrode.In addition, most of the current models only account for the physical processes in the bulk electrode and the solid crystal because the investigated electrode materials tended to have larger crystal sizes (>50 nm in diameter) which do not readily form agglomerates.One exception is the work of Dargaville and Farrell, which simulated the performance of a lithium iron phosphate battery and included an agglomerate length-scale. 19In that work, the authors used insight from the experimental literature to assume the FePO 4 crystals formed porous agglomerates. 30,31They concluded from simulations that the agglomerates only impacted the electrochemical performance at high rates of discharge.In contrast, transmission electron microscope images of Fe 3 O 4 electrodes indicate that the nanocrystals, which typically have diameters of 8 to 32 nm, form tightly packed agglomerates with small void spaces. 15These observations suggest that the Fe 3 O 4 agglomerates could have a more significant impact on the electrochemical performance of the electrode.
This work seeks to investigate the performance-limiting processes of a magnetite electrode using a mathematical model that was developed with insight from experimental discharge and voltage recovery data.Recent voltage recovery experiments performed by Zhu et.al. have shown that electrodes fabricated with nano-crystalline magnetite take over 100 hours to reach an equilibrium voltage in response to current interruption. 11In previous work, it was suggested that these long voltage recovery times were caused by the relaxation of concentration distributions within the agglomerate and/or crystal lengthscales of the electrode. 32To further investigate this behavior, a multiscale mathematical model was proposed, which was validated against discharge and recovery data from electrodes comprised of 6 nm crystals.In the present work, a more in depth study is conducted using the multi-scale model.The performance-limiting processes within the agglomerate and crystal length-scales of a Li/Fe 3 O 4 electrode are simulated and compared to experimental discharge and voltage recovery data for electrodes composed of 6, 8, and 32 nm nanocrystals.Results provide information on how the ion transport on each length-scale impacts the electrochemical performance.In addition, an extension to the model is proposed which accounts for the influence of distributions in agglomerate size.

Method of Approach
Experimental.-Nanocrystallinemagnetite (6 and 8 nm), Fe 3 O 4 , was synthesized using a co-precipitation approach, utilizing aqueous solutions of iron (III) chloride hexahydrate, iron (II) chloride hexahydrate, and base according to a previously reported method. 11,12arger sized nanocrystalline magnetite, ∼32 nm, was purchased from Alfa Aesar.X-ray diffraction data was collected using a Rigaku Smart Lab diffractometer with Cu Kα radiation.The crystallite sizes of the Fe 3 O 4 powders were calculated by applying the Scherrer equation to the FWHM of the (311) peak. 33An instrumental broadening correction was applied using a LaB 6 standard.
Electrodes were prepared using magnetite, carbon, and polyvinylidene fluoride binder coated onto an aluminum foil substrate.Electrochemical tests were performed using two electrode coin-type experimental cells with lithium metal anodes and 1 M LiPF 6 in dimethyl carbonate:ethylene carbonate electrolyte.The electrodes were comprised by weight of 90% Fe 3 O 4, 5% acetylene carbon black, and 5% PVDF.Discharge was conducted with no preconditioning under a C/200 rate to 0.5, 1.0, 1.5, 2.0 and 2.5 electron equivalents per Fe 3 O 4 and then allowed to rest under open circuit conditions for up to 30 days.A total of 15 coin-cells were used, one for each depth of discharge and each crystal size (6, 8, and 32 nm).Good agreement was observed between the discharge curves for cells comprised of the same crystal size (see Figure 1a in Ref. 32).The magnetite electrodes had a thickness of 51 ± 4 μm and an active mass loading of 4.2 ± 0.3 mg cm −2 .All electrochemical testing was conducted at 30 • C.
Cross sectional TEM images of the Fe 3 O 4 electrodes were acquired by embedding the electrode samples in an epoxy resin.A Reichert-Jung UltracutE ultramicrotome was used to slice 80 nm sections of the embedded electrodes for TEM analysis.Sections were viewed with a FEI Tecnai12 BioTwinG 2 transmission electron microscope.Digital images were acquired with an AMT XR-60 CCD Digital Camera system.The public domain Java image processing program Image J was used to determine the size and agglomerate distributions from the TEM cross sectional images. 34Examples of the TEM images used in the analysis can be found in. 15,32deling.-Simulationswere conducted using a multi-scale model developed by the authors that incorporates the transport of lithium in the agglomerate and crystal length-scales based on Fickian diffusion. 32On the crystal scale, transport processes were also simulated using concentrated solution theory, [24][25][26][27] but significant improvements (when compared to experiments) were not found.At the surface of the crystals, the reaction-kinetics is modeled using the Butler-Volmer equation, which was formulated for a lithium-intercalation reaction.The reaction rate constant was approximated from current interrupt data from electrodes with magnetite crystals with diameters of 6, 8, and 32 nm.Changes in thermodynamic potential due to changes in the degree of lithium intercalation are accounted for by fitting experimental data to an equation for the equilibrium potential.The experimental data for the equilibrium potential was obtained by using the maximum voltage during voltage recovery experiments at various depths of discharge.The 32 nm experiments were not used to determine the equilibrium potential because, even after 700 hours, the voltage during recovery had yet to reach a plateau (see Figure 2).The equilibrium potential equation, which was derived from thermodynamic principles, resembles the Nernst equation but includes large corrections. 26The governing equations and boundary conditions for the model are given in Table I.Definitions of the variables can be found in Appendix B. The detailed description of the model formulation along with model parameters, constants, assumptions, and the equilibrium potential equation can be found in Ref. 32.
)] - In this study, simplifications to the multi-scale model were made in order to develop crystal-only and agglomerate-only models (Figure 1).Results from all three models were compared to experimental data in order to understand which length-scales were responsible for the observed trends in electrochemical performance.The crystalonly model was developed by assuming a fast diffusion coefficient of lithium on the agglomerate scale (10 −6 cm 2 s −1 ), which, under the current experimental conditions, yielded negligible variations in lithium-ion concentration throughout the agglomerate.Likewise, the agglomerate-only model was developed by assuming a fast diffusion coefficient for lithium on the crystal scale (10 −12 cm 2 s −1 ).

Results and Discussion
Voltage recovery experiments.-Figures2a and 2b show the voltage recovery data for experiments conducted with electrodes comprised of crystals with diameters of 8 and 32 nm, respectively.Data for 6 nm diameter crystals can be found in. 32For each set of data, the cells were discharged at C/200 to 0.5, 1.0, 1.5, 2.0 and 2.5 electron equivalents per Fe 3 O 4 prior to the observed voltage recovery.During voltage recovery of electrodes with 8 nm crystals, the voltage of electrodes that were discharged to low electron equivalents (0.5, 1.0, and 1.5) goes through a maximum before falling to an equilibrium value.A similar behavior is observed for electrodes with 6 nm crystals that were discharged to low electron equivalents. 32For all other sizes and levels of discharge, the maximum voltage occurs at the end of the recovery experiment.
Recovery times greater than 200 hours are observed for all five depths of discharge and all crystal sizes.These long recovery times could be associated with a slow phase transition occurring within the material.However, this is unlikely because the phase transition of magnetite from an inverse cubic spinel (Li x Fe 3 O 4 ) to a rock-salt like phase (LiFeO 2 ) does not occur until between x = 2.8 and x = 4.0 (for x in Li x Fe 3 O 4 ). 14Instead, these long times are likely caused by the slow relaxation of concentration profiles to a uniform value.A similar assumption provides the basis for determining mass transport parameters in solids (i.e., diffusion coefficients) using the galvanostatic intermittent titration technique (GITT). 35,36Based on this assumption, the long voltage recovery times indicate large non-uniformities in the lithium concentration within the electrode at the end of discharge.This suggests a poor utilization of the active material, where only a fraction of the material participates in the reaction.
The relaxation of concentration profiles in the magnetite electrode can occur within three length-scales: the bulk electrode, agglomerate, and crystal scales. 15In order to determine which length-scale is responsible for the long recovery time (and poor active material utilization), a comparison of the mass transport time-constants was conducted.To accomplish this, the relaxation time of the concentration profile on each length-scale was estimated using the time-constant, τ, which characterizes the time required for the concentration in the system to relax after a step change in concentration at one boundary.This time-constant is defined as follows: where is the characteristic length and D is the mass transfer diffusion coefficient.In our previous work, the time-constants for each length-scale were determined based on the physical properties of the magnetite electrodes used in the voltage recovery experiments.The resulting analysis suggested that mass transport on the agglomerate and/or crystal length scales could be responsible for the slow voltage recovery. 32mparison of models to experimental data.-Tofurther explore whether mass transport effects on the agglomerate and/or crystal length-scales are responsible for the long voltage recovery times, simulations from three different models (i.e., crystal-only, agglomerateonly, and multi-scale) were compared to experimental data.To accomplish this, the mass transport diffusion coefficients were first obtained by fitting each model to the 6 nm experiments.The resulting values were kept constant for the 8 and 32 nm simulations in order to gauge how each model predicted the observed trends in experimental data.For the agglomerate-only and multi-scale models, the agglomerate radius (r agg ) was 1.05 μm, which was obtained by taking the average of all agglomerate sizes from transmission electron micrographs of the cross-section of two electrodes fabricated with 8 and 32 nm crystals, respectively.This value was used for all simulations, unless specified otherwise.
To the best of the authors' knowledge, there are no reported values of lithium-ion diffusion coefficients in magnetite crystals (or agglomerates of crystals) reported in the literature.Therefore, the mass transport diffusion coefficients for each model were determined by fitting the simulated voltage recovery times to experimental data.The diffusion coefficients were selected to give the best agreement over the entire range of simulations (i.e., for discharges to 0.5, 1.0, 1.5, 2.0, and 2.5 electron equivalents per Fe 3 O 4 ).For the fitting procedure, the recovery time was defined as the time it takes the voltage to reach   Table II contains the diffusion coefficients obtained for each model using the fitting procedure.For the multi-scale model, the diffusion coefficient in the agglomerate (D agg ) was set equal to D agg from the agglomerate-only model.The diffusion coefficient in the crystal (D x ) was determined by selecting the lowest possible value that did not have an impact on the simulated discharge or voltage recovery.Obtaining the diffusion coefficient in this manner yielded the best agreement for the multi-scale model with the discharge and voltage recovery data, including trends in electrochemical performance with changes in crystal size.
Note that the agglomerate diffusion coefficient used in the simulations is within the range of experimentally reported values for the solid-state lithium diffusion coefficient in commercial lithium-ion materials (10 −8 to 10 −18 cm 2 s −1 ). 37This value is much lower than the agglomerate diffusion coefficient (∼10 −6 cm 2 s −1 ) used by other authors to simulate lithium transport through loosely-packed agglomerates of FePO 4 . 19The low diffusion coefficient likely results from the tight packing of the nano-crystals within the Fe 3 O 4 agglomerates, which has been observed using transmission electron microscopy. 15urthermore, assuming the crystals are close-packed, the maximum and minimum size of the void spaces for ion transport in the agglomerate can be determined from the size of the octahedral and trigonal void spaces, respectively.Figure 4 shows how the size of the crystals impacts the size of the void spaces in the agglomerate.The upper and lower bounds of the highlighted region were determined from a geometric analysis of the void spaces.The size of the void space is determined from the diameter of the largest sphere capable of fitting in the void (see equations in Figure 4).For the magnetite experiments, the maximum crystal radius was 16 nm, which suggests that the largest void spaces range from ∼5 to 13 nm.In this range, it is likely that ionsurface (as opposed to ion-ion or ion-solvent) interactions dictate the rate of mass transport through the agglomerate.Therefore, diffusion coefficients in the agglomerate are expected to be significantly lower than those obtained by using porosity/tortuosity corrections, which inherently assume ion-solvent interactions dominate.
In addition, the agglomerate diffusion coefficient is five orders of magnitude higher than the solid-state lithium diffusion coefficient used in the multi-scale simulations.9][40] This suggests that the mechanism responsible for diffusion through the agglomerate may be similar to the mechanisms associated with grain boundary and surface diffusion.Additionally, this trend agrees with the observations of Wang et.al., who tracked lithium transport and conversion in FeF 2 nanoparticles using in-situ transmission electron microscopy and concluded that diffusion along the surface of the nanoparticles was much quicker than the diffusion in the bulk material. 41igure 5 shows the maximum voltage change during the voltage recovery ( V max ) for all three experiments (6, 8, and 32 nm crystals) and all three models (crystal-only, agglomerate-only, and multi-scale).
For each set of experimental data, there were slight variations in the applied current density due to variations in the active mass.The current densities for the 6, 8 and 32 nm datasets were 4.6 ± 0.1, 4.6 ± 0.1, and 4.6 ± 0.2 mA g −1 (mean ± standard deviation), respectively.To account for these variations, three simulations were conducted for each crystal size.In the figure, each simulation curve was obtained using the average current density with the error bars corresponding to the V max obtained from simulations using the maximum and minimum current densities.
According to the experimental data, there is almost no difference in V max for the 6 and 8 nm crystals and a significant increase in V max for the 32 nm crystals.The results of Figure 5 indicate that only the multi-scale model is able to predict these trends (Fig. 5c).For instance, the crystal-only model (Fig. 5a) over-predicts the change in V max due to changes in crystal size.It also cannot simulate discharges past an average of 0.4 electron equivalents per Fe 3 O 4 for the 32 nm crystals because, at this point, the model predicts that the surface of the crystal will be fully lithiated (8 electron equivalents).In addition, the agglomerate-only model (Fig. 5b) significantly under-predicts the increase in V max for the 32 nm crystals.
A similar behavior is observed for the predicted discharge curves.Figure 6 shows the discharge curves for all three experiments and all three models.The discharges were conducted to a cutoff voltage of 1.5 V with current densities of 4.4, 4.7, and 4.5 mA g −1 for the 6, 8, and 32 nm cases (simulations and experiments), respectively.Figure 6a shows that the crystal-only model over-predicts the changes in discharge time due to changes in crystal size.It also predicts a change in the discharge time between the 6 and 8 nm data which is not observed experimentally.Figure 6b shows that the agglomerate-only model does not predict the observed decrease in discharge time for the 32 nm crystals.In fact, all three simulations for the agglomerate model produce the same result.Only the multi-scale model (Fig. 6c) is able to predict the similar discharge time for the 6 and 8 nm crystals and the decrease in discharge time for the 32 nm crystals.
The results of the comparison of the models suggest that the long voltage recovery times of the Fe 3 O 4 electrodes are caused by the relaxation of concentration distributions on both the agglomerate and crystal length-scales.These concentration distributions arise from mass transport limitations within the electrode.For the electrodes comprised of 6 and 8 nm crystals, recovery is caused by concentration   relaxation on the agglomerate scale.This explains why there is no variation in the discharge or recovery behavior between the two crystal sizes.For the electrodes composed of 32 nm crystals, recovery is caused by concentration relaxation on both the agglomerate and crystal scales.The crystal-scale becomes a factor when going from 8 to 32 nm crystals because, according to Eq. 1, this corresponds to a 16-fold increase in the mass-transport time constant.
Multi-scale model results.-In the previous section, it was concluded that the multi-scale model provides the best agreement with the Fe 3 O 4 voltage recovery experiments.This section provides a more thorough analysis of the multi-scale simulation results.Figure 7 compares the voltage recovery after discharges to 1.5, 2.0, and 2.5 electron equivalents per Fe 3 O 4 for the 8 and 32 nm data.Good agreement is observed between the experimental and simulated results for the 8 nm data.For the 32 nm data, fairly good agreement is observed for the final voltage, but there are discrepancies between the recovery times to reach the final voltage.The origin of this discrepancy can be identified through an analysis of the predicted concentration distributions.
For instance, Figures 8 and 9 contain predicted concentration distributions within the electrodes composed of 8 and 32 nm crystals during voltage recovery after discharge to 1.5 electron equivalents.Figure 8 shows the predicted distributions of the average solid-state lithium concentration within the agglomerate, and Figure 9 shows the predicted distributions of the solid-state lithium within the crystal at the agglomerate surface.In the figures, c agg is the concentration of lithium-ions in the agglomerate, c 0 is the bulk concentration of lithium-ions in the electrolyte, and c x is the concentration of intercalated lithium in each crystal.r and r are the radial positions within the agglomerate and crystal, respectively, and r agg and r x are the radii of the agglomerate and crystal, respectively.In Fig. 8, there is little difference between the 8 and 32 nm simulations.Both show an equilibration of the solid-state lithium concentration in the agglomerate within 200 hours.Because the model assumes there is no direct crystal to crystal exchange of solid-state lithium, the predicted equilibration is caused by lithium transport within the void spaces of the agglomerate.Near the surface of the agglomerate, the solid-state lithium is oxidized to produce mobile lithium-ions, which diffuse to the center of the agglomerate and subsequently reduce back into solid-state lithium.The rate of reduction toward the center is controlled by the mass transport of lithium through the crystal.This is shown by the inset in Fig. 8a, which provides the simulated relaxation profiles of the mobile ions in the agglomerate.A similar result is observed for the 32 nm simulations.
In addition, the agglomerate distributions provide information about the utilization of the active mass.For both simulations, at the end of discharge/start of recovery (t = 0), the reaction was only able to penetrate ∼20% of the crystal radius.Using this value, the percent volume of the agglomerate which participated in the discharge, ν active , can be calculated using the following: where r agg is the radius of the agglomerate (1.05 μm for these simulations) and r is the penetration depth of the reaction into the agglomerate (0.21 μm).The results of this calculation indicate that only 48.8% of the active material in the agglomerate was utilized during discharge.
In contrast to the agglomerate-scale distributions, there are strong differences between the distributions of solid-state lithium within the crystal at the surface of the agglomerate.For the 8 nm simulations (Fig. 9a), there is little spatial variation within the crystal, which indicates that mass transport within the crystal has a negligible impact on the voltage recovery.The decreases in concentration over time are due to the relaxation behavior on the agglomerate scale.For the 32 nm data (Fig. 9b), large variations within the concentration of the solid-state lithium are observed.At the start of relaxation, the concentration at the surface of the crystal is over 3× higher than the concentration at the center.These distributions reinforce the conclusion that crystalscale effects only impact the voltage recovery for experiments with the large, 32 nm crystals and not those conducted with the 6 or 8 nm crystals.It also suggests that the discrepancy between the simulations and experiments at 32 nm (Fig. 7b) is because other phenomena such as phase change within the crystal may contribute to the large recovery times.
In addition to phase change, other factors may explain the discrepancy in the voltage recovery times for the cases with the 32 nm crystals.In the following two paragraphs, three such factors are discussed which were investigated by making adjustments to the multi-scale Figure 10.Experimental and simulated agglomerate distributions for electrodes comprised of 8 nm crystals.A similar distribution was observed for electrodes with 32 nm crystals. 15del.For brevity (and because none of the factors showed significant improvements in the predicted results), no simulated data are shown.First, a decrease in the agglomerate scale diffusion coefficient due to variations in the geometry and packing of the nano-crystals may partially explain the discrepancy in voltage recovery times.To test this hypothesis, simulations were conducted with a decreased agglomerate scale diffusion coefficient (2× lower or 1.15 × 10 −13 cm 2 s −1 ).Slight improvements were observed in the agreement between the rise time of the simulations and experiments.However, significant discrepancies between the maximum voltage changes during recovery ( V max ) were observed.In an attempt to get better agreement with V max , the diffusion coefficient in the crystal scale (D x ) was also adjusted for the 32 nm simulations.However, even over a wide range of D x (1.0 × 10 −15 to 1.0 × 10 −18 cm 2 s −1 ), good agreement between the experiments and simulations for both the voltage recovery time and V max could not be obtained.Based on these findings, we concluded that a decrease in the agglomerate scale diffusion coefficient is not likely to be the sole reason for the discrepancies in rise time.
In addition, we tested the hypothesis that the presence of electrochemically inactive crystals could explain the discrepancy in the voltage recovery time.Simulations were conducted assuming 60% and 85% of the crystals were electrochemically active.Both sets of simulations yielded negligible improvements in the agreement between the rise time and negatively impacted the agreement with V max .Therefore, we also concluded that inactive crystals would not likely explain the discrepancy.Finally, we investigated the use of concentrated solution theory to describe the mass transport within the crystal, whereby the differences in chemical potential provide the driving force for mass transport.To accomplish this, the governing equation for mass conservation within the crystal (see Table I) was reformulated as shown in Refs.24-27.A new solid-state diffusion coefficient was obtained by fitting the reformulated set of equations to experimental data using the procedure outlined earlier in the text.Slight variations in the predicted performance were observed with concentrated solution theory; however, the agreement with experimental data was not improved.
Impact of agglomerate distributions.-Theanalysis of the voltage recovery data using the multi-scale model indicated that the formation of agglomerates impacts the electrochemical performance of the magnetite electrodes.These simulations were conducted using a single, average agglomerate size; however, electrodes typically contain a range of agglomerate sizes. 15To study the impact of agglomerate size distributions on the predicted performance, the multi-scale model was adjusted to include a representative distribution of agglomerate sizes.
Figure 10 shows the experimental distribution of agglomerates composed of 8 nm crystals obtained from transmission electron micrographs of the cross-section of a fully-fabricated battery electrode.A similar distribution for agglomerates composed of 32 nm crystals was also observed. 15Along with the experimental data, Figure 10 includes the representative distribution of agglomerates used in the simulation.The sizes (i.e., diameters) and number fractions (i.e., frequencies) were selected to best match the experimental data.These values were incorporated into a multi-agglomerate model, which is capable of simultaneously solving for the coupled concentration and potential distributions in all three representative agglomerate sizes.The full mathematical formulation of the multi-agglomerate model is available in Appendix A.
In order to understand the impact of the agglomerate distribution, two case studies were conducted, which compared the simulated results from a multi-agglomerate and a single-agglomerate model.For the single-agglomerate model, the agglomerate diameter was set equal to the average agglomerate size reported earlier in the manuscript (2.1 μm).In the first study, separate diffusion coefficients were obtained for the single-agglomerate and multi-agglomerate simulations by fitting both of the models to experimental data for a discharge at 4.7 mA g −1 to a cutoff of 1.5 V.The results of this study are shown in Fig. 11a.To get good agreement for both models, the multi-agglomerate diffusion coefficient (D m-agg ) is ∼2× higher than the single-agglomerate diffusion coefficient (D s-agg ).This indicates that, for the magnetite electrodes, a failure to incorporate the agglomerate distribution in the multi-scale model impacts the fitted agglomerate diffusion coefficient by a factor of ∼2.
In the second study, the diffusion coefficient in the agglomerate (D agg ) was obtained by fitting the multi-agglomerate model to experimental data.This value was then used in the single-agglomerate simulation.Figure 11b shows the results of this study.For the singleagglomerate simulation, a ∼50% increase in capacity is observed when compared to the experimental and multi-agglomerate simulation results.The increase in capacity can be explained by an increase in the utilization of the active mass.For instance, Figure 12 shows the predicted concentration distributions of solid-state lithium in the agglomerate for the single-agglomerate and multi-agglomerate models during the simulations in Fig. 11b.For all four agglomerates (one from the single-agglomerate model and three from the multi-agglomerate model), the simulated distributions indicate that the reactions only occur near the surface of the agglomerate.This is due to mass transport limitations of lithium ions through the agglomerate.
In Figure 12, the percent volume of each agglomerate which participates in the reaction, ν active , was calculated using Eq. 2. For the single agglomerate model, ν active is representative of the total active mass utilization in the electrode.For the multi-agglomerate model, the total utilization can be determined using the following equation.[3]   where f k is the number fraction of agglomerates with radius k.Evaluation of Eq. 3 results in a total active mass utilization of 47.3% for the multi-agglomerate simulation.The utilization predicted in the singleagglomerate model (71.9%) represents a 52% increase in the total utilization when compared to the multi-scale model.This accounts for the ∼50% increase in capacity predicted in Fig. 11b.The low utilization predicted by the multi-agglomerate model results from the low utilization of the large agglomerates in the distribution.Although the largest agglomerates have a number fraction of less than 10% (f 3 = 7.5%), they account for over 50% of the active mass.Therefore, even a few large agglomerates in the electrode can negatively impact the capacity of the battery.

Conclusions
We report here an analysis of the ion transport limitations occurring within a lithium-magnetite electrochemical cell using a combined experimental and theoretical approach.A multi-scale mathematical model, which accounted for mass transport in the agglomerate and crystal length-scales, was used to analyze experimental discharge and voltage recovery data.It was concluded that the long voltage recovery times of the magnetite electrodes were caused by the relaxation of concentration distributions, which developed as a result of mass transport limitations within the electrode.For electrodes comprised of 6 and 8 nm crystals, the mass transport limitations were shown to mostly occur within the agglomerate length-scale.For electrodes composed of 32 nm crystals, mass transport limitations were shown to occur in both the agglomerate and crystal length-scales.Therefore, the observed decrease in the discharge capacity between 8 and 32 nm was attributed to the addition of crystal-scale transport limitations.
In addition, the impact of a representative agglomerate size distribution on simulation results was studied using an expanded version of the multi-scale model.Inclusion of a representative agglomerate distribution indicated that variations in agglomerate size could impact the values of the fitted diffusion coefficients by a factor of ∼2.The inclusion of a small number fraction of large agglomerates was shown to significantly decrease the predicted capacity, which indicates a possible direction for improving magnetite electrode performance.

Appendix A: Multi-Agglomerate Model
The multi-agglomerate model was developed using the same assumptions and governing equations reported in Ref. 32.The agglomerate size distribution was simulated using three representative agglomerate sizes.The lithium concentrations (c agg,k and c x,k ) and voltage (ϕ 1,k ) distributions in the three agglomerates were solved for simultaneously by defining the variables using dimensionless groups and applying coupled boundary conditions between the agglomerates.To accomplish this, the following dimensionless groups were employed: where the subscript k, denotes the 1 st , 2 nd , or 3 rd agglomerate size in the representative distribution.The reaction rate, i rxn,k , is defined using the dimensionless Butler-Volmer equation: The equilibrium potential, U k , is obtained by fitting a thermodynamic description of the system to experimental voltage recovery data. 32Definitions for all other variables can be found in Appendix B.
The dimensionless groups in Eq.A1 are used to transform the original governing equations into dimensionless equations, resulting in the following expressions for mass and charge conservation in each agglomerate, k: Charge in agglomerates : [A5] The multi-agglomerate simulations were only conducted for the experiments with 8 nm crystals.Therefore, to decrease the solving time, it was assumed that the concentration of solid-state lithium in each crystal increased uniformly (i.e, no spatial variations of lithium within the crystals).The validity of this assumption is demonstrated by the predicted concentration distributions in Fig. 8a.This assumption makes it possible to solve for the solid-state concentration without solving for the mass transport in the crystal.Instead, the conservation of mass for the solid-state lithium was determined using the following equation: The dimensionless form of Eq.A6 is written as follows: Mass in crystals : The solution of these equations (Eqs.A2-A5 and A7) was obtained using the following boundary conditions: In Eq.A10, the potentials of all three agglomerates are set equal because it is assumed there are no spatial variations of potential within the bulk electrode.This assumption is valid for the small applied current (C/200) and thin electrodes (50 μm thick) used in the present experiments.The final conservation of charge boundary condition is obtained by setting the sum of the electronic current at the surface of all the agglomerates equal to the applied current.where f k is the number fraction of agglomerates of size k in the electrode and both sides of Eq.A11 have units of amps.
Note that the value of the dimensionless distance, r (Eq.A1), does not depend on the agglomerate size k.Because of this, it is possible to simultaneously solve the domain equations in all three agglomerate sizes using the same finite-difference grid.To accomplish this, the finite-difference method is used to discretize all nine governing equations in dimensionless space and dimensionless time (Eqs.A4, A5 and A7 for all three agglomerate sizes).At each time step in real time, t, the resulting block, tri-diagonal matrix containing all nine independent, dimensionless variables is solved using the BAND(J) algorithm. 42pendix B: List of Symbols a specific surface area (cm 2 cm −3 ) c agg lithium concentration in the agglomerate (mol cm −3 ) c 0 bulk concentration in the electrolyte (mol cm −3 ) c x solid-state lithium concentration (mol cm −3 ) c x,max maximum solid-state lithium concentration (mol cm −3 ) D agg diffusion coefficient in the agglomerate (cm

Figure 1 .
Figure 1.Schematic of the transport processes occurring on the crystal and agglomerate length scales, which provide the foundation for the crystal-only, agglomerate-only, and multi-scale models.

Figure 2 .
Figure 2. Voltage recovery data for Fe 3 O 4 electrodes comprised of crystals with an average diameter of a) 8 nm and b) 32 nm.

Figure 3 .
Figure 3.Comparison of experimental recovery time to agglomerate models with different Li + diffusion coefficients.Best fit diffusion coefficient was used for all agglomerate model simulations.
90% of its maximum value.As an example, the results of the fitting procedure for the agglomerate model are shown in Figure3.The figure compares the experimental recovery time to simulations using the best fit diffusion coefficient (D agg ), 0.5 × D agg , 2 × D agg , and 10 × D agg .The results indicate that fitting the simulations to the recovery time can provide a reasonable estimate for the diffusion coefficient, at least within an order of magnitude.The simulations using 0.5 × D agg could not be completed past 2.05 electron equivalents due to mass transport limitations (i.e., onset of a mass-transfer limited current was predicted).

Figure 4 .
Figure 4. Range of void space sizes expected for an agglomerate of closepacked nano-crystals.Size of the void space is determined from the diameter of the largest sphere capable of fitting in the void.Octahedral and trigonal packed spheres provide the upper and lower bounds of the agglomerate void space, respectively.

Figure 5 .
Figure 5. Maximum voltage during voltage recovery.Comparison of a) crystal-only, b) agglomerate-only, and c) multi-scale models to experimental data.Multiple simulations were conducted for each crystal size to account for slight variations in the experimental current density.The maximum and minimum V max are displayed as error bars.Current densities: 6 nm (4.4 to 4.8 mA g −1 ), 8 nm (4.4 to 4.8 mA g −1 ), 32 nm (4.5 to 4.9 mA g −1 ).

Figure 7 .
Figure 7.Comparison of experimental voltage recovery to the multi-scale model for electrodes composed of a) 8 and b) 32 nm magnetite crystals.

Figure 8 .
Figure 8. Distribution of the average intercalated lithium in each crystal (i.e., c x,avg ) throughout the agglomerate during voltage recovery after discharge to 1.5 electron equivalents per mole Fe 3 O 4 .Plots are for simulations of electrodes with a) 8 and b) 32 nm crystals.Inset in a) shows the distribution of lithiumions in the agglomerate (c agg ) during the same recovery.Similar predictions of c agg are observed for the simulations with 32 nm crystals.Symbols: r is the radial position within the agglomerate, r agg , is the radius of the agglomerate, and c 0 is the concentration of lithium-ions in the bulk electrolyte.

Figure 9 .
Figure 9. Distribution of lithium in the crystal (i.e., c x (r )) located at the surface of the agglomerate during voltage recovery after discharge to 1.5 electron equivalents per mole Fe 3 O 4 .Plots are for simulations of electrodes with a) 8 and b) 32 nm crystals.Inset in b) shows profiles at early times for simulations with 32 nm crystals.Symbols: r and r are the radial positions within the agglomerate and crystal, respectively, and r agg and r x are the radii of the agglomerate and crystal, respectively.

Figure 11 .
Figure 11.Comparison of simulations conducted using a single, average agglomerate size (single agglomerate) or a representative distribution of three agglomerates (multi-agglomerate). a) Diffusion coefficient is adjusted for both simulations to fit data.b) Simulations are conducted with the same diffusion coefficient.The electrode in the experiment was 50 μm thick with an active mass loading of 4.0 mg cm −2 and was discharged at 4.7 mA g −1 to a cutoff voltage of 1.5 V.

Figure 12 .
Figure 12.Simulated distributions of the solid-state lithium (c x ) throughout the agglomerates during discharge of the a) single-agglomerate and b-d) multiagglomerate simulations.Capacities are in reference to Figure 11b.

Table II . Diffusion coefficients used to fit models to experimental voltage recovery for electrodes with 6 nm crystals. Diffusion Coefficient (cm 2 s −1 )
2 s −1 ) D x diffusion coefficient in the crystal (cm −2 s −1 )