Analysis of Li-Ion Battery Electrochemical Impedance Spectroscopy Data: An Easy-to-Implement Approach for Physics-Based Parameter Estimation Using an Open-Source Tool

The quantitative analysis of electrochemical impedance spectroscopy (EIS) data is important for both characterization and prognostic applications in many electrochemical systems. Here we describe an open-source platform, the ImpedanceAnalyzer, for easy-to-use physics-based analysis of experimental EIS spectra. To demonstrate the use of the platform, we explore the basic capabilities of the pseudo two-dimensional (P2D) battery model to predict publicly available experimental EIS data from a 1500 mAh commercial lithium-ion (LiCoO 2 /graphite) cell. An a priori computed dataset of 38,800 P2D-based impedance spectra simulations, covering a wide range of frequencies (1 mHz to 100 kHz) and model parameters, enables a straightforward least squares matching approach for analyzing experimental spectra. We ﬁnd an average error of 1.73% between the best-matching computed spectrum from the 38,800 member library and the experimental spectrum being analyzed. Our analysis shows there is signiﬁcant opportunity to improve the ﬁt between experimental data and physics-based impedance simulations by a combination of a larger computed dataset, local optimization, and further additions to the model physics. The approach and open source tools developed here can be easily extended to other electrochemical systems.

Electrochemical impedance spectroscopy (EIS) is a powerful tool for investigating a wide variety of electrochemical systems. 1-3 EIS spectra separate individual electrochemical processes by their characteristic timescales, enabling both qualitative and quantitative analysis of electron transport, 4,5 reaction rates and mechanisms, 6,7 intercalation processes, 8 mass transport, 9,10 and electrode structure. 11,12 The noninvasive nature of EIS also makes impedance measurements useful in prognostic applications such as fuel cell health estimations 13,14 or prediction of remaining useful lifetime in batteries. 15,16 Qualitative analysis of EIS spectra generally involves assessing the shape of Nyquist plot features to determine the relative importance of different physicochemical processes. 17,18 In contrast, quantitative analysis relies on fitting a model to the data in order extract values for specific thermodynamic, transport, and/or kinetic parameters. Most experimental datasets are analyzed quantitatively using an equivalent circuit analog. Fitting an equivalent circuit to EIS data is straightforward using standard least squares regression techniques. 19,20 A good fit can often be found with a relatively simple equivalent circuit, particularly if non-ideal elements like the constant phase element are used. Moreover, many simple equivalent circuits, like the Randles' circuit, 21 have physically interpretable parameters based on linearized electrochemical processes. However, as more complex equivalent circuits are derived and utilized, the lumped parameters can lose their direct physical interpretability and the structure of the equivalent circuit analogs themselves can be degenerate. 22 An alternative to equivalent circuits for quantitative analysis of EIS data is to directly fit the data with a physics-based mathematical model of the electrochemical system. Many years of electrochemical modeling research have laid the groundwork for the physics-based analysis of impedance in a wide variety of fields including corrosion, 6 hydrodynamic systems, 23,24 fuel cells, 25,26 and lithium-ion batteries. [27][28][29] Parameter estimation by fitting a physics-based model to experimental EIS data is complicated by the fact that electrochemical models often contain a combination of coupled differential equations, algebraic equations, and dozens of unknown parameters. For complex physics-based models, convergence to a global best fit is rarely assured, even with an excellent initial guess for the unknown parameters. As a result, parameter estimation methods often need to rely on many independent measurements to drastically narrow the number of fitted parameters. 30 In short, today there is gap between the desire to use physics-based models to estimate parameters in EIS, and the actual (routine) use of physics-based models for parameter estimation from data.
Here we demonstrate an easily implemented and extendable approach for leveraging sophisticated physics-based models of EIS spectra to estimate parameters from experimental data. The parameter estimation approach used in this work relies on error minimization between experimental data and a large library of a priori simulated impedance spectra. One benefit of a dataset-based approach is that it always converges and the resulting parameter estimates are guaranteed to be reasonable if the original dataset was constructed from physically reasonable parameters. Another benefit of a dataset-based approach is that global sensitivity analysis can be used to understand the contribution of different parameters to the model variance. 31 On the other hand, a dataset-based approach will generally not provide the parameters with the lowest possible error between model and experiments.
The work and software tool presented here is extendable to a wide range of electrochemical systems, though our first implementation is using the Doyle-Fuller-Newman pseudo two-dimensional (P2D) lithium-ion battery (LIB) model as the basis for analyzing EIS experiments. 32,33 The earliest uses of LIB physics-based impedance models have been to inform the analysis of experimental EIS data. For example, Doyle et al. 28 used physics-based EIS simulations to show that the low frequency portion of a LIB impedance spectrum appears (qualitatively) to be interpretable as a Warburg impedance, but using a Warburg plot can produce erroneous (quantitative) diffusivity estimates. Subsequently, additional physics has been added to the original P2D model to aid in the interpretation of EIS data. For example, a surface oxide model was added to the positive electrode particles by Dees et al. 30 to understand the increase in interfacial impedance with aging of LiNi 0.8 Co 0.15 Al 0.05 O 2 -based (NCA) positive electrodes. Abraham et al. 34 extended the model further to interpret the changing impedance response at different voltages in NCA electrodes.
Despite the widespread use of the P2D model for simulating LIBs (2 of the top 4 most cited papers in the history of the Journal of the Electrochemical Society (as of November 1 st , 2017), 35 the model has seen limited use for quantitative analysis of experimental impedance data. In this work, we present the ImpedanceAnalyzer, an open-source, web-based analysis platform aimed at making physics-based models as easy to use as equivalent circuits for quantitative analysis of EIS experimental data. For the first implementation, we explore the basic capabilities of the original P2D impedance model to predict publicly available experimental EIS data from a commercial cell. We discuss the implications of this work in terms of the next steps for expanding the model dataset, extending the models physics (surface oxide layers, surface-electrolyte interphase (SEI) layers, etc.), and adding gradient-based parameter estimation to the basic Sobol' parameter sampling method 36 that forms the initial backbone of the ImpedanceAnalyzer. As an open source tool, any of these modifications can be added by other researchers to improve performance. Overall, we believe the open source tools and approach presented here have the capability to combine knowledge generated from decades of physics-based impedance modeling research with the deep materials and chemistry insight of experimentalists to accelerate progress in the field.

The pseudo two-dimensional (P2D) lithium ion battery model.-
We have previously detailed the isothermal P2D model and efficient frequency domain computational approach applied here to prepare a library of simulated lithium-ion battery EIS spectra. 37 Briefly, the P2D model is a set of partial differential equations describing the one-dimensional, volume-averaged distribution of lithium ions and potential in the solid and solution phases across a positive electrode, separator, and negative electrode cell sandwich. In each of the porous electrodes, lithium ion intercalation is assumed to be governed by Fickian diffusion into spherical particles coupled to faradaic charge transfer via Butler-Volmer kinetics. The non-faradaic capacitive currents assume a simple Helmholtz double-layer model. The full set of governing equations, boundary conditions, and parameter dependencies are transformed into the frequency domain by assuming a steady periodic solution form for all dependent variables driven by a single-frequency sinusoidal modulation of applied current. We apply a Volterra series to capture the amplitude-dependence for the linear and weakly nonlinear higher harmonics that result from the singlefrequency input current perturbation. Further information showing the model equations used for EIS spectra here are presented in the Appendix AI of our prior work. 37 The generic Coefficient Form PDE physics module in COMSOL v4.4 was used here for all computations. As noted in earlier work, 37 the highest density of nodes were placed at the interfaces between the solid and solution phases and at the electrode/separator interfaces, based on results from preliminary analysis. For all results presented here, the nodes are distributed in a geometric sequence with the node spacing at the electrode/separator interface 25X smaller than at the center of the separator or at the current collector interface. All 38,800 computed spectra in our EIS library had the same mesh with 250 nodes across the cell sandwich and 8000 elements in the particles. Because of the size of the computed EIS library, we did not use mesh refinement to validate that every individual spectrum was converged, as is normal when a small number of simulations (compared to 38,800) are used. We discuss our strategy for validation of the computational results in §2. 3. For a given set of physicochemical and geometric parameters and single input perturbation frequency, computing the impedance response took approximately 4 seconds on a Dell Precision T1500 with an Intel Core i7 CPU @ 2.80 GHz and 8GB RAM using a Windows 7 Professional 64-bit operating system. Parameters were updated and the results were saved to disk using the LiveLink for MATLAB version R2013b. The 38,800 member library of EIS spectra (with each spectrum having 25 frequencies) required approximately 2 CPU-months to compute, though much of it could be done in parallel.
Creating a large dataset of simulated spectra.-The P2D impedance model as implemented here is parameterized with 26 physicochemical and geometric coefficients. To adequately capture the wide variety of impedance spectra one might encounter for the diversity of lithium-ion battery chemistries, a dataset of different simulations was synthesized using a wide, but physically meaningful, range of parameters. The parameters and their ranges are shown in Table AI in the Appendix. Parameters whose ranges spanned more than two orders of magnitude were sampled logarithmically. Sampling the parameter space of the several dozen inputs to the P2D model necessitates a large number of simulations. The Sobol' sampling sequence 36 was used to efficiently explore the high dimensional space and the SALib 38 and savvy 39 python packages were used to generate the sampling sequences.
An important, but subtle, issue in establishing model parameters is that impedance measurements represent a "local" probe of the battery at the state-of-charge being assessed. The most relevant case for understanding the meaning of "local" is the treatment of electrode thermodynamics. The governing equations for the linear electrochemical impedance (or the weakly nonlinear electrochemical impedance, for that matter 37 ) do not rely on knowledge of "global" relationship between open circuit voltage and intercalated lithium concentration. Instead, what matters is the local gradient. Thus, Table AI shows the gradient of open circuit potential with intercalated concentration as a parameter influencing the impedance. The local open circuit voltage gradient values for each electrode are constrained by realistic bounds found in the literature for multiple chemistries. As local variables, the thermodynamic gradients in each electrode are taken as independent parameters to be sampled like all others in Table AI.

Incorporating physical knowledge into the validation of many computations.-
The dataset used in this work currently contains 38,800 impedance spectra each containing 25 frequencies. As the number of computations grows to cover a wide range of parameter and frequency space, deterministic methods for validating model accuracy such as adaptive mesh refinement or manually verifying the number of node points in high gradient regions become impractical. Leveraging the physical understanding of the electrochemical system provides a crucial insight into exploring the tradeoff between computational time and accuracy.
As an example, in any porous electrode, the competing effects of the solid-and solution-phase conductivities and reaction kinetics determine the distribution of current density within the system. 40 That is, for parameter sets with fast kinetics and low effective conductivities, the zone over which the lithium flux occurs is small, and high gradients form. The location of the reaction zone is driven by the ratio of solid-and solution-phase conductivities. Consequently, for the lithium-ion battery system, where electrode conductivities are typically higher than electrolyte conductivity, the region of significant lithium flux density in the porous electrodes tends to be shifted toward the electrode/separator interfaces. Similar considerations occur in the diffusive flux of the oscillating lithium ion concentrations within the solid electrode particles. The penetration depth is governed by Fickian mass transport in the system, where D s is the solid-phase diffusivity and ω is the perturbation frequency. The boundary layer thickness decreases for lower mass transport coefficients and higher frequencies. To account for these regions of high gradients in the concentrations, potentials, and current densities in the system, a high number of node points at the interface between the solid and solution phases and the electrode/separator interfaces were used. To establish confidence in the large EIS dataset we have produced, we explore the most challenging computational limit (high frequencies) and compare it with the easily calculable analytic limits of infinite frequency. In the limit of infinite frequency, the only contributor to the impedance response is the combined ohmic drop across each of the electrodes and separator. The internal resistance of the cell has the analytical solution, [2] where l i is the thickness of the region, and σ ef f,i and κ ef f,i are the effective solid-and solution-phase conductivities (given by Brugg where i and f, i are the void-and filler-fractions and Brugg is a Bruggeman-type tortuosity factor, respectively). Eq. 2 sets a theoretical lower bound for the computed electrochemical impedance for any set of parameters. Comparisons of high frequency computed solutions to Eq. 2 for the same parameter set is one way we will assess the numerical accuracy of our large computed dataset over the wide parameter range used.
Fitting computational spectra to experimental data.-To directly compare a computed EIS spectrum to an experimental EIS spectrum from a battery with unknown electrode area requires a self-consistent method to determine the superficial electrode area, A m sup , for every computed spectrum m (here, our computed library has 1 ≤ m ≤ 38,800). Superficial area is needed because the computations provide area normalized impedances, with units in − m 2 , whereas EIS experimental measurements have units in . For self-consistency, the area used to normalize EIS simulation results must produce a cell capacity that matches the experimentally determined capacity of the battery under test, C data (in Ah), which is a universally measured or reported value. Thus, for every computed spectrum m, the basic procedure is to calculate each electrode's capacity, C m pos and C m neg , according to where the subscript index i is either pos or neg, A m sup is the unknown superficial area (in m 2 ), and V i is the volumetric capacity (in Ah/m 3 ) of the material in each electrode. The self-consistent superficial area is the only unknown when the capacity-limiting electrode from each computed spectrum m is equated to the experimental capacity data by One last factor often found in experimental data, but absent from computed spectra, is contact resistances (R m contact in ). Uncompensated contact resistance leads to a shift along the real axis that will produce poor or erroneous fits to experiments. Thus, we use contact resistance as a free parameter to minimize the residual error E m between the experimental data and an area-scaled spectrum from the computed EIS dataset, where R m contact is the error-minimizing contact resistance for spectrum m, Z m (ω j ) and Z m (ω j ) are the real and imaginary components of the computed impedance values for spectrum m at frequencies ω j , Z data (ω j ) and Z data (ω j ) are the experimentally measured real and imaginary impedance data at frequencies ω j , and N is the total number of frequencies being fit. No weighting of the data with frequency is currently implemented.
Because R m contact is the only fit parameter to minimize error, it can be aphysical. We eliminate any spectrum from consideration if the best-fitting R m contact is negative and appreciable (here appreciable is taken to mean a magnitude greater than 10% of the high frequency limit). The residual error E m for all simulated spectra can then be ranked from lowest to highest error to determine the top matches for an experimental spectrum.

Results
Flexibility of physics-based modeling.-Due to the large ranges over which the physicochemical and geometric parameters were sampled, the resulting dataset of simulated impedance responses contains spectra of a wide variety of shapes in the Nyquist diagram. Figure 1 shows some of the variation among the 38,800 computed spectra demonstrating the flexibility of the P2D model to capture distinct impedance responses. Figures 1a and 1b show the Nyquist representation of the simulated spectra with one overlapping (Fig. 1a) or two separable (Fig. 1b) kinetic arcs. The arcs are related to the interfacial charge transfer resistance and double-layer capacitance at the two electrodes. When the two kinetic impedance responses are similar in magnitude, but differ in characteristic time constants, separate arcs (or a flat long arc) are visible. When either a single electrode has extremely facile kinetics (low charge transfer resistance) or the electrodes have similar kinetic time constants, a single arc defines the mid-frequency impedance response. Consequently, a wide range of high-to mid-frequency responses are seen in the simulated dataset without the requirement on introducing unphysical or nonideal circuit elements like the constant phase element. Similarly, despite the typical treatment of fitting a single Warburg element, the low-frequency response also encompasses a diverse range of features. Figure 1c shows an impedance spectrum with a small low-frequency response while Figure 1d shows an impedance spectra with a comparatively large contribution in the same frequency range. In the low-frequency regime, solid-and solution-phase mass transport as well as thermodynamics (derivatives in the open circuit potential) interact to generate distinct variations in the impedance responses. As an example, Figure 1e and Figure 1f show the effect of an impedance spectrum with a relatively flat low-frequency feature and a steeply sloped feature, respectively. It should be noted that there is no need for the spectra in Figure 1 to be scaled by superficial area or shifted by inclusion of uncompensated contact resistance since we are not comparing to experimental results, and therefore, the figure units are − m 2 .
Assessing the quality of the simulation library.-Preliminary studies were carried out to determine a meshing strategy that worked for a wide range of parameters and frequencies. The goal of the preliminary work was to balance the number of mesh nodes and their spatial distribution, with the computational time required to compute a converged solution. Given the wide range of parameters used, and the 38,800 unique parameter combinations explored at 25 frequencies (970,000 computations of the governing equations and boundary conditions), it was not easy to evaluate the accuracy of every computation. Nonetheless, we had a strategy to ensure the numerical approach and meshing was appropriate. In particular, we explored the hardest-to-converge limit of the governing equations, high frequencies, by comparing the simulated high-frequency data to the analytical high-frequency limit given by Eq. 2. Figure 2 shows a comparison of the simulated high-frequency real impedance at 10 5 Hz and the analytically predicted value of the ohmic resistance. An initial indicator of computational quality is seen immediately; Eq. 2 is the theoretical lower bound for the real impedance, and we see that all simulations lie on or above the diagonal line, with none below. There are two possible reasons for the simulated impedance at 10 5 Hz to be above the analytically-predicted ohmic (real) resistance: (i) the simulated frequency (10 5 Hz) is insufficient to reach a purely ohmic response for the set of parameters or (ii) numerical errors associated with inadequate meshing. Understanding the origin of the upward deviations we see in Fig. 2 at low cell resistances is aided by coloring the markers to indicate the computed phase of the calculated impedance at 10 5 Hz for all 38,800 simulations. We see a consistent trend of higher phase angles resulting in greater deviation from the diagonal line which is consistent with the deviation being associated with 10 5 Hz as too low of a frequency to reach the real axis. We can explore this effect statistically. The average deviation from the predicted value for all 38,800 high frequency points is 12.9%, while the Figure 2. Comparison between the 38,800 simulated data points for the highest frequency used here, Z (10 5 Hz), and the analytic ohmic resistance for the same parameters, R , predicted . The marker color indicates the simulated phase angle at 10 5 Hz. Inset shows that some combinations of parameters require a much higher simulated frequency to achieve a purely ohmic response. error for all spectra with less than a 5.7 • phase angle at 10 5 Hz (68% of the spectra) is 5.0%. In short, the closer a point at 10 5 Hz is to a purely real number, the more accurately Eq. 2 predicts the value. This may seem obvious, but there is no reason to believe that inadequate mesh refinement would produce this systematic behavior across such a diverse set of parameters.
To further test the idea that the meshing and computations are adequate, and deviations arise predominantly from the selection of 10 5 Hz as the highest frequency, we further evaluated some of the largest deviating points. The inset in Fig. 2 shows the decrease in error resulting from increasing the highest simulated frequency from 100 kHz to 100 MHz, while keeping the grid mesh fixed. We see that all three points move toward the diagonal line, further validating that it is primarily the frequency range used, for certain parameter combinations, not numerical error.
Based on these high frequency and other tests, we have confidence that the systematic deviations seen in Figure 2 are not primarily due to numerical error. Moreover, given that we primarily tested the most difficult-to-compute frequency limit of the library, we are confident that meshing-produced numerical errors are at least an order of magnitude below the 5.0% systematic deviation seen for the low phase angle points in Fig. 2, for all frequencies. Of course, adaptive methods for increasing the density of node points within the high gradient regions would lower the numerical error further, albeit at the significant cost of increased computational time, but there appears little reason to pursue that based on these results. Figure 3. To fit an experimental spectrum, a user would browse for a file on their computer containing a comma delimited file where the first column contains frequency, the second contains the real impedance, and the third contains the imaginary impedance. The type of analysis to be performed is then selected and any additional input (battery capacity for scaling according to Eq. 4, for example) is entered before clicking the "Go" button to transmit the data to the server for the spectrum matching and ranking process described above.

Using the ImpedanceAnalyzer: An example physics-based fit of a LiCoO 2 /graphite cell spectra.-The input panel for the web-based, ImpedanceAnalyzer tool is shown in
To demonstrate the process, an example experimental impedance spectrum is taken from the University of Maryland's Center for Advanced Life Cycle Engineering (CALCE) Battery Data Archive. 41 The spectrum shown here comes from the initialization impedance measurements on a 1500mAh LiCoO 2 /graphite battery cell (Cell 41 in the PLN Initialization Dataset). 42 The impedance response of the pouch cell was measured from 1.64 kHz to 12.5 mHz. The Nyquist representation of the raw experimental data is shown in Figure 4a. The spectrum consists of a depressed semicircle in the frequency range associated with interfacial processes and a low-frequency tail in the region associated with mass transport and thermodynamic processes. As a part of the spectrum ranking process, the experimentally sampled frequencies are quadratically interpolated to match the simulated frequencies ( Figure 4b) and then the process described above determines the simulated spectra from the dataset with the lowest residual error (Figure 4c). The resulting closest match for this example spectra and the P2D dataset described above has a run index, m, of 6230, a superficial area, A sup , of 300.3 cm 2 , a contact resistance, R contact , of 4.21 m , and an average residual error, E, of 1.73%. In this case, the positive electrode provides the limiting capacity. The parameters which generate the closest match are returned to the user in the web tool and are shown in Table AII in the Appendix. It should be noted that the experimental dataset used here has the data truncated at a moderate frequency and, thus, the spectrum contains negligible influence from the inductance often found experimentally at high frequencies.
For datasets that contain higher frequencies, there may need to be a simple inductive element added to the contact resistance parameter in Eq. 5 to capture these experimental artifacts.
Exploring the remaining solutions.-One of the benefits of the dataset-based approach taken in this work is that the nearby "nearly- matching" solutions can also be explored. Quickly visualizing the many spectra that are close to matching the experimental spectrum enables an experimentalist to begin piecing together a better understanding of the interacting physicochemical processes in the complex electrochemical system. Figure 5 shows the residual errors calculated using Eq. 5 for the top 50 simulated spectra as a function of the spectra's rank. Multiple results with the same residual error occur when the only change in the input parameter set is a parameter that the linear impedance response is insensitive to (such as the charge transfer coefficients, α a and α c , for example, where only their sum appears in the governing equations and boundary conditions). The steepness of the drop-off at low rankings (best fitting spectra) indicates that there is still potentially significant benefit to increasing the number of simulations in the dataset.  To further demonstrate the different variations by which the simulated spectra fail to match the experimental data and how we can learn from them, the three spectra in Figure 6 show the starred 3 rd , 7 th , and 10 th closest matching spectra in the dataset from Figure 5. For example, while Figure 6a (m = 34560) matches well in the lowfrequency region, it fails at matching in the high-frequency regime. Conversely, the spectrum in Figure 6b (m = 4932) qualitatively matches, but has the timescales for kinetics shifted, while the spectra shown in Figure 6c (m = 32014) has a qualitatively similar response in the high-frequency regime, but fails at low frequencies.
Interactively exploring these nearby solutions is made easy in the ImpedanceAnalyzer via an "Explore P2D" modal where a user can mouse over the residual points and interactively see the corresponding spectrum and parameter values. Clicking on a residual point allows the comparison of multiple spectra and a download of the selected parameters is made available. The three parameter sets for the spectra in Fig. 6 are presented in Table AIII. We see that the qualitative improvement in the low-frequency response of the 3 rd best spectra is likely due to the decreasing solid-phase diffusion coefficient in the negative electrode, D s,neg , while further decreasing the diffusion coefficients many more orders of magnitude leads to the very poor match at low-frequencies for the 10 th best spectrum. Additionally, decreasing the magnitudes of the open circuit potential derivatives leads to a smaller low-frequency tail between the 3 rd and 7 th best matching spectra. For the mid-frequency response, the lower negative exchange current density, i 0,neg , leads to similar time constants between electrodes resulting in a single, narrower arc for the 3 rd ranked spectrum with respect to the best match. Decreased double-layer capacitances (increasing the characteristic frequency of both electrodes) broadens the mid-frequency response of the 7 th ranked spectrum into a single, flatter kinetic arc, while increasing the positive electrode capacitance (lowering the characteristic frequency) and decreasing the negative electrode capacitance (increasing the characteristic frequency) removes the separation between the two electrode arcs in the 10 th best fit.
The limiting electrode (i.e. the electrode which has the computed area of 1500 mAh according to Eq. 4) in each of the simulations shown in Figure 6 is the positive electrode except for the 10 th best Figure 7. Nyquist plot comparing the experimental data with a two parameter locally optimized P2D simulated spectrum. A 20% increase in the double layer capacitance, C dl,neg , and a 50% decrease in the solid diffusion coefficient, D s,neg , resulted in a spectrum with a 1.25% error (33% lower than the best matching spectra in the 38,800 member dataset). match. Interestingly, due to the symmetry of the P2D model, it is possible to switch the labels of pos and neg parameters and compute an identical spectrum. Additionally, the estimated superficial area for each of these nearly matching spectra vary more than 100 cm 2 . Consequently, including additional information about the battery (including any known physicochemical or geometric ranges of parameters) will be useful for improving the matching process by filtering the dataset. Adding this feature is currently under development.

Matched spectra as an initial guess for local optimization.-Sam-
pling across a large range of physical parameters enables a fit such as that shown in Figure 4c, but, because we used parameters associated with a wide range of different battery chemistries, states-of-charge, etc., the 38,800 member EIS library is not a dense covering of the parameter space. Moreover, for any finite number of simulations, the global nature of the sampling process makes it likely that the resulting match is not a best-fit to the experimental spectrum. However, if there is low experimental noise, the selected spectrum makes a good starting point for further local optimization of the physics-based model. In the case of the spectrum shown in Figure 4c, the largest contributions to the remaining 1.73% error are in the low-frequency tail as well as the kinetic arc at high frequencies. Using conventional understanding of impedance spectra, the width of the high frequency arc can be attributed to the charge transfer resistance associated with the intercalation kinetics in both electrodes, while the time constants of the arcs (dictating their separation and shape) are also influenced by the double-layer capacitances. Thus, to change the shape of the kinetic arcs without altering the total width, the double-layer capacitances can be varied. Additionally, because the time constants for kinetics and mass transport are separated by more than an order of magnitude, changes in either parameter group only affect the highfrequency arc and low-frequency tail, respectively. Figure 7 shows the consequence of a simple iterative reduction of residual error by varying the double-layer capacitances as well as the solid-phase diffusion coefficients and open circuit potential derivatives. For a 20% increase in double-layer capacitance, C dl,neg , and a 50% decrease in the negative electrode solid diffusion coefficient, D s,neg , the residual error decreased to 1.25%. Directly incorporating sophisticated local optimization in the vicinity of several highly ranked solutions found by global sampling is a functionality currently being explored for the ImpedanceAnalyzer.

Implications and Concluding Remarks
The physics-based nature of the P2D and other electrochemical models enables a wide set of physically relevant and interpretable impedance responses to be captured by simply varying the combinations of physicochemical parameter inputs. With a physics-based approach, it is transparent what phenomena are included in the model and what assumptions underpin the model. In contrast, equivalent circuits are reduced order models of (normally) complex physical and chemical processes, sometimes leading to ambiguous or incorrect interpretation of results. 28 While it is normally possible to fit some form of an equivalent circuit to an experimental spectrum, linking the fit to the underlying phenomena for anything but the simplest electrochemical systems remains a fraught exercise.
In the work described here, the P2D battery model was used to demonstrate the dataset-based approach for parameter estimation. As noted in the introduction, the high citation rate of the P2D model makes it the de facto "standard" continuum model for capturing the dynamics of the lithium-ion system. Moreover, significant research efforts have focused on adding additional physics to the system making it a good candidate for demonstration. That said, the open-source, dataset-based approach described here enables an opportunity to start statistically comparing electrochemical models as different physical processes are included or removed. Combining the sampling scheme with methods for performing a global, variance-based sensitivity analysis 43 of the simulated dataset can provide information on the parameters most responsible for variations in the impedance spectra at a given frequency. It has been shown that parameter identifiability is an important consideration in fitting electrochemical models to experimental data. 44 Comparing the parameter sensitivity as more complex physical interactions are sequentially added to models can provide statistical insight into the tradeoff between a more descriptive physical representation of the system and potentially worse parameter identifiability as additional parameters are added. Additionally, comparing an experimental spectrum to the nearest matching spectra within a single dataset can enable a bootstrap approach to establishing confidence intervals in the estimated parameters.
Several additional benefits arise from the dataset-based approach described here. The computational time is spent upfront during the generation of the dataset, significantly reducing the time an experimentalist waits to find a best matching spectrum. The time needed to identify the best matching fit to experimental data is unchanged by the complexity of the model. Only the size of the library affects the estimation time. Of course, a more complex model will take longer to compute the library of spectra; however, separating the model complexity from the time required to estimate the parameters becomes increasingly important as the quest for deeper insight into complex electrochemical systems grows. Furthermore, the open-source nature of the platform enables the inclusion of additional electrochemical data such as higher order harmonic responses, 37 spectra at multiple states (depths-of-discharge or temperatures), as well as (dis)charge curves or cyclic voltammetry.
The power of open-source software like this lies in the opportunity for subsequent contributors to add features, additional simulated or experimental data sets, etc., thereby evolving the tool to a more useful form. To that end, the ImpedanceAnalyzer is simply the start of a platform to which others are encouraged to contribute. The code is openly available on GitHub and users of the software can cite the version released at the time of this manuscript. 45 Other open software and open data products being actively built by the electrochemical data science community are also available for those interested in collaborative development. 46 Acknowledgments Partial support of this research is provided by the U.S. Department of Education grant P200A120023, an NSF IGERT grant DGE-1258485, and the Boeing Endowment for Excellence. Data and the analysis code for this paper can be found online. 47