**3. Analytical fitting of electrochemical parameters**

The proposed analytical modeling of electrochemical reactions is founded on four building blocks. Figure 1 illustrates the structure of the modeling methodology.


#### **3.1 The mathematical expression for the proposed reaction and transport models**

Once the reaction model is defined, the transport model must be included. Using an RDE as the working electrode for the LSV experiments, the transport equations can be simplified. In addition, assuming that the electrolyte is a diluted solution, the migration term in the transport model can be neglected. This section provides the basic equations for mass and

Fig. 1. The four building blocks of the modeling methodology.

4 Electrochemical Cells

convective-diffusion equation for a rotating disk electrode have been solved rigorously for the steady state (Levich, 1962; Slichting, 1979). The axial symmetry of the configuration of the RDE reactor and the uniform current distribution allow a one-dimensional description. Moreover, at sufficient flow rate (when natural convection can be ignored), the hydrodynamics in diluted solution are not influenced by changes in concentrations due to electrochemical reactions. The mathematical problem can thus be solved more easily. Levich reduced the equation of convection transport to an ordinary differential equation (Albery & Hitchman, 1971; Levich,

To model an electrochemical reaction and determine its mass and charge transfer parameters quantitatively, an electrochemical data fitting tool has been developed in our research group. From an analytical approach, it is designed to extract a quantitative reaction mechanism from

The proposed analytical modeling of electrochemical reactions is founded on four building

**The results of the experimental study** These are the current-potential couples defining the polarization curve. The mean of multiple experiments that are performed under identical conditions is the experimental data for the modeling. The standard deviation of the

**The mathematical expression for the proposed model** The proposed model is based on well-considered reaction and transport models for the studied reaction. The mathematical expression of the reaction-transport model is derived from the basic equations that describe what happens during an electrochemical reaction. It has the following form: *current = function (potential, experimental parameters, model parameters)*, where the *experimental parameters* describe the experimental conditions, like e.g. temperature, rotation speed of the RDE, concentration, . . . , and the *model parameters* are the unknown parameters that need to be quantitatively determined, like e.g. rate constants, transfer coefficients, .... **The fitting procedure** In this block the differences between experimental and theoretical data are minimized. A weighted least squares cost function is formulated. The Gauss-Newton and Levenberg-Marquardt method are implemented to minimize this cost function and eventually provide the parameter values which best describe the data. Moreover, the

**A statistical evaluation** If a statistical evaluation of the fitting results demonstrate a good description of the experiment by the model, a quantitative reaction mechanism is obtained. If, on the other hand, no good agreement between experiment and model is achieved, a

Once the reaction model is defined, the transport model must be included. Using an RDE as the working electrode for the LSV experiments, the transport equations can be simplified. In addition, assuming that the electrolyte is a diluted solution, the migration term in the transport model can be neglected. This section provides the basic equations for mass and

new mechanism has to be proposed and the previous steps should be repeated.

**3.1 The mathematical expression for the proposed reaction and transport models**

1962; Slichting, 1979).

polarization curves.

**3. Analytical fitting of electrochemical parameters**

experiments is used in the fitting procedure.

blocks. Figure 1 illustrates the structure of the modeling methodology.

standard deviations of the estimated parameters are also calculated.

charge transfer, which can be found in numerous textbooks (Bamford & Compton, 1986; Diard et al., 1996; Newman, 1973; Pletcher, 1991; Thirsk & Harrison, 1972; Vetter, 1967).

Consider a uniformly accessible planar electrode, immersed in an electrolyte that contains electroactive species and an excess of inert supporting electrolyte. At the surface an electrochemical reaction is taking place, which has *P* partial heterogeneous electrochemical or chemical reactions with *Nv* electroactive species in the electrolyte or in the electrode material and *Ns* electroactive species present in an adsorbed phase on the electrode surface. *N* is the total number of electroactive species involved in the reaction: *N* = *Nv* + *Ns*. The *j th* step of the reaction can be written as:

$$\sum\_{i=1}^{K\_j} r\_{ij} \mathbf{X}\_i \leftrightarrow \sum\_{i=1}^{N} p\_{ij} \mathbf{X}\_i \pm n\_j \mathbf{e} \tag{1}$$
 
$$\mathbf{K}\_j^{\prime}$$

with *rij* and *pij* the stoichiometric coefficients, the index *i* refers to the considered species, the index *j* to the partial reaction. *Kj* and *K*� *<sup>j</sup>* are the rate constants. For an electrochemical reaction, they depend on the electrode potential. *nj* is the number of electrons exchanged in the *j th* partial reaction. For an electrochemical reaction *nj* is preceded by a plus sign if the reaction is written in the sense of the oxidation and by a minus sign if written in the sense of the reduction. For a chemical reaction *nj* equals zero.

The global reaction is described by the relations that connect the electrode potential *E* to the faradaic current density *if* , the interfacial concentrations of the volume species *Xi*(0) and the surface concentrations *Xi* of the adsorbed species1. Under steady-state conditions and using an RDE, the general system of equations describing this electrochemical systems is given by:

• the rate (*vj*) expressions for each partial reaction:

$$w\_j = K\_j \prod\_{i=1}^{N} X\_i^{r\_{ij}} - K\_j' \prod\_{i=1}^{N} X\_i^{p\_{ij}} \tag{2}$$

or

species. *X*∗

of the RDE.

**3.2 The fitting procedure**

the electrochemical reaction.

function.

*vXi*

function of the mass transport rate constant and the bulk concentration:

1986; Pintelon & Schoukens, 2001; Press et al., 1988; Sorenson, 1980).

*Ilim* = 0.620*nFSD*2/3

equation thus applies to the totally mass-transfer-limited condition at the RDE.

for diffusion and convection of species *Xi*, given by:

*i* = *Nv* + 1, . . . , *N*

*<sup>i</sup>* is the bulk concentration of species *Xi*. *mXi* is the mass transport rate constant

*Xi <sup>ν</sup>*−1/6*ω*1/2*X*<sup>∗</sup>

were *JXi* is the molecular flux (expressed in mol/m2 s) of *Xi*, equal to the number of moles of *Xi* going per unit of time across a unit plane, perpendicularly oriented to the flow of the

<sup>27</sup> Modeling and Quantification of Electrochemical Reactions in RDE (Rotating Disk Electrode) and IRDE (Inverted Rotating Disk Electrode) Based Reactors

with *DXi* the diffusion coefficient of *Xi*, *ν* the kinematic viscosity and *ω* the rotation speed

At the maximum diffusion rate, the concentration at the electrode is zero and the limiting current is achieved. The Levich equation provides an expression of the limiting current as a

with *S* the electrode surface and *n* the number of electrons exchanged in the reaction. This

Many comprehensive textbooks about parameter estimation and minimization algorithms are available. The development of the fitting procedure for the analytical modeling of electrochemical reactions is founded on a few of them (Fletcher, 1980; Kelley, 1999; Norton,

Given a set of observations, one often wants to condense and summarize the data by fitting it to a 'model' that depends on adjustable parameters. In this work the 'model' is the current-potential relation, describing the polarization curve, which is derived from the basic laws for mass and charge transfer (given in section 3.1). The fitting of this mathematical expression provides the values of characteristic parameters (rate constants, transfer coefficients, diffusion coefficients), resulting in a quantitative reaction mechanism for

The basic approach is usually the same: a *cost function* that measures the agreement between the data and the model with a particular choice of parameters is designed. The cost function usually defines a distance between the experimental data and the model and is conventionally arranged so that small values represent close agreement. The parameters of the model are then adjusted to achieve a minimum in the cost function. This is schematically illustrated in Figure 2. It shows an imaginary experiment, modeled with the well-known Butler-Volmer equation, which describes the rate of an electrochemical reaction under charge transport control. By changing the values of the transfer coefficients and the rate constants, the distance between model and experiment varies for each data point. The cost function takes all these distances into account. The *estimates* or *best-fit-parameters* are the arguments that minimize the cost

*mXi* <sup>=</sup> 0.620*D*2/3

(*t*) = 0 (6)

*Xi <sup>ν</sup>*−1/6*ω*1/2 (7)

*<sup>i</sup>* (8)

*Xi* is the concentration of the electroactive species. In electrode kinetics, concentrations are usually employed rather than activities and thus the rate constants include the product of activity coefficients. *Kj* and *K*� *<sup>j</sup>* are the potential dependent rate constants of the partial reaction, given by:

$$K\_{\circ} \text{ or } K\_{\circ}^{\prime} = k\_{ox,j} \exp\left(\frac{\mathfrak{a}\_{ox,j} n\_{j} FE(t)}{RT}\right) \tag{3}$$

$$K\_{\circ} \text{ or } K\_{\circ}^{\prime} = k\_{red,j} \exp\left(\frac{-\mathfrak{a}\_{red,j} n\_{j} FE(t)}{RT}\right)$$

with *αox*,*<sup>j</sup>* + *αred*,*<sup>j</sup>* = 1. F is the Faraday constant (96485 C/mol), R is the ideal gas constant (8.32 J/mol K) and T is the absolute temperature (*K*). *kox*,*<sup>j</sup>* is the rate constant of the partial reaction *j* in the sense of the oxidation and *kred*,*<sup>j</sup>* is the one in the sense of the reduction. *αox*,*<sup>j</sup>* and *αred*,*<sup>j</sup>* are the transfer coefficients in the sense of the oxidation and reduction respectively. They are supposed to be independent of the electrode potential (Diard et al., 1996). If adsorbed species are involved in the electrochemical reaction the rate constants *Kj* and *K*� *<sup>j</sup>* may depend on the coverage.

• the relations that connect the faradaic current density to the rates of the partial electrochemical reactions:

$$\dot{a}\_f = F \sum\_{j=1}^{P} s\_{\varepsilon j} v\_j \tag{4}$$

where *sej* = ±*nj*. The + or − sign is fixed by the following convention: an oxidation current is counted positive, a reduction current negative.

• the relations expressing the transformation rate of the electroactive species at the interface electrode/electrolyte, and the continuity relation which expresses that this transformation rate is equal to the interfacial mass transport flux of the species *Xi*:

$$\begin{aligned} v\_{X\_i} &= f\_{X\_i} \\ &= -m\_{X\_i} [X\_i^\* - X\_i(0)] \\ &= \sum\_{j=1}^P s\_{ij} v\_j \\ i &= 1, \dots, N\_v \end{aligned} \tag{5}$$

<sup>1</sup> To simplify the notation the symbol *Xi* is used to describe the chemical species in a reaction and its concentration *Xi* = [*Xi*]

or

6 Electrochemical Cells

surface concentrations *Xi* of the adsorbed species1. Under steady-state conditions and using an RDE, the general system of equations describing this electrochemical systems is given by:

*Xi* is the concentration of the electroactive species. In electrode kinetics, concentrations are usually employed rather than activities and thus the rate constants include the product

*<sup>X</sup>pij*

*αox*,*jnjFE*(*t*)

−*αred*,*jnjFE*(*t*) *RT* )

*<sup>j</sup>* are the potential dependent rate constants of the partial

*<sup>i</sup>* (2)

*RT* ) (3)

*sejvj* (4)

*vXi* = *JXi* (5)

*N* ∏ *i*=1 *Xrij <sup>i</sup>* − *K*� *j N* ∏ *i*=1

*<sup>j</sup>* = *kox*,*<sup>j</sup>* exp (

*<sup>j</sup>* = *kred*,*<sup>j</sup>* exp (

with *αox*,*<sup>j</sup>* + *αred*,*<sup>j</sup>* = 1. F is the Faraday constant (96485 C/mol), R is the ideal gas constant (8.32 J/mol K) and T is the absolute temperature (*K*). *kox*,*<sup>j</sup>* is the rate constant of the partial reaction *j* in the sense of the oxidation and *kred*,*<sup>j</sup>* is the one in the sense of the reduction. *αox*,*<sup>j</sup>* and *αred*,*<sup>j</sup>* are the transfer coefficients in the sense of the oxidation and reduction respectively. They are supposed to be independent of the electrode potential (Diard et al., 1996). If adsorbed species are involved in the electrochemical reaction the rate constants

• the relations that connect the faradaic current density to the rates of the partial

*P* ∑ *j*=1

where *sej* = ±*nj*. The + or − sign is fixed by the following convention: an oxidation

[*X*∗

*<sup>i</sup>* − *Xi*(0)]

• the relations expressing the transformation rate of the electroactive species at the interface electrode/electrolyte, and the continuity relation which expresses that this transformation

= −*mXi*

*i* = 1, . . . , *Nv*

<sup>1</sup> To simplify the notation the symbol *Xi* is used to describe the chemical species in a reaction and its

= *P* ∑ *j*=1 *sijvj*

*if* = *F*

*vj* = *Kj*

*Kj* or *K*�

*Kj* or *K*�

*<sup>j</sup>* may depend on the coverage.

current is counted positive, a reduction current negative.

rate is equal to the interfacial mass transport flux of the species *Xi*:

• the rate (*vj*) expressions for each partial reaction:

of activity coefficients. *Kj* and *K*�

reaction, given by:

*Kj* and *K*�

electrochemical reactions:

concentration *Xi* = [*Xi*]

$$\begin{aligned} v\_{X\_i}(t) &= 0\\ \dot{\mathbf{r}} &= N\_{\overline{v}} + \mathbf{1}\_i \dots \mathbf{N} \end{aligned} \tag{6}$$

were *JXi* is the molecular flux (expressed in mol/m2 s) of *Xi*, equal to the number of moles of *Xi* going per unit of time across a unit plane, perpendicularly oriented to the flow of the species. *X*∗ *<sup>i</sup>* is the bulk concentration of species *Xi*. *mXi* is the mass transport rate constant for diffusion and convection of species *Xi*, given by:

$$m\_{X\_l} = 0.620 D\_{X\_l}^{2/3} v^{-1/6} w^{1/2} \tag{7}$$

with *DXi* the diffusion coefficient of *Xi*, *ν* the kinematic viscosity and *ω* the rotation speed of the RDE.

At the maximum diffusion rate, the concentration at the electrode is zero and the limiting current is achieved. The Levich equation provides an expression of the limiting current as a function of the mass transport rate constant and the bulk concentration:

$$I\_{\rm lim} = 0.620 n FSD\_{X\_l}^{2/3} \nu^{-1/6} \omega^{1/2} X\_l^\* \tag{8}$$

with *S* the electrode surface and *n* the number of electrons exchanged in the reaction. This equation thus applies to the totally mass-transfer-limited condition at the RDE.

#### **3.2 The fitting procedure**

Many comprehensive textbooks about parameter estimation and minimization algorithms are available. The development of the fitting procedure for the analytical modeling of electrochemical reactions is founded on a few of them (Fletcher, 1980; Kelley, 1999; Norton, 1986; Pintelon & Schoukens, 2001; Press et al., 1988; Sorenson, 1980).

Given a set of observations, one often wants to condense and summarize the data by fitting it to a 'model' that depends on adjustable parameters. In this work the 'model' is the current-potential relation, describing the polarization curve, which is derived from the basic laws for mass and charge transfer (given in section 3.1). The fitting of this mathematical expression provides the values of characteristic parameters (rate constants, transfer coefficients, diffusion coefficients), resulting in a quantitative reaction mechanism for the electrochemical reaction.

The basic approach is usually the same: a *cost function* that measures the agreement between the data and the model with a particular choice of parameters is designed. The cost function usually defines a distance between the experimental data and the model and is conventionally arranged so that small values represent close agreement. The parameters of the model are then adjusted to achieve a minimum in the cost function. This is schematically illustrated in Figure 2. It shows an imaginary experiment, modeled with the well-known Butler-Volmer equation, which describes the rate of an electrochemical reaction under charge transport control. By changing the values of the transfer coefficients and the rate constants, the distance between model and experiment varies for each data point. The cost function takes all these distances into account. The *estimates* or *best-fit-parameters* are the arguments that minimize the cost function.

generates a sequence of points, noted as *<sup>θ</sup>*(1), *<sup>θ</sup>*(2), *<sup>θ</sup>*(3), . . ., or {*θ*(*k*)}. The iteration mechanism

<sup>29</sup> Modeling and Quantification of Electrochemical Reactions in RDE (Rotating Disk Electrode) and IRDE (Inverted Rotating Disk Electrode) Based Reactors

where *�<sup>i</sup>* is the relative error in the minimization algorithm, which is defined by the user. That means that the fitting stops when the parameter values from one iteration to another do not

Using singular value decomposition (SVD) techniques, the Gauss-Newton algorithm can be solved without forming the product *J*(*θ*(*k*))*<sup>T</sup> J*(*θ*(*k*)) so that more complex problems can be solved because the numerical errors are significantly reduced. By SVD the matrix *J* is transformed into the product *J* = *U*Σ*V<sup>T</sup>* with *U* and *V* orthogonal matrices: *UTU* = *I* and *VTV* = *VV<sup>T</sup>* = *I*. Σ is a diagonal matrix with the singular values on the diagonal. This leads

Under quite general assumptions on the noise, some regularity conditions on the model *F*(*u*0(*l*), *θ*) and the excitation (choice of *u*0(*l*)), consistency <sup>2</sup> of the least squares estimator is proven. Asymptotically (for the number of data points going to infinity) the covariance

The Levenberg-Marquardt algorithm is a popular alternative to the Gauss-Newton method. This method is often considered to be the best type of method for non-linear least squares

)*TCy J*(*θ*(*k*)

))−<sup>1</sup> *J*(*θ*(*k*)

problems, but the rate of convergence can be slow. The iterative process is given by:

⎛

*s*1 *s*2

1+*λ*<sup>2</sup> 0 0 0 *<sup>s</sup>*<sup>2</sup> *s*2 2+*λ*<sup>2</sup> <sup>0</sup> 0 0 ...

with *s* the singular values and *λ* is called the Levenberg-Marquardt factor. As starting value

Some numerical aspects of these methods are also worth mentioning. Whether the Gauss-Newton or the Levenberg-Marquardt algorithm is used, the expression for equation

*<sup>θ</sup>* is strongly consistent if it converges almost surely to *<sup>θ</sup>*0: a.s. lim*Ndp*→<sup>∞</sup> <sup>ˆ</sup>

<sup>100</sup> is chosen. If the value of the cost function decreases after performing an iteration, a new iteration is performed with *λnew* = *λold* ∗ 0.4. If the cost function increases, *λnew* =

⎞

⎟⎟⎠

⎜⎜⎝

(*k*) *i*


*<sup>δ</sup>* <sup>=</sup> <sup>−</sup>*V*Σ−1*UTe* (11)

)*<sup>T</sup> J*(*θ*(*k*)

*UTe* (13)

))−<sup>1</sup> (12)

*θ*(*Ndp*) = *θ*<sup>0</sup> ,

)(*J*(*θ*(*k*)

stops when a convergence test is satisfied, and the following criterion is used here:

*θ* (*k*+1) *i*


matrix *CLS* of the estimated model parameters is given by:

with *Cy* = *cov*{*ny*} the covariance matrix of the noise.

1. Solve (*J*(*θ*(*k*))*<sup>T</sup> <sup>J</sup>*(*θ*(*k*)) + *<sup>λ</sup>*<sup>2</sup> *<sup>I</sup>*)*δ*(*k*) <sup>=</sup> <sup>−</sup>*J*(*θ*(*k*))*Te*(*θ*(*k*))

*λold* ∗ 10 is chosen and the old value of *θ* is maintained.

)*<sup>T</sup> J*(*θ*(*k*)

*δ* = −*V*

*CLS* = (*J*(*θ*(*k*)

change more than a low *�<sup>i</sup>* value.

to the following expression for *δ*:

2. Set *θ*(*k*+1) = *θ*(*k*) + *δ*(*k*)

to be solved is always written as:

with *θ*<sup>0</sup> the true (unknown) value *θ*.

or using SVD:

*λ* = *<sup>s</sup>*<sup>1</sup>

<sup>2</sup> An estimator ˆ

Fig. 2. The cost function as a measure of the agreement between the experiment and the model.

In most data fitting problems the match between the model and the measurements is quantified by a least squares cost function. Consider a multiple input, single output system modeled by *<sup>y</sup>*0(*l*) = *<sup>F</sup>*(*u*0(*l*), *<sup>θ</sup>*) with *<sup>l</sup>* the measurement index, *<sup>y</sup>*(*l*) <sup>∈</sup> **<sup>R</sup>**, *<sup>u</sup>*(*l*) <sup>∈</sup> **<sup>R</sup>**1×*M*, and *<sup>θ</sup>* <sup>∈</sup> **<sup>R</sup>***nθ*×<sup>1</sup> the parameter vector. The aim is to estimate the parameters from noisy observations of the output of the system: *y*(*l*) = *y*0(*l*) + *ny*(*l*). This is done by minimizing the sum of the squared errors, *e*(*l*) = *y*(*l*) − *y*0(*l*):

$$\hat{\theta}\_{LS} = \underset{\theta}{\text{arg min }} V\_{LS}(\theta) \text{ with } V\_{LS} = \frac{1}{2} \sum\_{l=1}^{N\_{d\rho}} \varepsilon^2(l) \tag{9}$$

with *Ndp* the number of data points. *arg min<sup>θ</sup>* stands for: "find the argument *θ* that minimizes ...".

In this work polarization curves will be modeled. The current response (= output *y*) to a linearly changing potential (= input *u*0) is measured in *Ndp* different points (with index *l*) and the transfer coefficients, rate constants and diffusion coefficients (= *θ*) need to be estimated. The function *F* is a theoretical expression which relates the current to the potential and the model parameters.

The Gauss-Newton algorithm is very well suited to deal with the least squares minimization problem. The numerical solution is found by applying the following iterative process:


with *J*(*θ*) the *Ndp* × *n<sup>θ</sup>* Jacobian matrix. This is the matrix of all first-order partial derivatives of the error *e*. To improve the numerical conditioning of the expression 1 of the iterative process, the Jacobian can be scaled by multiplying each column with the corresponding parameter value from the previous iteration. As the Gauss-Newton method is an iterative method, it generates a sequence of points, noted as *<sup>θ</sup>*(1), *<sup>θ</sup>*(2), *<sup>θ</sup>*(3), . . ., or {*θ*(*k*)}. The iteration mechanism stops when a convergence test is satisfied, and the following criterion is used here:

$$|\left.\frac{\theta\_i^{(k+1)} - \theta\_i^{(k)}}{\theta\_i^{(k+1)}}\right| \le \epsilon\_i \text{ } \forall i \tag{10}$$

where *�<sup>i</sup>* is the relative error in the minimization algorithm, which is defined by the user. That means that the fitting stops when the parameter values from one iteration to another do not change more than a low *�<sup>i</sup>* value.

Using singular value decomposition (SVD) techniques, the Gauss-Newton algorithm can be solved without forming the product *J*(*θ*(*k*))*<sup>T</sup> J*(*θ*(*k*)) so that more complex problems can be solved because the numerical errors are significantly reduced. By SVD the matrix *J* is transformed into the product *J* = *U*Σ*V<sup>T</sup>* with *U* and *V* orthogonal matrices: *UTU* = *I* and *VTV* = *VV<sup>T</sup>* = *I*. Σ is a diagonal matrix with the singular values on the diagonal. This leads to the following expression for *δ*:

$$\delta = -V\Sigma^{-1}\mathcal{U}^T e\tag{11}$$

Under quite general assumptions on the noise, some regularity conditions on the model *F*(*u*0(*l*), *θ*) and the excitation (choice of *u*0(*l*)), consistency <sup>2</sup> of the least squares estimator is proven. Asymptotically (for the number of data points going to infinity) the covariance matrix *CLS* of the estimated model parameters is given by:

$$\mathcal{C}\_{LS} = (f(\theta^{(k)})^T f(\theta^{(k)}))^{-1} f(\theta^{(k)})^T \mathcal{C}\_{\mathcal{Y}} f(\theta^{(k)}) (f(\theta^{(k)})^T f(\theta^{(k)}))^{-1} \tag{12}$$

with *Cy* = *cov*{*ny*} the covariance matrix of the noise.

The Levenberg-Marquardt algorithm is a popular alternative to the Gauss-Newton method. This method is often considered to be the best type of method for non-linear least squares problems, but the rate of convergence can be slow. The iterative process is given by:


or using SVD:

8 Electrochemical Cells

potential

Fig. 2. The cost function as a measure of the agreement between the experiment and the

In most data fitting problems the match between the model and the measurements is quantified by a least squares cost function. Consider a multiple input, single output system modeled by *<sup>y</sup>*0(*l*) = *<sup>F</sup>*(*u*0(*l*), *<sup>θ</sup>*) with *<sup>l</sup>* the measurement index, *<sup>y</sup>*(*l*) <sup>∈</sup> **<sup>R</sup>**, *<sup>u</sup>*(*l*) <sup>∈</sup> **<sup>R</sup>**1×*M*, and *<sup>θ</sup>* <sup>∈</sup> **<sup>R</sup>***nθ*×<sup>1</sup> the parameter vector. The aim is to estimate the parameters from noisy observations of the output of the system: *y*(*l*) = *y*0(*l*) + *ny*(*l*). This is done by minimizing the sum of the

*VLS*(*θ*) with *VLS* <sup>=</sup> <sup>1</sup>

with *Ndp* the number of data points. *arg min<sup>θ</sup>* stands for: "find the argument *θ* that minimizes

In this work polarization curves will be modeled. The current response (= output *y*) to a linearly changing potential (= input *u*0) is measured in *Ndp* different points (with index *l*) and the transfer coefficients, rate constants and diffusion coefficients (= *θ*) need to be estimated. The function *F* is a theoretical expression which relates the current to the potential and the

The Gauss-Newton algorithm is very well suited to deal with the least squares minimization problem. The numerical solution is found by applying the following iterative process:

with *J*(*θ*) the *Ndp* × *n<sup>θ</sup>* Jacobian matrix. This is the matrix of all first-order partial derivatives of the error *e*. To improve the numerical conditioning of the expression 1 of the iterative process, the Jacobian can be scaled by multiplying each column with the corresponding parameter value from the previous iteration. As the Gauss-Newton method is an iterative method, it

2

*Ndp* ∑ *l*=1 *e*

<sup>2</sup>(*l*) (9)

distance between BV 1 and experiment in data point i

 Butler-Volmer 1 Butler-Volmer 2 experimental data

current

squared errors, *e*(*l*) = *y*(*l*) − *y*0(*l*):

ˆ

1. Solve *<sup>J</sup>*(*θ*(*k*))*<sup>T</sup> <sup>J</sup>*(*θ*(*k*))*δ*(*k*) <sup>=</sup> <sup>−</sup>*J*(*θ*(*k*))*Te*(*θ*(*k*))

*θLS* = *arg min θ*

model.

...".

model parameters.

2. Set *θ*(*k*+1) = *θ*(*k*) + *δ*(*k*)

distance between BV 2 and experiment in data point i

$$\delta = -V \begin{pmatrix} \frac{s\_1}{s\_1^2 + \lambda^2} & 0 & 0\\ 0 & \frac{s\_2}{s\_2^2 + \lambda^2} & 0\\ 0 & 0 & \ddots \end{pmatrix} \mathcal{U}^T e \tag{13}$$

with *s* the singular values and *λ* is called the Levenberg-Marquardt factor. As starting value *λ* = *<sup>s</sup>*<sup>1</sup> <sup>100</sup> is chosen. If the value of the cost function decreases after performing an iteration, a new iteration is performed with *λnew* = *λold* ∗ 0.4. If the cost function increases, *λnew* = *λold* ∗ 10 is chosen and the old value of *θ* is maintained.

Some numerical aspects of these methods are also worth mentioning. Whether the Gauss-Newton or the Levenberg-Marquardt algorithm is used, the expression for equation to be solved is always written as:

<sup>2</sup> An estimator ˆ *<sup>θ</sup>* is strongly consistent if it converges almost surely to *<sup>θ</sup>*0: a.s. lim*Ndp*→<sup>∞</sup> <sup>ˆ</sup> *θ*(*Ndp*) = *θ*<sup>0</sup> , with *θ*<sup>0</sup> the true (unknown) value *θ*.

$$X\delta = e\tag{14}$$

from repeated measurements is used. It can be shown that in this case the covariance matrix

<sup>31</sup> Modeling and Quantification of Electrochemical Reactions in RDE (Rotating Disk Electrode) and IRDE (Inverted Rotating Disk Electrode) Based Reactors

(*J*(*θ*(*k*)

where *M* is the number of polarization curves that are measured. At least 11 experiments

*CWLS* <sup>=</sup> *<sup>V</sup>*Σ−2*V<sup>T</sup> <sup>M</sup>* <sup>−</sup> <sup>1</sup>

In this work the Gauss-Newton and Levenberg-Marquardt methods are implemented. The fitting starts with the Gauss-Newton method, but when the cost function is no longer decreasing after an iteration, it switches to Levenberg-Marquardt. The reason for starting with the Gauss-Newton method is that this method usually converges more rapidly than the

After the determination of the best-fit parameter values, the validity of the selected model should be tested: does this model describe the available data properly or are there still indications that a part of the data is not explained by the model, indicating remaining model errors? We need the means to assess whether or not the model is appropriate, that is, we need

With the best-fit-parameters, the modeled data is calculated and compared with the experimental data. Also, the difference between the experimental and calculated data is plotted. If the model describes the experiment appropriately, this difference should lie in the confidence band. This interval is defined by ± two times the experimental standard deviation of the current, calculated from a series of repeated experiments. When performing multiple

If the statistical evaluation does not show a good agreement between model and experiments, the model cannot be accepted indicating that any of the modeling steps is not well designed. In that case, we have to go back to the previous modeling steps: the experimental study, the proposed model (reaction or transport model) or/and the fitting procedure, and check their performance and the assumptions and simplifications that have been made to develop every step of the modeling methodology. Adjustments or possible alternatives must then be

An important message is that the fitting of parameters is not the end-all of parameter estimation. To be genuinely useful, a fitting procedure should provide (i) parameters, (ii) error estimates on the parameters, and (iii) a statistical measure of goodness-of-fit. When the third item suggests that the model is an unlikely match to the data, then items (i) and (ii) are probably worthless. Unfortunately, many practitioners of parameter estimation never proceed beyond item (i). They deem a fit acceptable if a graph of data and model "look good". This

)*TC*−<sup>1</sup> *<sup>y</sup> <sup>J</sup>*(*θ*(*k*)

))−<sup>1</sup> (21)

*<sup>M</sup>* <sup>−</sup> <sup>5</sup> (22)

of the estimates is given by (Schoukens et al., 1997):

After SVD of *C*−<sup>1</sup>

Levenberg-Marquardt algorithm.

**3.3 A statistical evaluation**

proposed.

*CWLS* <sup>=</sup> *<sup>M</sup>* <sup>−</sup> <sup>1</sup>

have to be performed to use this expression (Schoukens et al., 1997).

to test the *goodness-of-fit* against some useful statistical standard.

experiments, 95 % of the experiments are expected to fall in this interval.

approach is known as *chi-by-eye* and should definitively be avoided.

*<sup>y</sup> J* this equation becomes:

*M* − 5

Consequently, a change Δ*e* of *e*, also causes *δ* to change by Δ*δ*. The problem is said to be *well conditioned* if a small change of *e* results in a small change of *δ*. If not, the problem is *ill-conditioned*. It can be shown that

$$\frac{||\,\Delta\delta\;||\_2}{||\,\delta\;||\_2} \le \frac{s\_{\max}}{s\_{\min}} \frac{||\,\Delta\varepsilon\;||\_2}{||\,\varepsilon\;||\_2} \tag{15}$$

with || ||<sup>2</sup> the 2-norm given by:

$$|||\propto|||\_2 = \sum\_{i} \mathbf{x}\_i^2\tag{16}$$

and *κ* = *smax smin* , the ratio of the largest singular value to the smallest, is called the condition number of *X*. If *κ* is large, the problem is ill-conditioned.

In Eq. (9) all measurements are equally weighted. In many problems it is desirable to put more emphasis on one measurement with respect to the other. This is done to make the difference between measurement and model smaller in some regions. If the covariance matrix of the noise is known, then it seems logical to suppress measurements with high uncertainty and to emphasize those with low uncertainty. In practice it is not always clear what weighting should be used. If it is, for example, known that model errors are present, then the user may prefer to put in a dedicated weighting in order to keep the model errors small in some specific operation regions instead of using the weighting dictated by the covariance matrix.

In general, the weighted least squares estimator ˆ *θWLS* is:

$$\hat{\theta}\_{WLS} = \underset{\theta}{\text{arg min}} \, V\_{WLS}(\theta) \text{ with } V\_{WLS} = \frac{1}{2} e(\theta)^T W e(\theta) \tag{17}$$

where *<sup>W</sup>* <sup>∈</sup> **<sup>R</sup>***N*×*<sup>N</sup>* is a symmetric positive definite weighting matrix. All the remarks on the numerical aspects of the least squares estimator are also valid for the weighted least squares. This can be easily understood by applying the following transformation: *�* = *Se* with *STS* = *W* so that *VWLS* = *�T�*, which is a least squares estimator in the transformed variables. This also leads to the following Gauss-Newton algorithm to minimize the cost function:

$$\theta^{(k+1)} = \theta^{(k)} + \delta^{(k)}, \text{with } J(\theta^{(k)})^T W \!\!/ (\theta^{(k)}) \delta^{(k)} = -J(\theta^{(k)})^T W e(\theta^{(k)}) \tag{18}$$

Eq. (12) is generalized to (noticing that *W<sup>T</sup>* = *W*):

$$\mathcal{L}\_{\text{WLS}} = (f(\theta^{(k)})^T \mathcal{W} f(\theta^{(k)}))^{-1} f(\theta^{(k)})^T \mathcal{W} \mathcal{C}\_y \mathcal{W} f(\theta^{(k)}) (f(\theta^{(k)})^T \mathcal{W} f(\theta^{(k)}))^{-1} \tag{19}$$

By choosing *W* = *C*−<sup>1</sup> *<sup>y</sup>* the expression simplifies to:

$$\mathcal{C}\_{WLS} = (J(\theta^{(k)})^T \mathcal{C}\_y^{-1} J(\theta^{(k)}))^{-1} \tag{20}$$

It can be shown that, among all possible positive definite choices for W, the best one is *W* = *C*−<sup>1</sup> *<sup>y</sup>* since this minimizes the covariance matrix (Pintelon & Schoukens, 2001). These expressions depend on the covariance matrix of the noise *Cy*. In practice this knowledge should be obtained from measurements. In this work a sample covariance matrix obtained from repeated measurements is used. It can be shown that in this case the covariance matrix of the estimates is given by (Schoukens et al., 1997):

$$\mathbb{C}\_{WLS} = \frac{M-1}{M-5} (f(\theta^{(k)})^T \mathbb{C}\_y^{-1} f(\theta^{(k)}))^{-1} \tag{21}$$

where *M* is the number of polarization curves that are measured. At least 11 experiments have to be performed to use this expression (Schoukens et al., 1997).

After SVD of *C*−<sup>1</sup> *<sup>y</sup> J* this equation becomes:

$$\mathbb{C}\_{WLS} = V\Sigma^{-2} V^T \frac{M-1}{M-5} \tag{22}$$

In this work the Gauss-Newton and Levenberg-Marquardt methods are implemented. The fitting starts with the Gauss-Newton method, but when the cost function is no longer decreasing after an iteration, it switches to Levenberg-Marquardt. The reason for starting with the Gauss-Newton method is that this method usually converges more rapidly than the Levenberg-Marquardt algorithm.

### **3.3 A statistical evaluation**

10 Electrochemical Cells

Consequently, a change Δ*e* of *e*, also causes *δ* to change by Δ*δ*. The problem is said to be *well conditioned* if a small change of *e* results in a small change of *δ*. If not, the problem is

> || Δ*e* ||<sup>2</sup> || *e* ||<sup>2</sup>

≤ *smax smin*


In Eq. (9) all measurements are equally weighted. In many problems it is desirable to put more emphasis on one measurement with respect to the other. This is done to make the difference between measurement and model smaller in some regions. If the covariance matrix of the noise is known, then it seems logical to suppress measurements with high uncertainty and to emphasize those with low uncertainty. In practice it is not always clear what weighting should be used. If it is, for example, known that model errors are present, then the user may prefer to put in a dedicated weighting in order to keep the model errors small in some specific

operation regions instead of using the weighting dictated by the covariance matrix.

also leads to the following Gauss-Newton algorithm to minimize the cost function:

))−<sup>1</sup> *J*(*θ*(*k*)

, with *J*(*θ*(*k*)

*i x*2

*smin* , the ratio of the largest singular value to the smallest, is called the condition

*θWLS* is:

2

)*δ*(*k*) <sup>=</sup> <sup>−</sup>*J*(*θ*(*k*)

)(*J*(*θ*(*k*)

*VWLS*(*θ*) with *VWLS* <sup>=</sup> <sup>1</sup>

)*TW J*(*θ*(*k*)

)*TC*−<sup>1</sup> *<sup>y</sup> <sup>J</sup>*(*θ*(*k*)

*<sup>y</sup>* since this minimizes the covariance matrix (Pintelon & Schoukens, 2001). These expressions depend on the covariance matrix of the noise *Cy*. In practice this knowledge should be obtained from measurements. In this work a sample covariance matrix obtained

It can be shown that, among all possible positive definite choices for W, the best one is

)*TWCyW J*(*θ*(*k*)

where *<sup>W</sup>* <sup>∈</sup> **<sup>R</sup>***N*×*<sup>N</sup>* is a symmetric positive definite weighting matrix. All the remarks on the numerical aspects of the least squares estimator are also valid for the weighted least squares. This can be easily understood by applying the following transformation: *�* = *Se* with *STS* = *W* so that *VWLS* = *�T�*, which is a least squares estimator in the transformed variables. This


number of *X*. If *κ* is large, the problem is ill-conditioned.

In general, the weighted least squares estimator ˆ

*θWLS* = *arg min θ*

ˆ

*θ*(*k*+1) = *θ*(*k*) + *δ*(*k*)

*CWLS* = (*J*(*θ*(*k*)

By choosing *W* = *C*−<sup>1</sup>

*W* = *C*−<sup>1</sup>

Eq. (12) is generalized to (noticing that *W<sup>T</sup>* = *W*):

)*TW J*(*θ*(*k*)

*<sup>y</sup>* the expression simplifies to:

*CWLS* = (*J*(*θ*(*k*)

*ill-conditioned*. It can be shown that

with || ||<sup>2</sup> the 2-norm given by:

and *κ* = *smax*

*Xδ* = *e* (14)

*<sup>i</sup>* (16)

*e*(*θ*)*TWe*(*θ*) (17)

)*TWe*(*θ*(*k*)

)*TW J*(*θ*(*k*)

))−<sup>1</sup> (20)

) (18)

))−<sup>1</sup> (19)

(15)

After the determination of the best-fit parameter values, the validity of the selected model should be tested: does this model describe the available data properly or are there still indications that a part of the data is not explained by the model, indicating remaining model errors? We need the means to assess whether or not the model is appropriate, that is, we need to test the *goodness-of-fit* against some useful statistical standard.

With the best-fit-parameters, the modeled data is calculated and compared with the experimental data. Also, the difference between the experimental and calculated data is plotted. If the model describes the experiment appropriately, this difference should lie in the confidence band. This interval is defined by ± two times the experimental standard deviation of the current, calculated from a series of repeated experiments. When performing multiple experiments, 95 % of the experiments are expected to fall in this interval.

If the statistical evaluation does not show a good agreement between model and experiments, the model cannot be accepted indicating that any of the modeling steps is not well designed. In that case, we have to go back to the previous modeling steps: the experimental study, the proposed model (reaction or transport model) or/and the fitting procedure, and check their performance and the assumptions and simplifications that have been made to develop every step of the modeling methodology. Adjustments or possible alternatives must then be proposed.

An important message is that the fitting of parameters is not the end-all of parameter estimation. To be genuinely useful, a fitting procedure should provide (i) parameters, (ii) error estimates on the parameters, and (iii) a statistical measure of goodness-of-fit. When the third item suggests that the model is an unlikely match to the data, then items (i) and (ii) are probably worthless. Unfortunately, many practitioners of parameter estimation never proceed beyond item (i). They deem a fit acceptable if a graph of data and model "look good". This approach is known as *chi-by-eye* and should definitively be avoided.

0.04 *μ*m Al2O3 paste (Struers) for the Au electrode; 2) ultrasonic rinsing with deionized water followed by degreasing with chloroform, also in an ultrasonic bath (Elma model T470/H); 3) four cyclic voltammograms are performed before each experiment in order to remove oxide and trace contaminants from the metallic surface (Robertson et al., 1988). During the cyclic voltammetry measurements, the potential is swept between -0.35 V to 1.45 V vs Ag/AgCl at a scan velocity of 50 mV/s on the Au electrode (hexaammineruthenium (III)/(II) system) and between +0.55 and -0.45 V vs Ag/AgCl at a scan velocity of 10 mV/s, on the Pt electrode

<sup>33</sup> Modeling and Quantification of Electrochemical Reactions in RDE (Rotating Disk Electrode) and IRDE (Inverted Rotating Disk Electrode) Based Reactors

The simplest electrochemical reactions, which can be found among the different kinds of electrode processes, are those where electrons are exchanged across the interface by flipping oxidation states of transition metal ions in the electrolyte adjacent to the electrode surface (Bamford & Compton, 1986), i.e. an ET (electron transfer) mechanism. The electrode acts as the source or sink of electrons for the redox reaction and is supposed to be inert. The reduction of ferricyanide to ferrocyanide (Angell & Dickinson, 1972; Bamford & Compton, 1986; Bruce

(ferri/ferrocyanide system).

**4.2 Modeling the ferri/ferrocyanide redox reaction**

**4.2.1 Results of the experimental study**

The resulting voltammograms for the Fe(CN)−<sup>3</sup>

−4

−2

0

I (A)

at 2000 rpm.

2

4

<sup>6</sup> x 10−4

et al., 1994; Iwasita et al., 1983) is an example of such a mechanism:

*Fe*(*CN*)−<sup>3</sup>

<sup>6</sup> + *e*

A set of equivalent polarization curves of the ferri/ferrocyanide reaction is used as the experimental data, which is the input to the model methodology. As advised in (Tourwé et al., 2006), 11 voltammograms are measured under identical conditions for every rotation speed.

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 −6

E (V vs NHE)

Fig. 3. Voltammograms of the reduction/oxidation of 0.005 M ferri/ferrocyanide in 1 M KCl,

The analytical expression is derived for the reduction/oxidation of ferri/ferrocyanide (reaction (23)). The basic equations that describe this mechanism, when studied under

**4.2.2 The analytical expression for the current as a function of the potential**

<sup>−</sup> *Fe*(*CN*)−<sup>4</sup>

<sup>6</sup> /Fe(CN)−<sup>3</sup>

<sup>6</sup> (23)

<sup>6</sup> system are shown in Figure 3.
