EIS study on metal hydride electrodes using a porous model: Fitting methodology and SOC effects.

Journal of Energy Storage 27 (2020) 101067

Contents lists available at ScienceDirect

Journal of Energy Storage journal homepage: www.elsevier.com/locate/est

EIS study on metal hydride electrodes using a porous model: Fitting methodology and SOC effects. Marcos Martíneza, Erika Telizb, Carlos F. Zinolab, Verónica Díaza,

T

⁎

a University of the Republic, School of Engineering, Chemical Engineering Institute, Electrochemical Engineering Interdisciplinary Group, J.Herrera y Reissig 565, CP 11300, Montevideo, Uruguay b University of the Republic, School of Science, Fundamental Electrochemistry Lab, Electrochemical Engineering Interdisciplinary Group, Igua 4225, CP 11400, Montevideo, Uruguay

A R T I C LE I N FO

A B S T R A C T

Keywords: EIS Hydrides Batteries Electrochemical energy storage

In order to understand and improve the performance of Metal Hydride (MH) materials for Ni-MH battery's applications, Electrochemical Impedance Spectroscopy is used to provide interfacial information. However, it still offers significant challenges regarding methodology practices and results' interpretation due to the multiple interpretations arising from equivalent circuit fittings. The use of physicochemical models expresses in a better way the functionalities of charge/discharge mechanisms in MH. Moreover, for the impedance of MH porous electrodes there is little information from previous models, so the physical meaning of parametric figures and time constants is here emphasized. A complete fitting methodology is detailed and implemented via a specifically designed Octave program. An AB2 alloy of Zr0.3Ti0.7Cr0.7Mo0.3Ni composition is used to construct a porous electrode. Impedance spectra of the activated electrode are recorded and standard deviations of data are experimentally estimated. Validity of both the model and the methodology are thoroughly assessed. Changes in the system response when varying its state of charge are fully analyzed.

1. Introduction It is well known that in our modern societies the worldwide energy consumption keeps a clear growing tendency [1,2]. On top of that, fossil fuels from finite sources exhibit negative impacts on the environment; so many alternatives to those feedings are being examined. In this context, research and development in both generation and storage of renewable energy acquires the highest importance [3–8]. In light of this situation, investigation on metal hydrides (alloys capable of store hydrogen and form a metal hydride phase) is of great relevance [9–16]. These materials are studied to use them in nickel-metal hydrides (Ni-MH) batteries (electrochemical hydrogen storage) and for gas phase hydrogen storage. It is worthwhile noticing that other more innovative applications are also being considered [17,18]. Ni-MH batteries became a consolidated technology for diverse mobile applications or low power consumption devices [19–21], showing interesting characteristics regarding technical requirements, costs and environmental impacts. Besides that, MH constitutes one of the safest choices for storing gaseous hydrogen [22–25]. However, with the aim of spreading this technology out to other applications with higher performances, continuous research for materials with better mechanical

⁎

and electrochemical properties is mandatory. Electrochemical Impedance Spectroscopy (EIS) is an electrochemical technique which provides valuable information about interfacial properties of dynamic systems, such as in MH electrodes. EIS has already been used to study this type of electrodes through the application and development of different models and methodologies [26–30]. The basis of EIS is to measure the response of a system when it is perturbed with a low amplitude sinusoidal potential signal (potentiostatic EIS) or current (galvanostatic EIS) with continue change in the frequency in order to obtain an impedance spectrum. One of its biggest advantages is that it provides qualitative and quantitative information regarding pathway's contributions to the involved mechanism. However, it still offers various experimental and methodological challenges, many of them regarding interpretation of the results. When the technique is not used in the best way, there is a loss on quality and quantity of the obtained information. If it is not well understood, it can easily lead to misinterpretation of the results [31,32]. This work is related –continuing the line of research of our Group– to the study of the composition and electrochemical properties of Zrbased alloys (AB2 metal hydrides) in order to understand their behavior

Corresponding author. E-mail address: verodiaz@fing.edu.uy (V. Díaz).

https://doi.org/10.1016/j.est.2019.101067 Received 18 September 2019; Received in revised form 6 November 2019; Accepted 6 November 2019 2352-152X/ © 2019 Elsevier Ltd. All rights reserved.

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al.

especially its frequency dependence. Electrode activation was performed in a four-channel cycler device. Potentiostatic EIS measurements were performed using a Gamry Interface 1000™ Potentiostat/Galvanostat.

and then, to improve their electrochemical energy storage performance [12–16]. With the aim of improving the quality of the analysis of EIS results, it was sought to conduct a new work methodology including a formal support without losing physical meanings. In this respect, EIS data are able to be compared with results obtained from other studies on the same alloys. Besides consolidating the work methodology, it was decided to develop a model which takes into account the porosity of the electrode. One of the aims was to study the frequency range where porous effects are negligible, and other zones where they are imperative. While many models already exist that describe porous intercalation electrodes [33–35] (and the developed model is based on them), it was considered as obligatory to advance a model that clearly enables to fit the data with simple, intuitive and physically meaningful parameters. A fitting program was created using Octave programming language in order to test the different choices in fitting methodology and to customize the used model in a well understood manner. The adequacy of the selected model and fitting choices could be assessed. Changes in the system response when varying the electrode State Of Charge (SOC) are also studied, since it was found that this parameter was relevant to EIS results. In order to focus the attention into this methodology and modeling, it was decided to narrow the study to a single alloy, (Zr0.3Ti0.7Cr0.7Mo0.3Ni, named ``Ti7′') which was characterized in previous studies with useful kinetic properties for electrochemical hydrogen storage [16].

2.2. Physicochemical model A physicochemical model of the porous electrode is proposed here to obtain the equations describing the impedance of the system. Metal hydride electrodes store energy through the electrochemical adsorption and subsequent absorption of hydrogen into metallic particles. The charging process on the alloy begins with hydrogen adsorption on the metal surface (Volmer reaction). In alkaline medium, it can be depicted as follows:

H2 O + e− + Ms ⇌ MHads + OH−

(1)

where Ms is a free adsorption site at the surface. Under certain conditions, adsorbed hydrogen may give rise to Hydrogen Evolution Reaction (HER), following either Heyrovsky or Tafel mechanisms;

MHads + H2 O + e− ⇌ H2 + Ms + OH−

(2)

2MHads ⇌ 2Ms + H2

(3)

In the charged or discharged alloy, adsorbed hydrogen has to follow absorption and diffusion steps:

2. Experimental Section

MHads + Mss ⇌ Ms + MHabs, R0

(4)

MHabs, R0 ⇌ MHabs, r

(5)

where Mss is a free interstitial subsurface site, MHabs,R0 is an absorbed hydrogen atom in r = R0 (surface of the alloy particle), and MHabs,r is an absorbed hydrogen atom in r < R0 (bulk of the particle). Inside the bulk of the particle, nucleation and growth of one or various hydride phases are known to happen. Alloys to be studied with this model are complex multiphasic alloys, but in order to simplify the approach, the bulk of the particle is assumed to be composed of a single phase with a given hydrogen saturation concentration CH,max and a fixed apparent diffusion coefficient DH. The electrode under study consists of alloy powder particles, mixed and compacted with a conductive carbon powder. As a result, a porous conductive matrix is formed, in which the solution remains in intimate contact with alloy particles. To account for the porous nature of the electrode we introduce volumetric expressions for impedances (Ω cm3). The volumetric interfacial impedance Zi can be expressed by the following equation:

2.1. Experimental Design A Zr0.3Ti0.7Cr0.7Mo0.3Ni (called Ti7) alloy was synthesized by arc melting under high purity argon from pure metallic ingots of each component. The resulting alloy ``button'' was re-melted twice to ensure the sample homogeneity and it was mechanically pulverized and sieved for the electrode formation. The latter was prepared by mixing and compacting 100 mg of sample powders –with an initial particle size of 63–125 μm– with equal amounts of powder teflonized carbon (Vulcan XC-72). The mixture was pressed under 250-MPa at room temperature and a nickel wire was used as current collector. The electrochemical single compartment cell was mounted placing the working electrode (metal hydride electrode), counter electrode (nickel mesh) and reference electrode (Hg/HgO electrodes). All the experiments were carried out at room temperature and all potentials referred to the Hg/HgO electrode. The electrolyte, 6 M KOH solution, was prepared from reagent-grade potassium hydroxide and MilliporeR MilliQ Plus water. The electrode was cycled by a charge–discharge process at constant density currents of 80 mA g−1 and 26 mA g−1, respectively, with a charge time of 5 h and cut-off potential of −0.60 V during discharge. Potentiostatic EIS measurements were performed once the alloy had reached its maximum capacity after the activation. Prior to each EIS experiment, the electrode was completely charged and partially discharged, until reaching the desired SOC, then left at open circuit potential (EOCP) during 90 min in order to minimize the system drift. EIS spectra were recorded, at EOCP, in the 5.00E +04 – 1.00E-04 Hz frequency range, with a 5 mV amplitude, ten points per decade. After this, the electrode was completely discharged (until reaching a cut-off potential of −0.60 V) and the SOC was recalculated, considering both discharge times (after and before recording the EIS spectrum) and discharge capacity of the electrode. EIS spectra were recorded at SOCs (recalculated) of 92%, 77%, 73%, 62%, 57%, 46% and 30%. Multiple consecutive impedance spectra (n = 6) of the same electrode were recorded at SOC = 73%, with the objective of determining the standard deviation behavior for a given state of the system,

Zi−1 = Zdl−1 + ZF −1

(6)

where Zdl is the volumetric impedance of the double-layer and ZF is the volumetric faradaic impedance. Both of these elements arise from surface processes. Consequently, parameters that relate area to volume ratios (cm2/cm3) must be considered. 2.2.1. Faradaic impedance Faradaic impedance is related to the electron transference occurring at the surface of metallic particles. The faradaic impedance per unit of interfacial active area Zf (Ω cm2) follows the following relation (demonstrated and detailed in references [30] and [32]):

Zf = R ct +

1 jωCp +

1 R ev

+

1 R ab + Zw

(7)

where Rct is the charge transfer resistance (related to the Volmer step), Cp is the adsorption pseudocapacitance (related to the surface capability of storing adsorbed hydrogen), Rev is a HER resistance (related to the Heyrovsky step), Rab is an absorption resistance and Zw is the Warburg impedance (related to the hydrogen diffusion inside the alloy). Fig. 1 depicts the equivalent circuit for the faradaic impedance Under certain conditions, this equation can be simplified: 2

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al.

Fig. 1. Equivalent circuit for the faradaic impedance according to Eq. (7).

- If the Tafel reaction rate is negligible and experimental potentials are too low for the Heyrovsky reaction to take place, i.e. there is no HER altogether and Rev is taken as infinite. - If the absorption step rate is fast (fast enough so that it can be considered to be always in equilibrium), then the existence of Cp in parallel to Rab + Zw will not modify the impedance of the system [30]. - If both of the above conditions are met, then Rct and Rab become indistinguishable. It has been stated that most of times Rab is very small [32]. Usually, the whole resistance is reported as Rct [30,36–39].

Fig. 2. (Left) Diagram displaying electrode geometry and dimensions. (Right) Diagram indicating the direction in which the current density iT is considered to flow when iT > 0.

2.2.3. Structure and impedance of the porous electrode Structurally, the electrode is a cylinder of width 2 L with a current collector at the center (x = 0, between -L and L). Since the diameter is much larger than the width (D>>2 L), electron flow will be considered in only one dimension (x). The electrode porous matrix (alloy and carbon particles) is filled with electrolyte. Since both phases are conductive, the total current density iT, per geometric electrode area (A cm−2) is the sum of iL and iS, which are the current densities in the liquid and solid phases (see Fig. 2 The charge transference between phases at any point will be determined by the volumetric interfacial impedance Zi. Near the current collector (x = 0) all the current will flow through the solid phase (iT = iS) and at the end of the electrode (x = L) all the current will flow through the liquid phase (iT = iL). A potential distribution is encountered in both phases, related to the phase conductivities (κ and σ for liquid and solid phases, respectively). The following set of equations describe the potential distribution in solid (ΦS) and liquid (ΦL) phases, as previously stated by Castro et al. [34] (following earlier works [33,41]):

Thus, assuming these conditions, the faradaic impedance is reduced to:

Zf = R ct + Z w

(8)

The Warburg impedance expression depends upon the assumed geometry. For identical spherical particles of radium R0, it can be written as:

Zw =

⎡ σ′ ⎢ DH ⎢ jω coth ⎢ ⎣

1

(

R0 DH

jω

)−

DH R0

⎤ ⎥ ⎥ ⎥ ⎦

(9)

where σ’ (in Ω cm3 s−1) is a parameter that, under the assumptions considered beforehand, can be figured as:

∂v1 ⎞ σ ′ = R ct ⎛⎜− ⎟ ∂ ⎝ CH ,0 ⎠

dϕS (10)

iS σ

(13)

dϕL i =− L dx κ

(14)

iT = iL + iS

(15)

dx

where v1 is the Volmer reaction rate and CH,0 is the alloy hydrogen concentration at r = R0. The expression displayed in Eq. (10) was achieved by simplifying a previous one for σ’, demonstrated in reference [32]. Details for the simplification are exhibited in Appendix A.

=−

ϕM − ϕL = −Zi 2.2.2. Double layer impedance The double-layer impedance is usually modeled with a CPE (Constant Phase Element):

ZCPE =

1 T (jω) ϕ

Zp =

where T (in units of F sΦ−1 cm−2 or F sΦ−1) is a parameter related to the capacitance of the double-layer and Φ is the constant phase exponent (0 < Φ < 1). Since the CPE behavior is considered to be originated from capacitance dispersions with the frequency, distinct models to estimate an average double-layer capacitance have been proposed. The one reported by Brug et al. [40] is used here:

ET iT

(17)

where ET is the potential difference between current collector and electrolyte (ET = (ΦS)x=0 − (ΦL)x=L). The system can be simplified by considering that the conductivity of the solid phase is infinite (σ → ∞). This approximation may not be close to reality, but considering both conductivities they give rise to models whose solutions can be approximated in this way verified by results of simpler models [32]. Besides, using a more complex model would not add any more information of interest to the results. Thus, ΦS will not be a function of x and ET simplifies to

1− 1 ϕ

⎜

(16)

By solving these equations, along with the boundary conditions at x = 0 and x = L (mentioned beforehand), one can find another one for the electrode impedance per unit of geometric area Zp (Ω cm2):

(11)

1 1 1 ⎞ Cdl = T ϕ ⎛ + R ct ⎠ ⎝ Rs

diS di = Zi L dx dx

⎟

(12) 3

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al.

τx(s) is a time constant related to hydrogen absorption rates, defined

ET = (ΦS − ΦL)x=L== Zi (diL/dx)x=L. From the above set of equations, it can be found that d2iL/dx2 = (1/Ziκ)iL. Since the boundary conditions are iL = 0 in x = 0 and iL = iT in x = L, the solution for iL is

iL =

iT iT [e λx − e−λx ] = sinh(λx ) (e λL − e−λL) sinh(λL)

R

(18)

2

where λ2 = 1/(Ziκ). This solution can be used to find expressions for ET and Zp:

di ET = (ϕS − ϕL) x = L = Zi ⎛ L ⎞ = Zi λ iT coth(λL) ⎝ dx ⎠ x = L

Zp =

ET = Zi iT

1 1 ⎞ L coth(Λ0.5) coth ⎜⎛ L⎟ = κZi Λ0.5 κ ⎝ κZi ⎠

⎡ DH τx = ⎢ ⎢ ∂v1 ⎢ − ∂CH ,0 ⎣

(

(19)

(20)

- the lower limit of the data frequency range is not low enough. - the particles are too big or the diffusional coefficient is too low. - there is not enough low frequency data (number of points per decade), or the data quality is not good.

2.2.4. Total impedance of the system The measured total impedance Ztot(Ω) is the sum of the electrode impedance and the solution resistance:

Ztot = Rsol +

AG

= Rsol

2.3. Fitting methodology

L

Fitting methods for parameter determination and evaluation of results are mentioned in this section. Methodology is more thoroughly explained in Appendix B. References [32] and [31] provide a complete summary of proper EIS fitting practices. A program was developed (using Octave programming language) specifically designed for fitting impedance data. Octave version 4.4.0 was used. Prior to performing EIS experiments, it was sought to minimize systematic errors arising from experimental design. However, it was found to be very difficult to completely avoid these errors for this particular case. Altought Kramers-Kronig relations are often quoted as being able to evaluate the significance of systematic errors [31,32,42,43], it was not possible to integrate them properly in the methodology. The chosen approach was to obtain more measurements to diminish systematic errors, but never to assume a fully random error structure. Complex non-linear least-squares (CNLS) method was used to fit the data, as it is the most accepted regression technique for this purpose [32]. The Marquardt-Levenberg algorithm was used for this purpose. It was implemented via Octave’s leasqr function, included in the optim package. Version 1.5.2 of this package was used. The employment of different weighting choices was evaluated (unit, modulus, proportional) through analyzing standard deviation of the data, as exhibited in the Results Section. Modulus weighting was selected and used for carrying out the fitting processes. Goodness of fit was assessed through visual inspection of complex plane and Bode plots, and through the calculation of the reduced sum of squares Sv. Values of Sv for the selected fits of the data were assessed through analyzing standard deviation of them. Confidence intervals for model parameters (estimated from the regression procedure) were calculated and their validity discussed. The intervals were calculated from the parameters covariance matrix, provided within the results of Octave’s leasqr function. Simplifications on the model (explained in the Physicochemical Model Subsection) were also considered through analyzing the data. In addition, F-test was employed to estimate at which frequency ranges the practice of a porous model is statistically justified.

(21)

where AG is the total geometric area of the electrode (both faces). The impedance of the porous electrode can be simplified at lower frequencies (ω → 0):

coth(Λ0.5) 1 1 ≈ + Λ0.5 3 Λ

(22)

So, the total impedance results in:

Ztot = Rsol +

L Zi + 3 κ AG L AG

Ztot = Rsol + Rpor +

Zi L AG

(25) R2

where Λ has been defined as Λ = L /(Ziκ). It is noticeable that this expression has a great similarity with that developed by De Levie for cylindrical pores [41].

coth(Λ0.5) + AG κ Λ0.5

)

⎤ ⎥ ⎥ ⎥ ⎦

Finally, τdif(s) is a diffusional time constant expressed as τdif = D 0 . A H low τdif indicates that hydrogen diffuses easier to the center of the alloy particles. Sometimes, for given sets of impedance data, it is difficult or not possible to find the value of τdif for many reasons, mainly because:

2

Zp

D

here as τx = [ ct σ ′ H ]2 . An alloy with high τx will be able to absorb higher quantities of hydrogen at low frequencies. Under the considered assumptions (following simplification of σ’ detailed in Appendix A) it can simply be expressed as:

(23)

(24)

where Rpor is the ohmic contribution of the porous electrode, which becomes constant when ω → 0. This means that at the lowest frequencies Rsol and Rpor are indistinguishable. In this sense, one of our objectives was the setting of the frequency range at which the above approximation becomes acceptable for the studied data. 2.2.5. Model parameters and time constants The final expression for the impedance allows the outcome of seven parameters: two related to ohmic resistances (Rsol and Rpor), another pair related to the double-layer behavior (TCPE and ΦCPE), one being the charge transfer resistance (Rct, related to the Volmer step and the subsequent absorption) and the other two connected with hydrogen diffusion (σ’/DH0.5 and R0/DH0.5). It is important to note that when fitting real data to the model, values obtained from some parameters (TCPE, Rct and σ’/DH0.5) cannot be fitted independently from the total interfacial active area of the electrode. Thus, without knowing the area, their relative figures per unit of area will remain uncertain. In this case, the area dependence should be taken in account to prevent erroneous conclusions when comparing the fitted magnitudes for different electrodes. Values for the parameters also depend on the working conditions of the electrode (for example, on the electrode SOC). Use of time constants is useful to compare results independently of the electrode area. In order to do so, three time constants are proposed: τct, τx and τdif. τct(s) is defined as: τct = R ct Cdl being the charge transfer time constant, which is linked to the electrochemical reaction rate. A low τct could indicate favorable kinetic parameters or a high density of active sites.

3. Results and discussion 3.1. Model and fitting procedure analysis In order to select between different methodologies and to assess the adequacy of the model, experimental results were analyzed for SOC 4

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al.

measure and the next). Thus, the real standard deviation of each measure would be lower. It was found that if 1.0E-04 – 1.0E-02 data are excluded, estimated ε2 value will be of 9.8E-05, which is lower than the found Sν values.

Table 1 Experimental results for average measured impedance values and calculated standard deviation values. σ'

σ''

Hz

Average Z values (n = 6) Z' Z'' Ω Ω

|Z| Ω

Ω

Ω

1.0E+04 1.0E+03 1.0E+02 1.0E+01 1.0E+00 1.0E-01 1.0E-02 1.0E-03 1.0E-04

0.61 0.65 0.70 0.80 1.04 2.61 5.58 6.15 7.15

0.61 0.65 0.70 0.82 1.14 3.22 5.67 6.19 7.98

0.0068 0.0074 0.0074 0.0074 0.0089 0.0241 0.1247 0.1219 0.2181

0.0008 0.0003 0.0011 0.0012 0.0044 0.0463 0.0383 0.0292 0.2904

Frequency

0.02 0.03 0.06 0.13 0.47 1.89 1.03 0.67 3.55

3.1.2. Fitting results Fig. 3 presents Nyquist plot for the studied data where many zones can be differentiated: - At very high frequencies, where Z’ is lower, an inductive behavior can be observed, attributable to equipment and wires inductance. Next, a shoulder can be observed or, in other cases, a small semicircle. It has been reported that the resistance of this semicircle is independent of the mass electrode and related to the contact between electrode and current collector [37–39]. Impedance of these zones cannot be predicted with the developed model. - At somewhat lower frequencies, a zone is observed where Z’ and -Z’’ increase at the same rate for decreasing frequencies. This is characteristic of the porous behavior predicted by the model. - At middle frequencies a semicircle is found, resulting from the coupling of the double layer capacitance with the charge transfer resistance. Impedance contribution from the hydrogen diffusing in the alloy is still small. - At low frequencies, -Z’’ reaches a minimum and starts increasing again due to limitations imposed by hydrogen diffusion (Z’ also increases).

73%. 3.1.1. Standard deviation analysis Multiple impedance spectra were consecutively recorded at the same SOC (as described in Section 2) to calculate the standard deviation of the data at different frequencies and to compare its behavior with the one predicted by different weighting methods. Experimental results are shown in Table 1. Since experimental standard deviation values differ by orders of magnitude depending on the frequency, the use of unit weighting was discarded. Proportional and modulus weighting were evaluated, since they are popular and easily implementable choices. Standard deviation functions σi = ε Zi were estimated for both cases using the Microsoft Excel simple linear regression tool. Estimated values of ε were of 0.029 for proportional and of 0.021 for modulus weighting choices. Table 2 shows estimated standard deviations calculated by proportional and modulus weighting methods. It can be observed that neither of the two options can accurately predict experimentally calculated standard deviations. Coefficient of determination r2 was of 0.64 for modulus and of 0.59 for proportional weighting. Modulus weighting has a slightly better r2 value, although it highly overestimates σ’’ values. It was decided to discard proportional weighting since it gave too much weight to low –Z’’ data points, often causing the algorithm to converge to solutions with high systematic errors in other zones. Therefore, modulus weighting was chosen. Sν parameter (used to assess goodness of fit) theoretically approximates to ε2 when the quality of a fit is good. For modulus weighting, a ε2 value of 4.4E-04 was estimated, therefore higher values would be expected for Sν. However, Sν values as low as 1.1E-04 were obtained while fitting data for the same electrode at the same SOC. Nevertheless, it seems reasonable that Sν and ε2 values are of the same order. The difference may not be so odd, as the model used for predicting standard deviation was proved to be not very accurate. Another possible explanation is that experimental data at lower frequencies have been more affected by the system drift (systematic changes in the low frequency zone were visually detected between one

The absence of another semicircle between the semicircle zone and the diffusional zone indicates that the absorption step is never the rate determining step. This means that the absorption process is much faster than the adsorption one. Consequently, neither Rab or Cp should explicitly appear in the equivalent circuit derived from the model, as it was explained in Section 2.2. Since the model has many (seven) parameters, it is usually difficult to fit the entire model straightforward using all the relevant frequency range. This only can be done with seed values close to the parameter figures that minimize Sν. If only data from an arbitrary frequency range in the semicircle zone is fitted, a good fit is almost always easily achieved using a resistance R1 in series with a coupling of a resistance R2 and a CPE. This simple model has only four parameters. Three other parameters should be added to achieve the full model: one associated with the porous behavior and two associated with the hydrogen diffusion inside the alloy (Warburg impedance). It was found that a good strategy involves a first addition of a Warburg impedance. Since the model still does not account for porous behavior, all the high frequency range where porosity has an effect should be excluded. It is of interest to be able to define the frequency after which the porous and non-porous models are equivalent. This possibility was studied using f-tests, as it is explained in the following subsection. The data was found to have a porous behavior for frequencies higher than 0.8 Hz. Finally, points in the mid-high frequency zone can be added (>0.8 Hz for the data used in the current work) using a porous model. Seed value for Rpor can be visually estimated: it is approximately equal to the width of the zone of porous behavior in the complex plane plot. As more points are included in the fit, values of other parameters (found by the algorithm) get more defined. For the studied data, points were gradually added in the higher end of the frequency range, while trying to avoid including the zone that the model cannot explain. Sν was approximately constant until reaching a frequency of 1.2E+02 Hz. It was found that immediately afterwards, both Sν and systematic errors increased. Table 3 exhibits the fitting results before and after adding data affected by porosity (in the 1.2E+02 – 8.0E-01 Hz range, where the porous model had to be used). It can be observed that the sum of Rsol

Table 2 Comparison between standard deviation values calculated from experimental results (σexp) and estimated through modulus (σmod) and proportional (σprop) weighting methods. frequency Hz

σ'exp Ω

σ'mod Ω

σ'prop Ω

σ''exp Ω

σ''mod Ω

σ''prop Ω

1.0E+02 1.0E+01 1.0E+00 1.0E-01 1.0E-02 1.0E-03 1.0E-04

0.007 0.007 0.009 0.024 0.125 0.122 0.218

0.015 0.017 0.024 0.068 0.119 0.130 0.168

0.020 0.023 0.030 0.076 0.162 0.178 0.207

0.001 0.001 0.004 0.046 0.038 0.029 0.290

0.015 0.017 0.024 0.068 0.119 0.130 0.168

0.002 0.004 0.014 0.055 0.030 0.019 0.103

5

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al.

Fig. 3. Complex plane plot for the studied data. Left Panel shows the full spectrum and Right Panel depicts only data of the highest frequencies (marked with a dashed square on the left panel).

3.1.3. Study of porous impedance behavior using f-test Experimental data in the semicircle zone were fitted to a porous and to a non-porous model, both excluding the Warburg impedance. The upper limit of the frequency range was gradually increased and each data was fitted to the two models. For every pair of fits, a f-test was conducted to determine whether the porosity was justified. It was observed that when fitting at the high frequency zone, the number of the considered points influences the result of the f-test. As it is expected, if a smaller number of high frequency points are used, the ftest will not reject the hypothesis (meaning that the simplified nonporous model should be used). When increasing the upper limit of the fitting frequency range, a point is reached from which the f-test starts rejecting the hypothesis. This means that when fitting data within that frequency domain, the porous model significantly improves the fit (with a confidence level of 0.95). The observed limit frequency was consistent with different seed values for the parameters. Slight variations were found while varying the lower limit of the frequency range. Reported values varied between 0.8 and 1.6 Hz. It was decided to use 0.8 Hz for this data as it indicated a porous behavior in a set with less experimental points.

and Rpor values in the porous model is similar to the value of Rsol in the non-porous model. The other parameters are similar in both models. Parameter values obtained by the algorithm exhibit little deviation when varying seed values. However, it is possible to hold one or more parameters at fixed values. This way, a fit with different characteristics is achieved. It was found that different fits can be attained while maintaining the same goodness of fit (with a similar Sν and with systematic errors appreciable in different zones). Confidence intervals calculated from the regression procedure were found to underestimate some of these variations. A common difficulty found when fitting low frequency data is that a low variation in φCPE parameter produces considerable variations in other parameters, while maintaining a similar goodness of fit. This makes it difficult to choose between different fitting results. Usually, the semicircle zone has low systematic errors and the choice between different fits involves a compromise between systematic errors in different regions within the low frequency range. A final decision may be based on the numeric values of Sν and visual comparison. Over similar fits, those with the lesser number of fixed numbers (if possible, no fixed data) should be chosen Table 4 shows the fitting results for the parameters, only excluding data of the higher frequencies (>1.2E+02 Hz). Calculated confidence intervals are also shown (pj ± 2 σj), and their size is displayed in relation to the parameter values (2σj / pj in%). For comparison, an arbitrary fit was chosen, where φCPE was fixed to be equal to 0.870. It can be observed that most of its values fall outside of the confidence intervals (percentual variation to the non-arbitrary fit values is also shown). However, Sν is similar for both fits: 0.00016 and 0.00011, respectively. Plots for experimental and fitted impedances are shown in Fig. 4. None of the fits have visibly greater systematic differences than the other.

3.2. Analysis of SOC effects on system parameters Multiple EIS spectra were recorded varying the SOC on the same electrode, with the objective of observing the changes on the parametric behavior. Spectra were recorded at SOC values of 92%, 77%, 73%, 62%, 57%, 46% and 30%. Data was fitted to the model for each SOC using data in the 1.0E+02 – 1.0E-04 Hz frequency range. Changes were observed for the model parameters and time constants.

Table 3 Model parameter fitting results for the studied electrode at a SOC of 73%, corresponding to a porous and a non-porous model (in different frequency ranges).

non porous (8.0E-01 – 1.0E-04 Hz) Porous (1.2E+02 – 1.0E-04 Hz)

Rsol (Ω)

Rpor (Ω)

Rct (Ω)

TCPE (F sφ−1)

ΦCPE –

σ'/(DH0.5) (Ω s−0,5)

R0/(DH0.5) (s0.5)

0.871 0.624

– 0.243

4.76 4.83

0.424 0.427

0.853 0.840

0.067 0.064

92 86

6

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al.

Table 4 Algorithm fitting results for model's parameters and their confidence intervals (calculated from the regression procedure). For comparison, results from an arbitrary fit are also shown, along with their percentual variation with respect to the algorithm result fit parameter values. Rsol (Ω) Algorithm result fit pj 0.624 pj+/− 2 σj 0.622 0.626 (2σj / pj) (%) 0.4% Arbitrary fit (ΦCPE fixed in 0.87) p'j 0.629 (Δpj)/pj +0.8%

Rpor (Ω)

0.243 0.236 2.7% 0.263 +8.4%

Rct (Ω)

0.249

4.83 4.74 1.9%

ΦCPE –

TCPE (F sφ−1)

0.427 0.423 1.0%

4.93

4.62 −4.4%

0.424 −0.9%

0.432

0.840 0.832 1.0% 0.870 +3.6%

σ'/(DH0.5) (Ω s−0,5)

0.848

0.064 0.057 10.3% 0.072 +13.4%

R0/(DH0.5) (s0.5)

0.070

86 72 16.0%

100

104 +20.9%

does not have a significant variation. The τx shows its highest value on SOC 57% and it decreases in both directions as the SOC moves away from this value. As R0/DH0.5, the τdif has a high uncertainty and does not show any clear trend, although it could exhibit a significant increase as the SOC approaches 100%.

3.2.1. Results Figs. 5 and 6 show Bode and complex plane plots for each of the spectra measured and their corresponding fits. Spectrum measured at a SOC of 62% is maintained in both figures to allow a better comparison between data of higher and lower SOC. First of all, it can be observed that the impedance modulus |Z| increases its value as the SOC decreases, for frequencies lower than 1.0E-01 Hz. Differences are more evident for SOC lower than 73%. For higher SOC, the variation is less pronounced, although the trend remains clear in the semicircle zone. Table 5 shows the fitting results for the model parameters. Their trends are plotted in Fig. 7 for comparing each of the parameters. It can be observed the tendency of the parameters to remain constant for every SOC, such as Rsol and Rpor, which are linked to resistances of the system and TCPE and ΦCPE, related to the electrochemical double layer behavior. Parameters Rct and σ’/DH0.5 show similarities in their behaviors, with a decreasing trend when the SOC increases and a stabilization of their values as they approach SOC 100%. The R0/DH0.5 parameter has a less clear trend, even more when this figure is taken with greater uncertainty, given that it depends on data points of lower frequencies. Table 6 shows the calculated time constants and their trends are plotted in Fig. 8. The τct presents the same trend as Rct, given that Cdl

3.2.2. Discussion of SOC effects on system parameters The analysis of the results can change and may have diverse interpretations depending on the model. Different types of models could be distinguished: - The one that consider two phases α and β that coexist in equilibrium. The hydride β phase is formed at a given saturation concentration of hydrogen from α phase. - The other that only consider the existence of only one phase (like the one used to fit the data). If the real system behaves like this kind of scheme, hydrogen concentration inside the alloy would be directly proportional to the SOC. - Others that may consider more complex situations. A simple condition of equilibrium between phases α and β could be identified by a well-defined plateau in PCT and ECT curves. When

Fig. 4. Complex plane and Bode diagrams comparing experimentally measured spectra, algorithm fitting results (with all parameters set free) and algorithm fitting results setting parameter φCPE fixed at 0.87. 7

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al.

Fig. 5. Complex plane and Bode diagrams comparing experimental data and fits for the studied electrode at different SOC between 30% and 62%.

Fig. 6. Complex plane and Bode diagrams comparing experimental data and fits, for the studied electrode at different SOC between 62% and 92%.

characteristics, but also on the total area of the formed double layer. Its invariance means that area did not significantly change, even with the volume variation of the alloy particles with the SOC. Rct value decreases when the SOC increases until reaching a constant figure between 73% and 92%. A stable behavior of Rct could be explained by a coexistence of phases α and β. In this way, the active surface would be in contact with α phase at a concentration CH,0 with no significant variations due to phases equilibrium. On the other hand,

working with AB2 alloys (as in the present work) the plateau visible in PCT and ECT curves is not as defined as in AB5 alloys. The AB2 alloy used here is a quite complex multiphasic nature and the real characteristics of the system may be very difficult to model in a precise way. It can be observed from the results that Rsol, Rpor, TCPE and ΦCPE do not vary with the electrode SOC. This could have been expected beforehand. It is reasonable that resistances related to electrolyte and solid remain constant. Values of TCPE depend on double layer 8

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al.

Table 5 Reported values for each of the model's parameters at different electrode SOCs. SOC (%)

Rsol (Ω)

Rpor (Ω)

Rct (Ω)

TCPE (F sφ−1)

ΦCPE –

σ'/(DH0.5) (Ω s−0.5)

R0/(DH0.5) (s0.5)

92 77 73 62 57 46 30

0.66 0.64 0.62 0.63 0.63 0.63 0.63

0.25 0.25 0.24 0.24 0.25 0.25 0.25

4.2 4.6 4.8 6.3 7.7 9.1 13.0

0.43 0.43 0.43 0.43 0.43 0.43 0.43

0.83 0.83 0.84 0.83 0.83 0.84 0.84

0.06 0.06 0.06 0.07 0.08 0.10 0.16

144 95 86 87 110 112 108

Table 6 Reported values for each of the model time constants and for the average Cdl varying the electrode SOCs.

the Rct reduction when SOC increases requires another explanation. A possible hypothesis is that the SOC increase induces changes in the surface structure (due to alloy volume modifications). This would enable new active sites in the process. Since the alloy is of a complex multiphasic nature, changes on the phases equilibrium inside the alloy may also explain changes in Rct (due to changes in CH,0). The parameter τx shows its higher value between SOC of 46% and 62%, decreasing when the SOC approaches 100%. A decrease in the value could indicate:

SOC (%)

Cdl (F)

τct (s)

τx (s)

τdif (s)

92 77 73 62 57 46 30

0.34 0.34 0.34 0.35 0.35 0.35 0.35

1.4 1.6 1.7 2.2 2.7 3.2 4.6

4216 6019 5719 7769 9147 7620 6728

20,689 9014 7393 7635 12,009 12,450 11,624

either an increase in R0 or a decrease in DH occurs). In a situation of α and β phases' coexistence, the diffusion radius should decrease as the β nucleus gets bigger. In a single phase model, the diffusion radius is independent of the SOC. Thus, the variation of τdif would not be explained by a change in R0 according to these models. On the other hand, R0 could grow with the SOC due to the increase in the alloy volume (although variation of R0 due to volume changes should be less pronounced at higher SOCs). An increase in τdif could also be related to a decline in DH values (which would also explain the decrease in τx).

- lower values of apparent hydrogen diffusion coefficient DH. - higher sensitivity of the reaction rate to changes in surface concentration.

4. Summary and conclusions Concerning the used fitting methodology, we arrived to the following conclusions:

If CH,0 is considered to increase at a higher SOC, then it is possible that the larger hydrogen concentration interferes with the hydrogen absorption and/or diffusion. This could be either due to a decrease on DH or to a blockage of the active sites. Although the data is not determinant, an increase in τdif may be occurring when the SOC approaches 100% (which would mean that

- The designed Octave program was useful to fit the experimental data. A methodology was designed to make the process easier and to avoid arbitrariness. - The reduced sum of squares Sν was a useful indicator to assess

Fig. 7. Parameter values obtained from the Model vs. SOC plots. 9

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al.

Fig. 8. Model time constants vs. SOC plots.

goodness of fit, and its values were of the same order of the estimated ε2 value for modulus weighting in this particular electrode (being ε the relative error of the measured impedances). - The used algorithm (Marquardt-Levenberg as implemented in Octave’s optim package) was not always able to minimize the value of Sν. Thus, changes in parameters should be manually tested (changing and fixing one parameter) to ensure no significantly better fit can be found. Between different fits with similar goodness of fit, the one found by the program should be chosen because it is free of arbitrariness (but it should be noted that it is affected by systematic preferences of the used algorithm). - Specifically, changes in ΦCPE significantly altered the other parameters values without worsen the goodness of fit. Fixed changes in ΦCPE may be useful to explore the possible variability of the other parameters. - It was found that the model parameter standard deviation values (σj) calculated from the parameter covariance matrix (resulting from the regression algorithm used) underestimate the variability of the parameters. Thus, it follows that the hypotheses made to calculate them were not adequate. Consequently, they should not be used to assess the variability of the parameters with a known confidence level.

- Since no second semicircle was found before the diffusive low frequency region, it follows that the hydrogen absorption step is much faster than the Volmer step. Consequently, the expression for the system impedance should be simplified accordingly (Cp and Rab must not explicitly appear). - A frequency limit can be defined, below which it does not matter whether the model takes into account the porosity of the electrode or not (expressions become equivalent). For the studied data, f-test was used to find that frequency limit (it was found to be approximately 1 Hz). - For the studied data, it was found that when excluding high frequency points (above 0.8 Hz) and using a non-porous model, the found parameters did not have excessive variations. Only Rsol changes drastically (Rsol value in the non-porous model is similar to Rsol+Rpor value for the porous model). The following conclusions can be made from the analysis of SOC effects on the system parameters: - Ohmic resistance related parameters (Rsol and Rpor) and double layer related parameters (τCPE and ΦCPE) did not show any dependence on SOC. Invariance in τCPE means that area of the double layer did not change significantly in spite of those of alloy particles' volumes. - A declining Rct for increasing SOC values may be related to changes on the surface structure as a consequence of the alloy particles' expansion, or to changes on the phases equilibrium inside the multiphase alloy. A more stable value toward SOC 100% is consistent with a multi-phase equilibrium. - Lower values of τx towards 100% SOC could be related to a more difficult hydrogen diffusion through the alloy (lower value of DH), or to a blockage of subsurface interstitial sites (which would increase the absorption rate sensibility to changes in hydrogen concentration). - An increase in τdif may exist as SOC approaches 100%, which could also be related to a minor DH value. Diffusion radium R0 may also vary, but it could hardly explain an abrupt increase in τdif at a high SOC.

Regarding the weighting choice for the fitting algorithm: - None of the selected popular weighting alternatives (unit, modulus, proportional) accurately envisaged the standard deviation behavior for the studied data. - Modulus weighting predicted the standard deviation behavior more accurately than unit weighting. Proportional weighting was discarded because it was found to cause difficulties when fitting, related to the fact that it gives too much weight to low imaginary impedance (-Z’’) data points. In relation to the model and its adequacy to the studied data: - The chosen model was adequate to calculate the system impedance in most of the measured frequency range. Above a certain frequency, the model cannot explain the system impedance behavior (that frequency was found to be approximately 1.0E+02 Hz for the studied data).

Declaration of Competing Interest None.

Appendix A: Simplification of the σ’ expression The equation describing the parameter σ’, related to hydrogen diffusion, is demonstrated in [32] (in Ω cm3 s−1 units):

10

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al. 2

( )( ) ( ) ( ) ( )

∂r 1 ⎡− F ∂r0 ⎤ ⎡− F ∂v4 ⎤ σ1 ∂θH ∂η BE ⎦ ⎣ σx ∂X0 ⎦ σ′ = − 2 = − ⎣ 2 AD ∂r 0 ⎡−F ⎤ ⎡− F ∂v4 ⎤ ∂η ⎣ ⎦ ⎣ σ1 ∂θH ⎦

(A.1) −2 −1

where v4 is the rate of hydrogen absorption (mol cm s ), r0 = v1 + v2 is the sum of Volmer and Heyrovsky reaction rates, r1 = v1 – v2 – 2v3 is the rate of change in the surface coverage by hydrogen through Volmer, Heyrovsky and Tafel reactions, η is the overpotential, θH is the fractional surface coverage by adsorbed hydrogen, X0 = CH,0/CH,max is the dimensionless hydrogen concentration at the surface, σ1 is the charge necessary to reach one monolayer coverage by adsorbed H (C cm−2) and σx = FCH,max is the charge corresponding to the metallic saturation with hydrogen (C cm−3). In absence of H2 evolution, r0 = r1 = v1. Additionally, if the absorption process is much faster than the adsorption and diffusion processes, then it can be considered to be always at steady state. Thus, for given conditions, X0 is a function only of θH and so:

( )( ) = ∂v ( ) ∂X ∂v1 ∂θH

∂v4 ∂X0

1

∂v4 ∂θH

0

(A.2)

This enables a great simplification of (A.1):

(− ) = R * F ⎛− ∂v ⎞ σ ⎝ ∂X ⎠ ⎡−F ( ) ⎤ ⎣ ⎦ F σx

σ′ =

∂v1 ∂X0

∂v1 ∂η

1

ct

⎜

x

⎟

0

(A.3)

where the Rct parameter is introduced (demonstrated in [32]). Appendix B: Fitting methodology details B.1. Error structure At each frequency, ωi, the difference between real impedance (Zi) and fitted model impedance values (Zi,calc) is composed of various types of errors (random and systematic). Systematic errors arise from various sources, during and after the conduction of the experiment. On one hand, error sources through experiences can be associated with experimental drifts (meaning that the system did not achieve a steady state prior to measuring the impedances). On the other hand, while proccessing the experimental data, defects in the chosen model and fitting methodology add also sources of systematic errors [31,32]. Various strategies are usually employed in order to minimize systematic error values. Noteworthy, it is often stated that Kramers-Kronig relations are able to indicate the importance of systematic errors in a given impedance set of data [31,32,42,43]. However, it was found very difficult to integrate the use of those relations in the developed methodology. In our particular case, it has been found to be very difficult to completly avoid systematic errors (for example, necessary waiting times for totally eliminating experimental drift are unacceptable, or giving rise to unwanted processes like the self-discharge of the electrode). B.2. CNLS method In order to fit the data to the model, the Complex Non-linear Least-Squares (CNLS) method was used. It can estimate a set of model parameters using real and imaginary impedance data simultaneously (Z’ and Z’’, respectively). It finds the values for the parameters that minimize a weighted sum of squares (S) of the differences between experimentally measured impedances (Z’i,meas and Z’’i,meas) and impedance values calculable by the model (Z’i,calc and Z’’i,calc) for each frequency i of the N frequency points taken per spectrum: N

S=

∑ {ω ′i [Z ′i,meas − Z ′i,calc]2 + ω′′i [Z′′i,meas − Z′′i,calc]2 }

(B.1)

i=1

where ω’I and ω’’I are the statistical weights at each frequency. To carry out the process, a minimization algorithm and a set of statistical weights should be chosen. It is also necessary to choose adequate values for the adjustable parameters. In this case, the Marquardt-Levenberg algorithm was chosen, which was implemented via Octave’s leasqr function, included in the optim package (Version 1.5.2 was used). While selecting the weighting method, it would be desirable that the set of statistical weights were an estimation of the data variances (ω’i = 1/ σ’i2 and ω’’i = 1/σ’’i2, for each frequency). In that case, the minimum value of S should approach the minimum value of the χ2 statistic [31,32]. Other desirable characteristic for the weighting method is that it should enable the algorithm properly fitting the data in all the frequency range. If too much importance is given to some data points, then the information in other frequency ranges may be poorly fit to the model. Unit, modulus and proportional weighting usage was evaluated during the present work, as defined in [32]. Statistical weights can be calculated through these simple methods using only experimental data. They cannot estimate the magnitude of the data variance, but the sets of parameters that minimize S for these cases may be close enough to the one that minimize the χ2 statistic. B.3. Goodness of fit assessment Algorithms used to carry out the CNLS method do not always give a satisfactory outcome. It is always very important to have proper criteria in order to choose whether a fit is satisfactory or not. Besides, a first visual inspection (which is always a primary tool to detect systematic deviations) is desirable to have a quantitative method to compare different fits for a given set of data points. If it is possible to assume that the errors were purely random and if the CNLS method is carried out using data variance as statistical weights, then 11

Journal of Energy Storage 27 (2020) 101067

M. Martínez, et al.

the minimized sum of squares Smin should be close to the minimum value of the χ2 statistic. The χ2min figure has a mean value of ν with a standard deviation of 2ν [31], where ν is the number of degrees of freedom (ν = 2 N – m, for 2 N measured impedance values (real and imaginary) and m adjustable parameters used in the model). The value χv2min = χ2min/ν should approach unity for a good fit. If proportional or modulus weighting is assumed, then the standard deviations will be supposed to be equal to σi = ε Zi, where ε is the relative error of the measured impedances (expected to be constant).This means that the value Sv,min = Smin/ν should be equal to ε2χv2min and Sv = S/ν should approach the value ε2 if the quality of fit is good. In practice, ε is not constant (it varies with the frequency) and it changes in orders of magnitude depending on the weighting choice. If ε is not known previously for a given system, then Sv cannot be used independently to evaluate the goodness of fit, but it can be used to compare different fits for a given set of data points (for which, the best fit will be the one with less value of Sv). Thus, a less value of Sv means that the goodness of fit is better. The value of Sv could also be used to compare the goodness of fits on different sets of data points, when ε can be assumed to be the same. If variances of the system (for each frequency) are independently estimated, then ε2 can also be assessed. Then, values of Sv can be compared with ε2 (variances should not be calculated from the regression procedure for this purpose [31]). B.4. Confidence intervals for parameters If errors related to the model parameter evaluation are assumed to be normally distributed, and their standard deviation σj is calculated, then confidence intervals for each parameter will be defined. If the assessed value for a parameter j is pj, then the probability that its real value lies within pj ± σj is of 67%. Probability that it lies within pj ± 2σj is of 95.4%. Standard deviations σj can be calculated from the regression procedure, although fitting errors are expected to be normally distributed. The regression also needs to be supposed to be linear in the region of χ2min Under these conditions, variances (σj2) of the parameters are the diagonal elements of the covariance matrix of each parameter. Commercial programs usually provide confidence intervals based on this method. Octave’s leasqr function provides the covariance matrix within its results (altough calculation details are not presented). In the present work, confindence intervals are constructed in this way and their vailidity is discussed. B.5. Assessing statistical impact of parameters When fitting real data to a theoretical model, it is possible that some of the model parameters have no appreciable statistical impact on experimental values. In these cases, approximations must be made in order to exclude those parameters and simplify the model. If reliable confidence intervals are calculated, that parameter will not have statistical significance when its confidence interval includes zero [31,32]. Another tool to assess the statistical significance of a parameter is the Fisher-Snedecor Test (F-test). It can be used when one has two fits for the same data using two models: one with m parameters and the other with m + k parameters. The F-test can exam the hypothesis about fitting the data with the simplest model (the one with minimum parameters) with a goodness of fit not significantly lower than when fitting to the more complex model. If this hypothesis is not rejected by the F-test, there will be no reason to accept the more complex model and the simplest one should be retained. If the hypothesis is rejected, then it can be considered (with a given confidence level) that the more complex model improves significantly the goodness of fit and it is necessary to select the model.

[24] J. Bellosta von Colbe, J.-.R. Ares, J. Barale, M. Baricco, C. Buckley, G. Capurso, N. Gallandat, D.M. Grant, M.N. Guzik, I. Jacob, E.H. Jensen, T. Jensen, J. Jepsen, T. Klassen, M.V. Lototskyy, K. Manickam, A. Montone, J. Puszkiel, S. Sartori, D.A. Sheppard, A. Stuart, G. Walker, C.J. Webb, H. Yang, V. Yartys, A. Züttel, M. Dornheim, Int. J. Hydrog. Energy 44 (2019) 7780. [25] M. Lototskyy, I. Tolj, Y. Klochko, M.W. Davids, D. Swanepoel, V. Linkov, Int. J. Hydrog. Energy (2019) S036031991931554X, in press. [26] E.B. Castro, D.J. Cuscueta, R.H. Milocco, A.A. Ghilarducci, H.R. Salva, Int. J. Hydrog. Energy 35 (2010) 5991. [27] S. Real, M. Ortiz, E.B. Castro, A. Visintin, Int. J. Hydrog. Energy 33 (2008) 3493. [28] Y. Zhang, Z. Hou, Y. Cai, H. Shang, Y. Qi, D. Zhao, J. Iron Steel Res. Int. 24 (2017) 296. [29] L.O. Valøen, S. Sunde, R. Tunold, J. Alloys Compd. (1997) 253–254 656. [30] L.O. Valøen, A. Lasia, J.O. Jensen, R. Tunold, Electrochim. Acta 47 (2002) 2871. [31] M.E. Orazem, B. Tribollet, Electrochemical Impedance Spectroscopy, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2008. [32] A. Lasia, Electrochemical Impedance Spectroscopy and Its Applications, Springer, New York, New York, NY, 2014. [33] J.P. Meyers, M. Doyle, R.M. Darling, J. Newman, J. Electrochem. Soc. 147 (2000) 2930. [34] E. Castro, S.G. Real, A. Bonesi, A. Visintin, W.E. Triaca, Electrochim. Acta 49 (2004) 3879. [35] A.M. Svensson, L.O. Valøen, R. Tunold, Electrochim. Acta 50 (2005) 2647. [36] K. Young, B. Reichman, M.A. Fetcenko, J. Alloys Compd. 580 (2013) S349. [37] N. Kuriyama, T. Sakai, H. Miyamura, I. Uehara, H. Ishikawa, T. Iwasaki, J. Alloys Compd. 202 (1993) 183. [38] A.A. Volodin, R.V. Denys, C. Wan, I.D. Wijayanti, S. Suwarno, B.P. Tarasov, V.E. Antonov, V.A. Yartys, J. Alloys Compd. 793 (2019) 564. [39] C. Wan, R.V. Denys, M. Lelis, D. Milčius, V.A. Yartys, J. Power Sources 418 (2019) 193. [40] G.J. Brug, A.L.G. van den Eeden, M. Sluyters-Rehbach, J.H. Sluyters, J. Electroanal. Chem. Interfacial Electrochem. 176 (1984) 275. [41] R. De Levie, Adv. Electrochem. Electrochem. Eng. Interscience, New York, 1967, p. 329. [42] Y. Fernández Pulido, C. Blanco, D. Anseán, V.M. García, F. Ferrero, M. Valledor, Measurement 106 (2017) 1. [43] M. Urquidi-Mac Donald, J. Electrochem. Soc. 133 (1986) 2018.

References [1] X.D. Wu, J.L. Guo, X. Ji, G.Q. Chen, Energy Policy 127 (2019) 287. [2] G.Q. Chen, X.D. Wu, J. Guo, J. Meng, C. Li, Energy Econ 81 (2019) 835. [3] Z. Yang, J. Zhang, M.C.W. Kintner-Meyer, X. Lu, D. Choi, J.P. Lemmon, J. Liu, Chem. Rev. 111 (2011) 3577. [4] E.J. Cairns, P. Albertus, Annu. Rev. Chem. Biomol. Eng. 1 (2010) 299. [5] P. Tan, H.R. Jiang, X.B. Zhu, L. An, C.Y. Jung, M.C. Wu, L. Shi, W. Shyy, T.S. Zhao, Appl. Energy 204 (2017) 780. [6] D. Guyonnet, M. Planchon, A. Rollat, V. Escalon, J. Tuduri, N. Charles, S. Vaxelaire, D. Dubois, H. Fargier, J. Clean. Prod. 107 (2015) 215. [7] K. Habib, H. Wenzel, J. Clean. Prod. 84 (2014) 348. [8] K. Shinyama, Y. Magari, K. Kumagae, H. Nakamura, T. Nohma, M. Takee, K. Ishiwa, J. Power Sources 141 (2005) 193. [9] T. Iwai, S. Takai, T. Yabutsuka, T. Yao, J. Alloys Compd. 772 (2019) 256. [10] K. Morimoto, I. Nagashima, M. Matsui, H. Maki, M. Mizuhata, J. Power Sources 388 (2018) 45. [11] F. Ruiz, E. Castro, S. Real, H. Peretti, A. Visintin, W. Triaca, Int. J. Hydrog. Energy 33 (2008) 3576. [12] E. Teliz, R. Faccio, F. Ruiz, C.F. Zinola, V. Díaz, J. Alloys Compd. 649 (2015) 267. [13] E. Teliz, J. Diez, R. Faccio, F. Ruiz, F. Zinola, V. Díaz, J. Alloys Compd. 744 (2018) 583. [14] E. Teliz, V. Díaz, F. Piganelli, R. Faccio, C.F. Zinola, J. Electrochem. Soc. 165 (2018) A3389. [15] E. Teliz, J. Diez, R. Faccio, E. German, F. Zinola, V. Díaz, J. Alloys Compd. 737 (2018) 530. [16] E. Teliz, J. Diez, M. Martínez, P. Díaz, F. Pignanelli, R. Faccio, C.F. Zinola, V. Díaz, JOM 71 (2019) 1952. [17] H. Uesato, H. Miyaoka, T. Ichikawa, Y. Kojima, Int. J. Hydrog. Energy 44 (2019) 4263. [18] Q. Cheng, D. Sun, X. Yu, J. Alloys Compd. 769 (2018) 167. [19] K. Liu, W. Zhou, D. Zhu, J. He, J. Li, Z. Tang, L. Huang, B. He, Y. Chen, J. Alloys Compd. 768 (2018) 269. [20] L. Kong, X. Li, X. Liao, K. Young, Int. J. Hydrog. Energy 44 (2019) 29338. [21] L. Ouyang, J. Huang, H. Wang, J. Liu, M. Zhu, Mater. Chem. Phys. 200 (2017) 164. [22] X. Zhao, L. Ma, Int. J. Hydrog. Energy 34 (2009) 4788. [23] L. Tong, J. Xiao, P. Bénard, R. Chahine, Int. J. Hydrog. Energy 44 (2019) 21055.

12