Analytical transfer function to simulate the dynamic response of the finite-length Warburg impedance in the time-domain

Based on fundamental electrode theory, an analytical transfer function to simulate the frequency impedance spectrum of the finite-length Warburg (FLW) impedance and the dynamic potential response of the FLW impedance in the time-domain has been developed in this study. Parameters reported in the literature estimated from experimental measurements carried out in polymer electrolyte fuel cells (PEFCs) have been considered to validate the new analytical transfer function. The analytical transfer function representing the FLW impedance can be considered in different equivalent electrical circuit configurations to simulate a more accurate dynamic output voltage of an electrochemical power system under the effect of diffusion phenomena. A Simulink model based on the Randles circuit and the new transfer function representing the FLW impedance is constructed to simulate the dynamic output voltage of a PEFC during a current-interrupt incident. In addition, a Simulink model based on an electrical circuit configuration and the new transfer function representing the FLW impedance is constructed to simulate the dynamic output voltage of a Li-ion battery. This study establishes a wider scope to relate the electrochemical impedance spectroscopy to the dynamic output voltage response of electrochemical power systems.


Introduction
Electrochemical impedance spectroscopy (EIS) is an experimental technique that can be applied in-situ to characterise physical processes of electrochemical power systems in the frequency-domain. The finitelength Warburg (FLW) impedance has been applied with EIS measurements to characterise diffusion processes in electrodes [1]. The mathematical equation of the FLW impedance was derived by Emil Warburg [2] and relates diffusion of electroactive species through an electrode surface. The conventional equation of the FLW impedance Z FLW = considers a resistance R W and time constant τ W for the diffusion of electroactive species in the electrode. Fig. 1 shows the impedance spectrum generated by the equation of the FLW impedance and presented in a complex-impedance (Nyquist) plot. In the high frequency range, a 45-degree straight line is apparent as shown in Fig. 1. The period of the sinusoidal signal during EIS measurements at high frequency is shorter than the time required for an electroactive species to diffuse through the finite diffusion layer. Hence, an infinite diffusion process is related to the 45-degree straight line in the impedance spectrum of the FLW. The time constant τ s where the imaginary component of the impedance spectrum of the FLW impedance reaches its minimum value at − Z j (the summit of the impedance spectrum) in the Nyquist plot is a factor 2.53 smaller than the characteristic time constant τ W represented in the FLW impedance equation τ W = 2.53τ s [3].
Simulink (MathWorks®) is a MATLAB-based graphical programming environment and a powerful programming environment for modelling and simulating multidomain dynamical systems such as electrochemical power systems. Simulink allows the representation of differential equations through the use of graphical block diagramming tools and a customizable set of block libraries. A transfer function can be defined with poles and zeros in the Laplace-domain and can be modelled in Simulink environment [4]. It is possible to simulate the input response of a transfer function in the time-domain using Simulink. It is not straight forward to represent the conventional equation of the FLW impedance with poles and zeros in the Laplace-domain. Although some approximations have been reported in the literature to represent the FLW impedance as a transfer function with poles and zeros in the Laplacedomain [5]. Boukamp [3] and Montella [6] reported that the FLW impedance can be well represented through the Voigt circuit considering an infinite combination of parallel resistor-capacitors. Eddine et al. [7] reported that it is possible to approximate the FLW impedance in the Laplace-domain through fractional order modelling. A fractional model to represent the FLW impedance has also been reported by Iftikhar et al. [8]. Gabano et al. [9] reported an impedance model based on fractional order modelling to simulate the bounded diffusion impedance represented in the impedance response of lithium batteries. Rubio et al. [10] reported an electrical circuit configuration based on resistors and capacitors connected in a parallel configuration to approximate the impedance of the FLW in the Laplace-domain. The aforementioned approaches reported in the literature can approximate the impedance response of the FLW with the response of a transfer function comprising poles and zeros. These approaches only provide an approximate representation of the impedance spectrum of the FLW impedance in the frequency-domain and may not be accurate enough to represent the dynamic potential response of the FLW impedance in the time-domain.
The Voigt circuit with an infinite series combination of parallel resistor-capacitor circuits can provide an exact representation of the impedance spectrum of the FLW in the frequency-domain [3,6] and could simulate the dynamic potential response of the FLW impedance in the time-domain. Nevertheless, the construction of the Voigt circuit with an infinite series combination of parallel resistor-capacitor circuits in Simulink environment may increase the computational effort and make the simulation process less effective for different modelling applications. For instance, during the fitting of a time-domain electrochemical model of a polymer electrolyte fuel cell (PEFC) with the experimental dynamic output voltage of the PEFC measured during a current-interrupt test [10]. Another disadvantage in using the aforementioned approaches reported in the literature is the fact that the diffusion time constant τ W defined in the conventional mathematical equation of the FLW impedance is commonly associated with the time constant from the electrical configuration comprising resistors and capacitors [7]. A transfer function representing the FLW impedance should be able to reproduce the characteristics (e.g. 45-degree straight line at high frequencies) of the impedance spectrum of the FLW impedance in the frequency-domain. Failing to do so, the transfer function will not simulate an exact representation of the dynamic potential response of the FLW impedance in the time-domain and will not be able to accurately simulate the effect of diffusion mechanisms on the dynamic output voltage of an electrochemical system under the effect of mass transport limitations.
The aim of this study is to develop a new transfer function in the Laplace-domain which can simulate the impedance spectrum of the FLW impedance in the frequency-domain and the dynamic potential response of the FLW impedance in the time-domain. The new transfer function is derived from fundamental electrode theory and considers the same parameters of the conventional mathematical equation of the FLW impedance such as diffusion resistance R W and diffusion time constant The aforementioned studies [7][8][9] do not provide guidance on how to apply the new analytical expression of the FLW impedance for the simulation of dynamic output voltage of electrochemical power systems and only a description of a generic algorithm containing extensive calculations to apply the analytical expression is presented. This may complicate the reproduction of the reported results and complicate the application of the new analytical expression in a programmable language carried out by early researchers or engineers from industry. In this study, a detailed guidance on how to apply the developed modelling architecture in a MATLAB/Simulink environment for the simulation of the dynamic voltage response of fuel cells and batteries considering parameters estimated from EIS measurements is presented. First the developed analytical transfer function is implemented in a MATLAB/ Simulink environment to simulate and study the dynamic potential response of the FLW impedance in the time-domain. Thereafter, the analytical transfer function representing the FLW impedance is considered in different electrical circuit configurations constructed in a Simulink environment to simulate a more accurate dynamic output voltage of the PEFC and battery under the effect of diffusion phenomena. A Simulink model based on the new transfer function representing the FLW impedance and the Randles circuit is constructed to simulate the dynamic output voltage of a PEFC during a change in current-step. Parameters from the Randles circuit reported in the literature [10,11] and estimated from experimental measurements carried out in PEFCs are considered in the Simulink model for the simulation of the dynamic output voltage of the PEFC. In the last section of the manuscript, it is demonstrated that the new transfer function representing the FLW impedance can also be considered in an electrical circuit configuration using a Simulink environment to simulate the dynamic output voltage of a Li-ion battery. The modelling architecture constructed in a Simulink environment requires parameters such as diffusion resistance and diffusion time constant from FLW impedance, charge transfer resistance, double-layer capacitance and ohmic resistance that can be estimated from EIS measurements carried out in the PEFC or in the battery at different operating conditions. The Simulink models presented in this study can also simulate the impedance response of the PEFC and battery; therefore, this study demonstrates that modelling can play an important role to study the relation between EIS and the dynamic output voltage of electrochemical power systems operated at different conditions.

Electrochemical phenomena in the cathode electrode
The transfer function representing the FLW impedance is derived from fundamental electrode theory. The following assumptions will be considered in the mathematical treatment: 1. The electrochemical reaction occurring at the cathode electrode is represented as a simple reaction or a single-step electron transfer process. 2. The anodic current represented in the current-overpotential equation is neglected because it can have a minimum contribution at large negative overpotential (Tafel approximation) [12]. 3. A decoupling between reactant transport limitation and charge transfer process is considered because the EIS technique can decouple physical processes with different time constants and different frequency dependence in electrochemical power systems.
The current-overpotential equation neglecting the anodic current contribution and neglecting reactant transport limitations can be expressed as [12]: where b is the Tafel slope with units V or expressed with units V/dec if b with units V is multiplied by a factor of 2.3, i 0 is the exchange current density, η is the overpotential defined as η = E − E eq , E is the potential of the electrode, and E eq is the potential at equilibrium conditions. Eq. (1) is also known as the Tafel equation. The use of a low amplitude value in the sinusoidal signal during EIS measurements allows the application of a linear model in order to analyse the electrochemical impedance spectra. Therefore, a linear model can be derived using the Taylor series expansion around the steady-state value [13] of Eq. (1) as such: where i F ∼ is the oscillating faradaic current density, E ∼ is the oscillating potential, R C is the charge transfer resistance during the reduction reaction at the cathode, and η SS is the overpotential at steady-state conditions. The denominator of Eq. (3) (charge transfer resistance) represents the steady-state current density [14]. It is possible to include the sinusoidal ratio φ ∼ between the oscillating concentration of electroactive species in the electrode surface c δ ∼ and oscillating bulk concentration of electroactive species c b ∼ in the oscillating faradaic current density expressed in Eq. (2) as such: Under low current density operation, the ratio of concentration of electroactive species is approximately φ ∼ ≈ 1. At high current density operation, the ratio of concentration of electroactive species is φ ∼ < 1. A value of φ ∼ less than 1 is attributed to the fact that the rate of consumption of electroactive species in the electrode surface to maintain the desired current density is higher than the rate of electroactive species being concentrated in the bulk solution c δ ∼ < c b ∼ . The inclusion of the parameter φ ∼ in Eq. (2) as expressed in Eq. (4) was supported by the fact that no coupling between charge transfer resistance and reactant transport resistance can exist during EIS measurements in electrochemical systems as the EIS technique allows the decoupling of physical processes with different time constants and different frequency dependence. Eq. (4) will allow the derivation of the conventional equation of the FLW impedance which has been broadly used to study mass transport limitations in PEFCs [11,15] and modern batteries [16].

Diffusion mechanisms
The electrochemical reaction takes place in the interface between dissimilar materials (e.g. metal/solution) in the electrode. The study of the reactant transport in the electrode can be treated with certain simplifications and approximations. A linear drop in the concentration of electroactive species across the diffusion distance of the film or solution surrounding the electrode allows the study of the mass transport through Fick's First law. Fick's Second Law [12] is considered to model the diffusion of electroactive species from the bulk solution to the electrode surface with respect to time, as such: ∂c(x, t) ∂t (6) where c is the local concentration of electroactive species diffusing through the film surrounding the electrode, D is the effective diffusion coefficient and x is the diffusion distance for electroactive species to diffuse from the bulk solution to the electrode surface. At steady-state conditions, it is possible to represent Eq. (6) in the Laplace-domain considering the bulk concentration in the solution c b at initial condition t = 0 as such: where the bar located over a letter represents the Laplace transform s of a variable. The mathematical solution of Eq. (7) through the method of undetermined coefficients for a nonhomogeneous linear differential equation [12] results in: with c b = c b /s and λ 1,2 = ± ̅̅̅̅̅̅̅̅ s/D √ Evaluating boundary conditions in Eq. (8) at the electrode surface c(δ, s) =c δ and at the bulk solution c(0, s) =c b yields: Substituting Eq. (9) and Eq. (10) into Eq. (8) yields: It is possible to define the relation between faradaic current density, charge transferred and the consumption of electroactive species as such: where v is the reactant flux, z is the number of electrons consumed during the electrochemical reaction and F is the Faraday constant. From Fick's First Law it is possible to establish that the flux of reactant is proportional to concentration gradient.
During steady-state conditions, the current density at which the electroactive specie is consumed in the electrochemical reaction from Faraday's Law is equal to the diffusion flux from Fick's First Law. Substituting the Laplace form of Eq. (4), and Eq. (12) into Eq. (13) yields: Differentiating Eq. (11) with respect to the finite diffusion distance x and considering hyperbolic trigonometric identities for exponentials yields: Substituting Eq. (15) into Eq. (14) yields: Rearranging Eq. (16) yields: The use of a small amplitude sinusoidal perturbation through the EIS technique allows the use of a linear model to represent the impedance of the electrode. It is possible to consider a linearized relation of the potential by considering the Laplace form Ē = E/s and considering E = RT/ zF [12]. In addition, considering the Laplace form c b = c b /s, considering and the transformation between the Laplace space and Fourier space s = jω in the term on the right-hand side of Eq. (17) yields: where Eq. (19) is the conventional equation of the FLW impedance [17] and represents the impedance response for reactant transport limitation in the frequency-domain, j is the imaginary component of a complex number and ω is the angular frequency with defined as the resistance for the diffusion process of electroactive species [15] and defined as the time constant to diffuse electroactive species towards the electrode surface [18].

Electrochemical impedance in the electrode
Substituting Eq. (18) into Eq. (4) yields the oscillating faradaic current density in terms of charge transfer resistance and finite-length Warburg impedance as such: A double-layer structure for the interface between the electrode and the solution is presented in electrochemical systems. At this interface electrons will be collected at the surface of the electrode and ions will be attracted to the solution. A potential difference is present because of the charge distribution between electrons and ions at the electrode-solution interface. This potential difference has a determinant role in the charge distribution within the reactants as well as in the position and orientation of the reactants to form the desired product. This double-layer can behave like a capacitor C dl that is in parallel with the electrode reactions. The current passing from the electrode to the solution either can take part in the charge transfer reaction or can contribute to the charge in the capacitive effect. The total current density is defined as the sum between faradaic current density and non-faradaic current density [19] (charge capacitance current density) i The electrochemical impedance of an electrochemical system neglecting the anode contribution can be defined as the ratio between oscillating potential and oscillating current density as such: where R e represents the ohmic resistance in the electrolyte that separates the anode and cathode, and the second term on the right-hand side of Eq. (24) represents the impedance of the cathode (Eq. (23)). Eq. (24) is analogous to the impedance of the Randles circuit considering the FLW impedance [13,15] as shown in Fig. 2. The Randles circuit could also consider a constant phase element (CPE) instead of the capacitor [20]. A CPE improves the quality of the fit between EIS measurements and the impedance response of the Randles circuit. Orazem and Triboller [13] discussed the physical meaning of the CPE when applying it with EIS measurements to study the impedance of electrochemical systems with time-constant dispersion. The Randles circuit shown in Fig. 2 has been considered to characterise the impedance response of PEFCs [21] and batteries [22]. More complex electrical circuit configurations based on transmission-line models have also been reported [23,24] to characterise the frequency response of porous electrodes in electrochemical systems.

Analytical transfer function for Warburg impedance
The solution of the inverse Laplace transform of Eq. (17) by the Heaviside's Expansion Theorem [25] is detailed in Appendix A and results in: The ratio between the concentration of electroactive species at the electrode surface c δ and bulk concentration in the solution c b can be expressed as a function of a mass transport resistance in the time-domain and charge transfer resistance by considering a linearized relation of potential [12] E = RT/zF and rearranging Eq. (25) as such: where R C is the charge transfer resistance represented in Eq. (3) and R M is the mass transport resistance in the time-domain as such: where R W is the resistance for the diffusion of electroactive species expressed in Eq. (20) and τ W is the diffusion time constant as expressed in Eq. (21). At this point, it would be possible to consider that Eq. (27) represents the FLW impedance in the time-domain. However, this can only be corroborated by comparing the frequency-domain response of Eq. (27) with the impedance response generated by the conventional mathematical equation of the FLW impedance (Eq. (19)). The potential related to the mass transport resistance can be expressed by multiplying the faradaic current density i F with Eq. (27) as such: The Laplace transform s of Eq. (28) results in The transfer function of Eq. (30) considering ī F = i F /s can be expressed as:

Fig. 2.
Randles circuit to characterise the impedance response of an electrochemical system neglecting anode contribution.
The frequency response of Eq. (31) represented in the Nyquist plot shown in Fig. 3 is simulated in MATLAB software considering the command nyquist(sys), where sys is a transfer function defined with poles and zeros in the Laplace-domain. The value of the parameters R W = 6.35 mΩ and τ W = 0.779 s are taken from the study reported by Rubio et al. [10] and were estimated by fitting an electrical circuit with the dynamic voltage of a PEFC during a current-interrupt test. The impedance spectrum generated by Eq. (31) is shifted towards the negative real component of the complex-impedance plot as shown in Fig. 3. In Appendix A, an indeterminate solution resulted in the Inverse Laplace transform of Eq. (17) by considering the pole (s 1 ) 3/2 = 0 in the Heaviside's Expansion Theorem. As a first approximation, the pole (s 1 ) 3/2 = 0 that yielded the indeterminate solution in the inverse Laplace transform was not considered in the mathematical treatment (see Appendix A). Therefore, this consideration in the mathematical treatment may have shifted the impedance spectrum shown in Fig. 3 towards the negative real component of the complex-impedance plot.
The correct transfer function to simulate the frequency response of the FLW impedance can be obtained by summing the resistance for the diffusion process expressed in Eq. (20) with Eq. (31) as such: The frequency response of Eq. (32) is simulated in MATLAB software and considering the script shown in Appendix B. The summation expression in Eq. (32) considered a value of n = 30 because minimum (negligible) changes on the simulated impedance response resulted at n > 30. The effect of the parameter n expressed in the summation expression of Eq. (32) on the impedance spectrum is shown in Appendix B. The parameter R W represented in the first term on the righthand side of Eq. (32) shifts the high-frequency end of the impedance spectrum shown in Fig. 3 towards the origin of the real axis. The impedance spectra shown in Fig. 4 labelled "Eq. (19)" and "Eq. (32)" were generated by Eqs. (19) and (32) considering the parameters R W = 6.35 mΩ and τ W = 0.779 s with a range of frequencies from 10 MHz to 0.001 Hz. The low frequency limit of the impedance spectra converged with the value of R W = 6.35 mΩ. The simulation of the impedance spectrum generated by Eq. (19) was also carried out in MATLAB software. Fig. 4 shows that Eq. (32) can reproduce the same characteristics of the impedance spectrum generated by the conventional equation (Eq. (19)) of the FLW impedance. The transfer function represented in Eq.
(32) can simulate the 45-degree straight-line feature at high frequencies attributed to the infinite diffusion process.

Current input response of the Warburg impedance in the time-domain
Eq. (32) can simulate the impedance spectrum of the FLW impedance and can simulate the current-input response of the FLW impedance in a MATLAB/Simulink environment. The simulation process in MATLAB/ Simulink to simulate the current-input response of the FLW impedance in the time-domain is carried out in two steps:    Further information about the "LTI System block" can be found elsewhere [26].
A step-change in current with values I = 0 A at t = 0 s and I = 13 A at t = 1 s has been considered as inputs in the Simulink model shown in where the imaginary component of the impedance spectrum of the FLW impedance reaches its minimum value − Z j (the summit of the impedance spectrum) in the complex-impedance plot can be represented in the dynamic potential shown in Fig. 5b. The time constant τ s is a factor 2.53 smaller than the characteristic time constant τ W represented in the FLW impedance equation [3,27] as such: τ W = 2.53τ s . The dynamic potential response had reached 68.5 % of the steady-state value at the time constant τ s = 1 + 0.3079 s and at τ W = 1 + 0.779 s the potential dynamic response had reached 92.1 % of its final value as shown in Fig. 5b. If the FLW impedance is approximated by the impedance of an electrical circuit configuration comprising resistors and capacitors, the diffusion time constant would be approximated to the time constant from resistors and capacitors τ W = R W C W [7]. Under this consideration, the potential response of the FLW impedance would be equivalent to the response of a first order system and will reach 63.2 % of the steady-state value [4] at τ W = R W C W (diffusion time constant related to resistor R W and capacitor C W ). Approximating the FLW impedance with the impedance of a circuit considering resistors and capacitors may incur an erroneous interpretation of the diffusion time constant on the dynamic potential and incorrect calculation/simulation of the dynamic potential of the FLW impedance in the time-domain.
The diffusion coefficient D plays an important role in the diffusion process of electroactive species to reach the active zones of the electrode during the redox reaction. The diffusion coefficient D is inversely proportional to the diffusion resistance R W (Eq. (20)) and diffusion time constant τ W (Eq. (21)). Increasing the bulk concentration c b of electroactive species reduces the diffusion resistance R W (Eq. (20)) and the FLW impedance. Fig. 6 shows the effect of decreasing the diffusion parameters D and c b on the dynamic potential estimated by the analytical transfer function Eq. (32) simulated by the Simulink model shown in Fig. 5a. The diffusion coefficient D was reduced by 50 % which is equivalent to an increase in diffusion resistance R W and diffusion time constant τ W by a factor of 2 in Eq. (32). Fig. 6 shows that decreasing the diffusion coefficient D increases the steady-state potential and the time for the curve to reach the steady-state value. The bulk concentration c b was reduced by 50 % and this reduction is equivalent to an increase in diffusion resistance R W by a factor of 2 in Eq. (32) where Fig. 6 also shows that decreasing the bulk concentration c b increases the steadystate potential.
The transfer function represented in Eq. (32) provides an exact representation of the impedance spectrum of the FLW impedance in the frequency-domain as shown in Fig. 4 therefore, it can confidently simulate the input-current response of the FLW impedance in the timedomain as demonstrated in Fig. 5b. The transfer function represented in Eq. (32) can be considered in the configuration of the Randles circuit and thus can simulate the dynamic output voltage of a PEFC during a current-interrupt incident. This will be demonstrated in the next section.

Current-input response of the Randles circuit in the timedomain
The impedance of the Randles circuit in the frequency-domain is represented in Eq. (24). Correa et al. [28] simulated the dynamic output voltage of a PEFC considering a steady-state model (polarisation curve model) of a PEFC with the Randles circuit shown in Fig. 2 but neglecting the FLW impedance. The authors attributed the dynamic output voltage of the PEFC to the double-layer capacitance in the electrodes and the effect of the diffusion time constant on the dynamic output voltage of the PEFC was neglected in the study. Fig. 7 shows an equivalent electrical circuit to simulate the voltage output of a PEFC during a step-change in current. The electrical circuit shown in Fig. 7 neglects the contribution of the anode on the voltage output of the PEFC because the hydrogen oxidation reaction is a less complicated and more facile reaction sequence than the oxygen reduction reaction (ORR) [29,30]. However, under certain operating conditions such as start-up events and long-term operation, anode losses can be increased [31,32] and should not be neglected in the electrical circuit configuration representing the PEFC. The electrical circuit shown in Fig. 7 can be considered to simulate the dynamic output voltage of a PEFC during a step-change in current [33] or during a current-interrupt test [34]. R e represents the ohmic resistance of the PEFC and considers the electron and proton opposition across the different layers of the PEFC, R C represents the charge transfer resistance during the ORR, C dl represents the double-layer capacitance in the cathode catalyst layer (CCL), and Z W represents the opposition of reactant transport from the gas channel to the electrode surface which can be represented by Eq. (32).
The parameters of the electrical circuit shown in Fig. 7 can be estimated from EIS measurements carried out in a PEFC.
If Kirchhoff's laws are applied to the electrical circuit shown in Fig. 7, the voltage output of the PEFC V cell can be calculated through the following equations: where E r is the reversible potential of the PEFC and can be represented through the open-circuit voltage equation [35], I cell is the current demand of the PEFC, E F is the overpotential (losses) during the ORR in the cathode and represents the potential in the parallel configuration of the electrical circuit considering the charge transfer resistance R C , transfer  function representing the FLW impedance Z W (Eq. (32)), and doublelayer capacitance C dl . One of the advantages of the Simulink environment is the fact that time-domain differential equations representing the transient response of a system can be defined in the Laplace-domain to simulate the inputresponse of the system in the time-domain. Fig. 8 shows the Simulink architecture to represent Eqs. (33) and (34) in the Laplace-domain.

Voltage output of a PEFC during a current-interrupt incident
The Simulink model shown in Fig. 8 can simulate the dynamic output voltage of the PEFC considering a current-step as an input. Prior to the execution of the Simulink model shown in Fig. 8, the values of the parameters E R , R C , C, R e and Z W have to be defined and loaded in the workspace of MATLAB. Previously in Section 3.1 it was discussed that the solution of the transfer function (Eq. (32)) representing the FLW impedance is loaded in the workspace of MATLAB when the script shown in Appendix B is executed. The simulation of the current-input response of Eq. (32) in Simulink environment has been demonstrated in Fig. 5.
Before carrying out the simulation of the dynamic output voltage of the PEFC using the Simulink model shown in Fig. 8, it is compulsory to demonstrate that the Simulink model is able to reproduce the same impedance spectrum in the frequency-domain generated by Eq. (24) (Impedance of the Randles circuit). The Simulink model considers the transfer function representing the FLW impedance defined in Eq. (32). Eq. (24) considers the conventional equation of the FLW impedance defined in Eq. (19). The considered parameters for the simulation of the impedance spectra generated by Eq. (24) and generated by the Simulink model (Fig. 8) are: E R = 0.92 V, R e = 0.0139 Ω, R C = 0.0053 Ω, C = 0.989 F, R W = 0.00635 Ω, and τ W = 0.779 s. These parameters were reported by Rubio et al. [10] and were estimated by fitting an electrical circuit with the dynamic output voltage of a PEFC measured during a current interrupt test.
Zhivomirov [36] reported a MATLAB script to calculate the impedance frequency response of a two-terminal circuit using the DAQ-system NI USB-6211. The MATLAB script reported by Zhivomirov [36] was considered to calculate the impedance frequency response of the Simulink model shown in Fig. 8. The MATLAB script defines a sinusoidal signal input at different frequencies and calculates the Fourier transform of the time-domain input and output signals. The simulation of the impedance frequency response of the Simulink model shown in Fig. 8 and using the MATLAB script reported by Zhivomirov [36] is carried out through the following steps: 1. The MATLAB script sends a sinusoidal current input to the Simulink model (Fig. 8) at different frequencies and the resulting sinusoidal voltage output from the Simulink model is read by the MATLAB script. 2. The MATLAB script calculates the impedance response of the Simulink model in the frequency-domain considering the Fourier transform of the sinusoidal input and output signals at different frequencies. Fig. 9 shows a comparison between the impedance spectrum generated by Eq. (24) and the impedance spectrum generated by the Simulink model shown in Fig. 8 considering the MATLAB script from Zhivomirov [36]. The considered range of frequencies was from 1 kHz down to 0.02 Hz. The impedance spectrum labelled "Eq. (24)" can be simulated by defining Eq. (24) in MATLAB or by constructing the Randles circuit shown in Fig. 2 in ZView software (Scribner Associates Inc.). Two overlapped semicircles are represented in the impedance spectra shown in Fig. 9. The diameter of the first semicircle in the high-low frequency region is attributed to the value of the charge transfer resistance R C during the ORR. The diameter of the low frequency semicircle is attributed to the value of the diffusion resistance R W represented in the FLW impedance. In some cases, the separation of the semicircles cannot be visualised on the overall impedance spectrum and this effect has been attributed to a low value of the diffusion time constant τ W [37].
The low frequency resistance can be calculated with the sum of the ohmic resistance, charge transfer resistance, and diffusion resistance as such R e + R C + R W = 0.0255 Ω. The high frequency resistance represents the ohmic resistance R e = 0.0139 Ω.
The comparison between an experimental impedance spectrum and the simulated impedance spectrum generated by an impedance model can be assessed through the real and imaginary relative residuals of the simulated and measured impedance data. The comparison of the impedance spectra shown in Fig. 9 can be assessed through Eqs. (35) and (36) which represent the average of the real and imaginary relative residuals of the simulated data generated by the Randles circuit (Eq. (24)) and the Simulink model (Fig. 8). Eqs. (35) and (36) have been adapted and considered from the study of Song and Bazant [38]. . 8. Simulink architecture to simulate the voltage output of the PEFC considering the current as input. Er is the reversible potential, Re is the ohmic resistance, Rc is the charge transfer resistance, Zw is the transfer function representing the FLW impedance (Eq. (32)), and C is the double-layer capacitance.

Fig. 9.
Comparison between impedance response generated by the Randles circuit represented in the frequency-domain (Eq. (24)), and impedance response generated by the Simulink model (Fig. 8) using the MATLAB script reported by Zhivomirov [36].
where K is the number of the considered frequencies ω, Z randles is the simulated data generated by the Randles circuit (Eq. (24)), Z Rrandles and Z Irandles are the real and imaginary components of the simulated data generated by the Randles circuit (Eq. (24)) and correspond to the frequency ω k , where Z Rsimulink and Z Isimulink are the real and the imaginary components of the simulated data generated by the Simulink model (Fig. 8). The comparison between the impedance spectra generated by Eq. (24) and the Simulink model shown in Fig. 9 is evaluated by the average of the real (Eq. (35)) and imaginary (Eq. (36)) relative residuals. The sum of Eqs. (35) and (36) ΔZ = ΔZ R + ΔZ I resulted in 7.1688 × 10 − 4 . A good agreement between the impedance spectra is obtained when the real and imaginary relative residuals have a minimum value [38]. The dynamic voltage shown in Fig. 10 was simulated by the Simulink model shown in Fig. 8 and considering the following parameters: E R = 0.92 V, R e = 0.0139 Ω, R C = 0.0053 Ω, C = 0.989 F, R W = 0.00635 Ω, and τ W = 0.779. These parameters were reported by Rubio et al. [10] and were estimated by fitting an electrical circuit with the dynamic voltage measured during a current-interrupt test in a PEFC. The electrical circuit reported by Rubio et al. [10] considered an approximation of the FLW impedance through a parallel configuration between resistors and capacitors. Thereafter, the authors simulated the impedance spectrum generated by the proposed electrical circuit. A step-change in current was considered in the Simulink model shown in Fig. 8. A current value of I cell = 13 A was initially considered at t = 0 s and the voltage output of the PEFC at steady-state resulted in 0.587 V as shown in Fig. 10. Thereafter the current was stepped down I cell = 0 A at t = 5 s.
Stepping the current to a zero value can be analogous to a currentinterrupt test during PEFC operation. Fig. 10 shows the dynamic output voltage of the PEFC during a current-interrupt incident.
The simulation of the voltage output of the PEFC shows a verticalline feature as the current is stepped down to zero. The magnitude of the vertical-line feature of the dynamic voltage at t = 5 s can be attributed to ohmic losses in the membrane of the PEFC [39]. The potential attributed to the ohmic resistance resulted in I cell R e = 0.1807 V. The magnitude of the vertical-line feature of the PEFC voltage V cell is approximately 0.1807 V as shown in Fig. 10. The magnitude of the vertical-line feature of the PEFC voltage V cell is proportional to the ohmic resistance R e . The vertical-line feature is related to electrochemical mechanisms in the membrane (dielectric relaxation time [40]) of the PEFC with a small time constant. Buchi et al. [40] reported that the available time window for membrane resistance measurements using current-interrupt test is between 5 × 10 − 10 s and 10 − 8 s. Measurement equipment with high bandwidth is required to record the fast transient response of the PEFC during a current-interrupt test [41].
An exponential increase in voltage V cell until reaching the value of the reversible potential E R = 0.92 V is observed in Fig. 10 32)) is fitted to the experimental dynamic voltage measured from a current-interrupt test in a PEFC. This however is the aim of future work. The estimation of the parameters of the Simulink model from experimental data can be achieved with the use of the Simulink toolbox "Simulink Design Optimization" [42].

Voltage output of a PEFC with and without microporous layers
Malevich et al. [11] reported EIS measurements carried out in a 100 cm 2 H 2 /air PEFC with and without microporous layers (MPLs) attached to the gas diffusion layers (GDLs) and operated at a current density of 0.5 A/cm 2 . The authors fitted the Randles circuit considering FLW impedance with the EIS measurements. The parameters reported by Malevich et al. [11] are shown in Table 1 and are considered in the Simulink model represented in Fig. 8 to simulate the effect of the MPLs on the dynamic output voltage of the PEFC. The parameters shown in Table 1 demonstrate that it is not straight forward to determine the role of the MPLs on PEFC performance because the diffusion parameters estimated from EIS-PEFC measurements can be related to either the GDL or CCL. The MPL could have removed water from the GDL and increased the amount of liquid water in the CCL. Under this condition a reduction of the oxygen diffusion resistance R W in the GDL is expected as shown in Table 1. Nevertheless, a lower oxygen diffusion time constant τ W in the PEFC with MPLs compared with the PEFC without MPLs would be expected because the diffusivity of oxygen through the GDL is increased and the effective diffusion distance is reduced. The relation between oxygen diffusion resistance R W and oxygen diffusion time constant τ W for PEFCs with and without MPLs shown in Table 1 has been discussed in the study of Malevich et al. [11]. A change in double-layer capacitance in PEFCs considering MPLs compared with PEFCs without MPLs would be expected and Roy and Orazem [43] reported an increase of the Fig. 10. Simulated dynamic voltage generated by the Simulink model (Fig. 8). The dynamic voltage represents the voltage output of a PEFC during a currentinterrupt test. The current I cell is stepped down from 13 A to 0 A. interfacial capacitance in PEFCs with MPLs. Malevich et al. [11] reported the parameters shown in Table 1, but did not report the value of the double-layer capacitance C dl in the CCL; therefore as a first approximation, the capacitance C dl = 0.372 F reported by Rubio et al. [44] will be considered in the simulation of the voltage output of the PEFC without MPLs. The capacitance C dl = 0.372 F reported by Rubio et al. [44] was estimated by fitting the Randles circuit with EIS measurements carried out in a 100 cm 2 H 2 /air PEFC. An increase of capacitance 3C dl will be considered in the simulation of the voltage output of the PEFC with MPLs. The increase of capacitance by 3C dl is consistent with the results reported by Roy and Orazem [43] for a PEFC with MPLs compared with a PEFC without MPLs operated at 0.5 A/cm 2 . Fig. 11 shows the impedance spectrum generated by the Randles circuit shown in Fig. 2 represented by Eq. (24) and the impedance spectrum generated by the Simulink model shown in Fig. 8. The parameters shown in Table 1 and the range of frequencies from 2 kHz down to 0.02 Hz were considered in the simulated impedance spectra shown in Fig. 11. The MATLAB script reported by Zhivomirov [36] was considered to simulate the impedance spectra generated by the Simulink model shown in Fig. 8. As previously discussed, the MATLAB script reported by Zhivomirov [36] defines a sinusoidal signal input at different frequencies and calculates the Fourier transform of the time-domain input and output signals of the Simulink model shown in Fig. 8. The diameter of the first semicircle from the high-low frequency region represents the charge transfer resistance R C during the ORR, and the diameter of the second semicircle in the low frequency range represents the oxygen diffusion resistance R W . The impedance spectrum related to the PEFC with MPLs presents a lower magnitude than the impedance spectrum related to the PEFC without MPLs. The semicircles representing the charge transfer process and oxygen diffusion process can be masked and overlapped in a single impedance spectrum at a very low value of the oxygen diffusion time constant [37]. The sum ΔZ = ΔZ R + ΔZ I of the average of the real (Eq. (35)) and imaginary (Eq. (36)) relative residuals of the simulated impedance spectra shown in Fig. 11 resulted in ΔZ = 0.0027 for the case PEFC with MPL and ΔZ = 0.0031 for the case PEFC without MPL. Fig. 12 shows the output voltage of the PEFC simulated by the Simulink model shown in Fig. 8 considering the parameters shown in Table 1. The simulation results shown in Fig. 12 were carried out with a step-change in current from I = 50 A at t = 0 s to I = 0 A at t = 5 s. The simulated output voltage during a step-change in current at I = 0 A is analogous to the measured voltage of the PEFC during a currentinterrupt test. A lower output voltage of the PEFC at I = 50 A results when the Simulink model considers the parameters from the PEFC without MPL. This low value of the output voltage of the PEFC without MPL compared with the voltage of the PEFC with MPL is attributed to the magnitude of the impedance spectra as shown in Fig. 11. The profile of the dynamic voltage of the PEFC with and without MPLs is observed in Fig. 12 where the steady-state voltage at 5.5 s represents the reversible potential of the PEFC at I = 0 A. The simulated output voltage of the PEFC in Fig. 12 is a function of the time constant related to the doublelayer capacitance τ C = [R C + R W ]C dl and the oxygen diffusion time constant τ W . Fig. 12 demonstrates the importance of considering a transfer function of the FLW impedance in the electrical circuit representing the PEFC to simulate the phenomenological processes with different time constants in the dynamic output voltage of the PEFC. The Simulink model shown in Fig. 8 could also simulate the effect of electrochemical mechanisms such as Tafel slope b, diffusion coefficient D, diffusion distance δ, bulk concentration c b , double-layer capacitance C dl on the dynamic output voltage of the PEFC.

Other applications
The Simulink model shown in Fig. 8 could incorporate the dependency of the electrochemical parameters of the PEFC with different operating conditions such as pressure, temperature and flow rate. For instance, the oxygen diffusion coefficient in the GDL can be defined as a function of the operating temperature and pressure of the PEFC through the empirical relation reported by Iranzo et al. [45]. The electrochemical model represented in a Simulink environment and shown in Fig. 8 can also be implemented with a Simulink model representing the balance of plant of a PEFC system [46]. It would be possible to simulate the dynamic output voltage of the PEFC during the dynamic change in reactant pressure, reactant flow rate and PEFC temperature. It could be a valuable simulation tool for the development and optimization of different control strategies to ensure optimal PEFC performance at different operating conditions [47]. The Randles circuit comprising the FLW impedance can also be considered to study the phenomenological processes represented in EIS measurements in solar cells [48,49]. The Simulink model shown in Fig. 8 could also simulate and study the dynamic transient response of solar cells attributed to the diffusion of electrical carriers across the diffusive layers of the solar cell.

Voltage output of a Li-ion battery
The FLW impedance has been connected in series with a resistor/ capacitor parallel configuration to represent the impedance of a Li-ion battery [7]. The FLW impedance has also been considered in a transmission line model [16] and Randles circuit [50] to simulate the impedance response of a Li-ion battery. Different electrical circuit configurations have also been reported to represent the impedance response of Li-ion batteries [51][52][53]. Surya et al. [54] reported a method to estimate the state of charge (SoC) in Li-ion batteries using MATLAB/ Simulink software and the Randles circuit shown in Fig. 2 but neglecting the FLW impedance. A better prediction of the SoC in batteries could Fig. 11. Impedance spectra representing the PEFC with and without MPLs. The impedance spectra are generated by Eq. (24) and Simulink model (Fig. 8) considering the parameters from Table 1. result if the equivalent electrical circuit representing the battery considers the diffusion impedance of electrical chargers. It is possible to construct a Simulink model considering the transfer function of the FLW (Eq. (32)) and other electrical components to simulate the impedance response of a Li-ion battery in the frequency-domain and the dynamic output voltage of the Li-ion battery in the time-domain. Fig. 13a shows the electrical circuit configuration to represent the voltage output of the Li-ion battery. A CPE is considered in the electrical circuit configuration to simulate the impedance response of the Li-ion battery represented at the lowest frequencies during EIS measurements. The impedance of a CPE is represented as: Y represents a parameter related to the CPE with units s P /Ω, and the superscript P represents a parameter to account for roughness at the blocking interface of the electrode [55]. For the specific case if superscript P is equal to 1 in Eq. (37); Eq. (37) would represent the impedance of an ideal capacitor with Y having units F. Eq. (37) can be represented in the Laplace-domain by considering s = jω. The CPE is connected in series with the rest of the electrical components in the electrical circuit configuration representing the impedance of the Li-ion battery as shown in Fig. 13a. The Simulink model of the electrical circuit representing the impedance of the Li-ion battery is shown in Fig. 13b. The Simulink model comprises the transfer function representing the FLW impedance defined in Eq. (32), a transfer function to represent the impedance of the parallel connection between the charge transfer resistance R a and capacitance C a in the anode, and a Simulink block representing a fractional operator to account for the frequency dispersion when the parameter P from the CPE is different to 1, (Eq. (37) considering s P with P ∕ = 1 and s = jω). The fractional operator block considered in the Simulink model shown in Fig. 13b was taken from the Simulink library FOTF Toolbox developed by Xue [56]. FOTF Toolbox has been designed for the modelling analysis and the design of fractional-order control systems and can be downloaded from MathWorks [57]. The fractional operator block shown in Fig. 13b allows the simulation of the impedance response of the CPE.
In a previous study [16], EIS measurements were carried out in a Li-Po battery at an open circuit voltage (OCV) of 11.44 V with a range of frequencies from 20 kHz to 0.002 Hz. Fig. 14 shows the experimental impedance spectrum of the Li-Po battery [16]. The electrical circuit configuration shown in Fig. 13a was constructed in ZView software (Scribner Associates Inc.) and was fitted to the experimental impedance spectrum of the Li-Po battery as shown in Fig. 14. ZView software considers the conventional equation of the FLW impedance represented in Eq. (19). The resulting parameters from the fitting process carried out in ZView software were R e = 0.0187 Ω, R a = 0.0031 Ω, C a = 3.19 F, R c = 0.0057 Ω, C c = 0.229 F, R W = 0.0156 Ω, τ W = 17.87 s, Y = 892.2 s P Ω − 1 , and P = 0.7. The aforementioned parameters were considered in the Simulink model shown in Fig. 13b to simulate the output voltage of the Li-Po battery in the time-domain. Before proceeding with the simulation of the voltage output of the Li-Po battery, the MATLAB script reported by Zhivomirov [36] was considered to simulate the impedance response of the Simulink model shown in Fig. 13b in the frequency-domain. The MATLAB script reported by Zhivomirov [36] was previously considered in Section 4.1. The MATLAB script defines a sinusoidal signal input at different frequencies and calculates the Fourier transform of the timedomain input and output signals of the Simulink model shown in Fig. 13b. Fig. 14 shows that the Simulink model shown in Fig. 13b can simulate the measured impedance spectrum of the Li-Po battery and also shows the contribution of the FLW impedance on the impedance spectrum. A disagreement between the experimental spectrum and the simulated impedance spectrum from the Simulink model is present at the lowest frequencies as shown in Fig. 14. This disagreement in the  Comparison between experimental impedance spectrum of the Li-Po battery, simulated impedance spectrum generated by the Simulink model (Fig. 13b) and fitted impedance spectrum considering ZView software.
lowest frequency region is mainly attributed to the fractional operator block shown in Fig. 13b representing the Laplace operator 1/s P of the CPE. The fractional operator block in Simulink was developed by Xue [56] and allows the user to select different filters such as Oustaloup and Matsuda-Fuki filters as well as the approximation order of the Laplace operator s P . Perhaps a correct selection of the filter and approximation order in the fractional operator block could improve the simulation results from the Simulink model shown in Fig. 14 at the lowest frequencies (impedance of CPE). The comparison between the experimental impedance spectrum and the impedance spectrum generated by the Simulink model shown in Fig. 14 is evaluated by the average of the real (Eq. (35)) and imaginary (Eq. (36)) relative residuals. The sum of Eqs. (35) and (36) ΔZ = ΔZ R + ΔZ I resulted in 0.025. The resulting value of ΔZ could be reduced by improving the simulated response of the fractional operator block (CPE) represented at the lowest frequencies of the impedance spectrum as shown in Fig. 14.
The aim of this study is not to discuss the values of the parameters from the electrical circuit representing the impedance of the Li-ion battery as different electrical circuit configurations can reproduce the same set of EIS measurements. For instance, in the previous study [16], the experimental impedance spectrum shown in Fig. 14 was analysed with an electrical circuit configuration consisting of a transmission-line model, FLW impedance and blocked-diffusion Warburg impedance. The aim of this study is to demonstrate that the developed transfer function representing the FLW impedance (Eq. (32)) can be considered in the electrical circuit configuration representing the Li-Po battery to simulate the impedance response in the frequency-domain and dynamic output voltage in the time-domain as demonstrated in Figs. 14 and 15. The dynamic output voltage of the Li-Po battery is shown in Fig. 15 and is simulated with the Simulink model shown in Fig. 13b. A constant current-input of 2 A is considered in the Simulink model where a decrease in the voltage output of the battery starting at the OCV of 11.44 V is shown in Fig. 15. The dynamic output voltage of the Li-Po battery does not reach a steady-state value. This behaviour is related to the fact that the impedance of the battery is increased with decreasing frequency during EIS tests. The time constant τ W from the FLW impedance related to the diffusion of lithium ions in the solution-phase of the cathode can be represented in the dynamic output voltage of the Li-Po battery. Following the methodology in MATLAB/Simulink reported by Surya et al. [54], the Simulink model shown in Fig. 13b considering the transfer function representing the FLW impedance (Eq. (32)) could be a valuable tool to provide a better prediction of the SoC in the Li-ion battery.

Conclusions
Based on electrode theory, an analytical transfer function to represent the FLW impedance has been derived. The analytical transfer function can simulate the impedance spectrum of the FLW impedance in the frequency-domain. The developed transfer function representing the FLW impedance can be implemented in a Simulink environment to simulate the current-input response of the FLW impedance in the timedomain. Parameters reported in the literature estimated from experimental measurements carried out in PEFCs have been considered to validate the new analytical transfer function. This study demonstrated that it is possible to simulate the dynamic output voltage of a PEFC during a step-change in current using a Simulink model based on the Randles circuit and considering the developed transfer function representing the FLW impedance. In addition, this study demonstrated that it is possible to simulate the dynamic output voltage of a Li-Po battery using a Simulink model based on an electrical circuit configuration considering the developed transfer function for the FLW impedance. The transfer function representing the FLW impedance can be considered in different electrical circuit configurations to simulate the dynamic output voltage of electrochemical power systems such as fuel cells, batteries, and solar cells. It is therefore possible to obtain an insight into the effect of the diffusion mechanisms on the potential dynamic response of electrochemical systems.

Declaration of competing interest
No conflict of interests exists.

Data availability
Data will be made available on request.

Appendix B. MATLAB script for simulation of the impedance spectrum of the transfer function representing the FLW impedance
The effect of the parameter n expressed in the summation expression of Eq. (32) on the impedance spectrum is shown in Fig. B1.