1.1 Why model voltammetry?

Voltammetry, the measurement of current as a function of applied voltage, is a sensitive technique for the investigation of chemical kinetics - the current contains quantitative mechanistic information, yet the experimental design is relatively simple and economical. In addition it is extremely useful for analytical chemistry, finding numerous applications in sensors1, since it provides a direct interface between chemistry and electronics.

The difficult part of a voltammetric experiment is extracting the chemical information from the current-voltage curve. Apart from very simplistic analysis (such as the visual identification of the relative number of electrons transferred from the height of voltammetric 'waves'), the measured current cannot be directly interpreted. As with many experiments, a model must be constructed to predict the current for a given set of conditions and a postulated chemical mechanism. An electrochemical model is concerned with the concentration distributions of chemical species (and possibly the potential distribution if an excess of background electrolyte is not used). If the concentration distributions of all the chemical species can be simulated, the current flowing at the 'working electrode' may be calculated fairly easily.

The majority of the work conducted in this area2 has been done for systems where the concentration distribution can be described in one spatial dimension. Where electrode geometries giving rise to two-dimensional distributions have been addressed, only simple 'model' reactions have been simulated in most cases. In the few cases where simple homogeneous chemical processes have been simulated, (explicit) approximations were necessary to incorporate the chemical steps, which break down for large rate constants.

The work presented in this thesis is concerned with the development of more accurate, efficient and powerful simulation methods for voltammetric experiments at hydrodynamic electrodes and microelectrodes. This is achieved by applying efficient methods for the solution of sparse linear systems equations including the multigrid and preconditioned Krylov subspace methods.

1.2 The four components of an electrochemical model

Three pieces of information are required to define the experimental system:

If all of these are known, the concentration distributions of the species throughout the experiment may be described mathematically as a set of coupled partial differential equations. The way these equations are perturbed during the course of the experiment and the boundary conditions required to solve them may also be deduced from these three pieces of information. The last component of the model is a method to solve this system of equations, often as a function of time as the concentration distributions evolve during the experiment.

1.3 The mechanism

An electrochemical mechanism is composed of two types of reaction:

The E and C notation3 introduced by Reinmuth and Testa in 1961 is a convenient form of mechanistic nomenclature denoting (heterogeneous) electron transfer steps as E, first-order homogenous (chemical) steps as C and second order homogenous steps as C2 etc. The Reinmuth and Testa notation does not distinguish between reactions with reversible and irreversible chemical steps, and there are examples of 'customised' notation in the literature45-6 (using Crev and Cirr).

1.3.1 Heterogeneous kinetics

Butler7 and Volmer8 derived a pair of equations for a redox couple (written here as a reduction):


which relate the forward and backward rate constants to the potential, E, a standard rate constant, koand a transfer coefficient, a (which is assumed to be potential-independent and takes values between zero and unity, describing the slope of the energy profile in the region of the transition state - a value of 0.5 is used throughout this work which implies the transition state lies midway between reactants and products) 8b:

; where Q = nF(E-Eo)/ÂT (1.2)

Such a couple is said to be quasi-reversible. A fully-reversible couple arises when ko (and hence kf and kb) tend to infinity. The ratio of the surface concentrations is then given by the Nernst equation:


The other limiting behaviour of this couple occurs when kb=0 thus it is (electrochemically) irreversible. Figure 1.1 shows how a voltammogram is shaped by electrochemical (heterogeneous) kinetics.

Figure 1.1: Voltammogram is shifted as the electron transfer becomes slower. A greater 'overpotential' is needed to overcome the kinetic barrier.

Heterogenenous processes also include dissolution/crystallisation and surface-bound chemical reactions6,9,10. Heterogeneous kinetics manifest themselves as boundary conditions on the solution of the material balance equation, namely they define the relative concentrations of species at the electrode surface.

1.3.2 Homogeneous kinetics

An electron transfer often initiates a cascade of chemical reactions by producing a reactive radical anion/cation. The mechanism can be described mathematically by a rate equation for each species. The rate law of the overall sequence is probed by an electrochemical experiment. Table 1.1 shows some common electrochemical reaction mechanisms with their rate laws. All are shown with irreversible chemical steps, apart from the CE mechanism where an initial irreversible C step is meaningless.

Table 1.1: Common electrochemical mechanisms
ProcessStepsRate Law of homogeneous steps
EC ;
EC2 ;
EC2E ;
EC' ;

1.4 The electrode geometry

The electrode geometry dictates the way material moves to and from the electrode surface, a process known as mass transport. There are three modes of mass transport:

Each of these three components may described mathematically. The effect of all three modes of mass transport may be summed giving the generalised Nernst-Planck equation:


where z (in this case) is the charge, c the concentration, E the potential, v is a velocity vector and the other symbols are defined in the glossary.

This describes how the amount of material at a given point (the concentration) varies through time due to diffusion, convection and migration. This is one major component of the electrochemical model. It depends only on the electrode geometry, the symmetry of which defines the Laplacian operator. Often a background 'supporting' electrolyte is used in excess to eliminate migration effects from the experiment. Since it is in a large excess, it is (on average) the supporting electrolyte ions that move to equalise the potential gradients and not the species being studied. If this is the case, only a convective-diffusion equation is necessary to describe the mass-transport.

1.4.1 Common geometries

Most commonly used electrodes consist of a conducting metal/glassy carbon or semiconducting surface embedded in an insulating wall. When the conducting surface is a rectangle or disc of a few millimetres, this is known as a 'planar' electrode. Diffusion to this surface is effectively planar (the effects of the edges are negligible), hence the Laplacian operator is one-dimensional:


Two other electrode geometries where diffusion occurs in only one spatial dimension are the hemispherical and hemicylindrical electrodes. The Laplacian operator for the former is:


and for the latter is:


The hemisphere can be achieved experimentally via a small drop of mercury positioned over a smaller conducting disc. A hemicylinder can be fabricated by half-embedding a wire in an insulator. Clearly these are not easily or accurately fabricated practical geometries, but these electrodes are useful for theoretical studies since the low dimensionality of the mass-transport equation under diffusion-only conditions allows it to be solved analytically. It also turns out that there are mathematical analogies between these and more robust practical geometries11.

The diffusion to a spherical electrode is convergent (it is effectively diffusion to a point) - it has the special property that the flux of material reaching the surface due to the convergent diffusion reaches a steady-state value12:


Obviously the same is true for a hemispherical electrode except that the current is half of that at the spherical electrode. The time it takes for the electrode to reach steady-state after it is 'switched on' depends on its size - longer for a larger electrode.

Closely related is the spherical Dropping Mercury Electrode (DME)13. This consists of a mercury drop, suspended from a pipette which grows in size until its weight overcomes its surface tension and falls from the pipette, where another begins to form. The current is usually recorded at the moment before the drop falls. The current does not reach a true steady-state since the drop size is constantly increasing and thus the time taken to reach steady-state is also increasing.

1.4.2 Microelectrodes

If one or more electrode dimensions are reduced from the usual few millimetres to a few microns, the electrode is known as a microelectrode (or ultramicroelectrode)14-16. The generally accepted definition2,15 seems to be that a microelectrode has at least one dimension of < 20mm. Perhaps the simplest, and certainly the most popular, microelectrode is the microdisc electrode which can be made by sealing a fine wire in glass/plastic and polishing the end. These are available commercially17 in gold, platinum and graphite with radii as small as 0.6mm. Another popular geometry is the microband. This may be fabricated in an analogous way to the microdisc - by sealing a thin layer of foil between two blocks and polishing the edge. Alternatively microbands may be fabricated by lithography18,19 (deposition of a thin film of metal onto a masked silicon wafer). This produces electrodes which are slightly elevated above the insulator on which they are deposited. Atomic Force Microscopy (AFM) and electrochemical investigations20 have shown that the approximation of treating these electrodes as if they were totally flush to the insulator introduces very little error.

The small size of microelectrodes allows some novel applications such as in vivo electroanalysis21. Such electrodes can be used to stimulate or monitor the response of individual nerve cells allowing experimental research into the central nervous system in both medicine22 and zoology23. Microelectrodes also have some useful properties due to the increased importance of mass transport to the edges of the electrode:


For both the microdisc and microband electrodes, the Laplacian operator is two-dimensional. For the microdisc:


and for the microband:


Note that the microband is still 'semi-infinite' (the edge effects are negligible) in the z co-ordinate.

1.4.3 Hydrodynamic electrodes

If convection is applied to the solution in the vicinity of the electrode, the rate of mass transport is increased - this means that faster reactions can be studied. Moreover if the fluid flow is well-defined (laminar rather than turbulent), the rate of mass transport can be controlled quantitatively by varying the solution flow rate. This is a much more convenient way of varying the rate of mass transport than changing the working electrode for one with different dimensions. It also introduces fewer random experimental errors; for example, by keeping the working electrode position constant, the Ohmic (IR) drop across the solution between counter and working electrodes remains constant. The current reaches a steady-state at hydrodynamic electrodes as fresh material is being transported to the electrode surface by convection. Rotating disc electrode (RDE)

Figure 1.2: Flow profile at a rotating disc electrode

The RDE was the first widely-used (solid) hydrodynamic electrode and is still the most popular25. A disc electrode is set in an insulating rod, which is rotated at a constant frequency in solution. The solution drag on the rotating surface results in a vortex as shown in Figure 1.2.

A major advantage when modelling a rotating disc electrode is that the convective-diffusion equation requires the simulation of only one spatial dimension:


The velocity component in the z co-ordinate is purely a function of z. For most systems vz can be approximated adequately by:

vz = CRDEz2 where CRDE = 8.032w3/2n-1/2(1.13)

In Chapter 7 the validity of this approximation is critically assessed.

The electrode and rod need to have perfect cylindrical symmetry so that the electrode does not wobble on its axis when it rotates. This does not present technical complications for electrodes with radii of several millimetres but would be a serious problem when fabricating rotating disc microelectrodes. This means that the most practical RDEs have radii of a few millimetres. Moreover the rotation speed is also limited by this constraint - a maximum practical rotation speed is about 50Hz.

One way around this problem would be to fabricate a rotating microring electrode. The diameter of the insulator inside the ring could be made large to reduce eccentricities in the rotation whist the ring itself could be made very thin to achieve high rates of mass transport. Although a microring (analogous to a microband) does not reach a steady-state under diffusion-only conditions, the convection due to rotation would ensure a steady-state response.

A few groups have successfully fabricated rotating microdisc electrodes26,27 but these have not been widely adopted, since accelerating the diffusional mass-transport rate effectively turns the RDE into a 'stationary' microdisc electrode (where the rate of diffusion is much greater than the rate of convection). Wall-jet electrode (WJE)

The wall-jet electrode28 is another disc-based arrangement. The electrode remains stationary, and a jet of solution is directed at the centre of the disc. This results in a cylindrically symmetric 'umbrella-shaped' flow profile as shown in Figure 1.3.

Figure 1.3: Flow profile at a wall-jet electrode

In terms of the rate of mass transport, the WJE is a considerable improvement over the RDE. Much higher rates of convection may be attained since the electrode remains stationary. In principle, microelectrodes may be used, though the flow jet must be of an appropriate size and precisely aligned on the centre of the disc.

Macpherson and Unwin29 have developed a 'Micro Jet Electrode' in which a jet of solution is directed at a microdisc electrode. The radius of the jet is several (4 in the published design) times larger than that of the disc (12.5 mm) so that the convection is effectively one dimensional and uniform across the electrode surface. For the flow rates employed (2x10-3 - 5x10-2 cm3s-1), the effects of convection are great enough to make edge effects negligible and the electrode effectively uniformly accessible. For smaller disc radii it is unlikely this approximation will hold. The validity of this limit is examined in Chapter 8.


Figure 1.4: (a) tubular and (b) channel electrodes The tubular and channel electrode (ChE)

These methods are very similar both conceptually and mathematically30. The solution flows down an insulating pipe which is cylindrical in the tubular system (Figure 1.4a) and rectangular in the case of the channel (Figure 1.4b). The electrode is set into the pipe so that its surface is flush20 with the surrounding insulator and thus does not perturb the laminar flow. A sufficient 'lead-in'31 is allowed before the electrode to enable a Poiseuille (parabolic) flow profile32 to develop (Figure 1.5).

Figure 1.5: Flow profile at a channel electrode

In the tubular system the electrode takes the form of a ring. In the channel, the electrode is a rectangular strip which stops short of the sides of the cell. The channel flow-cell is designed33 with the width in the y dimension much less than in the z dimension, so the flow velocity is approximately uniform across the electrode in the z co-ordinate. This reduces the convective-diffusion to a problem in two spatial dimensions, where the solution velocity in the x direction vx is given by:


where vo is the velocity down the centre of the channel. This is related to the solution volume flow rate (vf) and mean velocity (U) by:

and (1.15)

The channel and tubular electrodes have similar merits to the WJE but have a major experimental advantage over the WJE where the solution 'jet' has to be aligned on the centre of the electrode - this becomes increasingly difficult for smaller electrodes. The channel has the additional advantage that the rectangular working electrode may be fabricated using silicon-chip technology allowing microband electrodes to be made with a width as small as 1-2 microns34. Although tubular microelectrodes have been fabricated35, the tubular arrangement is less commonly used.

Recently Compton et al.36,37 modified the channel design using gas pressure-driven flow and a reduced aspect ratio cell to accelerate the flow velocity over the electrode whilst ensuring the flow profile remains laminar. This fast-flow cell results in a centre-line velocity of over 30ms-1, around 250 times faster than a conventional cell38. Cooper and Compton39 recently reviewed channel electrode voltammetry. Uniform vs. non-uniform accessibility

Under diffusion-only conditions, the concentration distribution above spherical, cylindrical and large planar electrodes is uniform across the electrode surface. The RDE has convection and diffusion acting in the same direction so, in this case also, the extent to which material escapes from the electrode surface, the diffusion layer thickness, is constant across the electrode.

In the tube/channel the convection is at its maximum at the centre of the pipe - resulting in a diffusion layer which spreads out more slowly further away from the electrode, as shown in Figure 1.6(a). In contrast, the jet velocity in the WJE is at a maximum in the centre of the electrode - the diffusion layer spreads out increasingly rapidly away from the centre of the disc as shown in Figure 1.6(b). Both the electrodes are termed as non-uniformly accessible since there is a variation of diffusion layer thickness along the electrode. However the non-uniformity is more pronounced in the case of the WJE. Under diffusion-only conditions, microdisc and microband electrodes are non-uniformly accessible since there is diffusion both into the plane of the electrode and to the edges of the electrode.

Figure 1.6: Simulated steady-state concentration profile of a simple one-electron transfer at (a) a channel electrode and (b) a wall-jet electrode. Parameters are: D=2x10-5cm2s-1, [A]bulk=1x10-6 molcm-3 and (a) Xe = 0.4cm, w = 0.4cm, 2h = 0.04cm, d = 0.6cm, vf = 0.01cm3s-1; (b) re = 0.2cm, rjet = 0.0345cm, n= 0.0089cm2s-1, vf = 0.05 cm3s-1.

The disadvantage of non-uniform accessibility is that the convective-diffusion and thus the resulting concentration distribution is inherently multi-dimensional whereas uniformly accessible electrodes may be modelled as problems in one spatial dimension and therefore more easily solved.

It is now convenient to introduce the step from the simulated concentration distribution to the current, related via the concentration gradient at the electrode surface:

I = nFJ where (1.16)

(J is the total flux reaching the electrode surface). The current may also be expressed in terms of an average flux {J}, which may be interpreted as an average linear (Nernst) diffusion layer thickness, d:

I = nFA{J} where (1.17)

where A is the electrode area and the electrode surface concentration, [C]0, is zero when the current is transport-limited. The dimensionless current (Nusselt number) is defined as:

(for the transport-limited case)(1.18)

(where le is the characteristic length of the electrode). For a non-uniformly accessible electrode, Nu is given by:


which simplifies to:


for a uniformly accessible electrode (where k is a constant).

Non-uniform accessibility has the effect of improving kinetic resolution since the electrochemical time scale varies along the electrode40,41. However this effect is not as dramatic as might be expected. Compton and Unwin42 showed that for electron-transfer and first-order homogeneous kinetics, the working curves and thus the kinetic discrimination of the ChE and RDE are virtually identical. In the concluding chapter, a short quantitative comparison of the kinetic resolution provided by common electrode geometries is presented, based on the working curves published in Chapters 6-9 of this thesis.

1.5 The experimental technique

The experimental technique controls how the mass transport and rate law are combined (and filtered) to form the overall material balance equation. Migration effects may be eliminated by the addition of supporting electrolyte; steady-state measurements eliminate the need to solve the equation in a time-dependent manner; excess substrate in a mechanism such as EC' can reduce the kinetics from second to pseudo-first order43. The material balance equations, with a given set of boundary conditions and parameters (electrode/cell dimensions, flow rate, rate constants etc.) define an I-E-t surface, which is traversed by the experimental technique.

Each experimental technique has an associated variable/scan parameter. This controls how the system is perturbed by the experiment, which can be visualised as the direction the I-E-t surface is crossed, as shown in Figure 1.7. Each of the common voltammetric techniques can be described by a slice through the I-E-t surface:

Steady-state voltammetry involves scanning the potential very slowly, or stepping the potential and waiting for the system to relax before recording the current, so that measurements are always taken at equilibrium. This corresponds to the sigmoidal I-E curve, Figure 1.7(b), (known as a 'wave' or 'voltammogram') as t => infinity.

The steady-state limiting current is the value at the top of this wave - experimentally when E is large enough to attain the limiting I value. This is often referred to as the steady-state transport-limited current, since electrochemical kinetics are effectively infinitely fast, the current is limited only by mass transport.

In the dropping mercury electrode (see 'common geometries' above), the current must be recorded before steady-state is attained. It is however recorded at a constant time so that values lie along the I-E curve shown in Figure 1.7(d). Provided this time is always constant, experiments may be conducted as if they were at steady-state. In fact the DME is often referred to as a pseudo-steady-state method.

Figure 1.7: I-E-t surface (for a spherical electrode) traversed by voltammetric experiments.

Potential step chronoamperometry takes place at a constant potential (usually large so that the limiting current is recorded). The potential is stepped from a value where no current flows to this value, and the current relaxes along an I-t curve, Figure 1.7(a), (known as a 'transient') towards its steady-state value.

In linear-sweep voltammetry (LSV) the potential is varied linearly with time (the waveform is often called a 'potential ramp'). This represents a diagonal path across the I-E-t surface, as shown in Figure 1.7(c). In contrast to the sigmoidal steady-state wave, the I-E curve (linear sweep voltammogram) goes through a maximum (known as a 'peak').

Figure 1.8: Cyclic voltammogram & potential waveform.

In cyclic voltammetry (CV), a linear sweep voltammogram is recorded and then the potential is reversed, resulting in an 'opposite' I-E curve. The pair of curves, shown in Figure 1.8, are known as a cyclic voltammogram (CV). The overall triangular E-t waveform may be applied several times, resulting in several 'cycles' being recorded - useful for recording the disappearance of reactive intermediates.

1.5.1 Obtaining kinetic information

The object of most voltammetric experiments is to record the current as a function of the rate of mass transport and/or electrochemical kinetics (which have an associated time scale), from which the rate constant(s) of the homogenous steps may be deduced. The variable parameters are tabulated for common experimental methods in Table 1.2. When comparing theory and experiment, either the whole curve may be 'fitted' to a simulated one, or fixed points may be examined and compared to theoretically predicted values.

Table 1.2 Common voltammetric techniques and the parameters which may be modulate the time scale.
MethodVariable ParameterControls
Steady-state voltammetryPotentialRate of electrochemical kinetics
Linear-sweep/cyclic voltammetryPotential and timeRate of electrochemical kinetics and mass transport
Potential-step chronoamperometryTimeRate of mass transport
Transport-limited current: wall-jetJet velocityRate of mass transport
Transport-limited current: channel flow cellFlow rate, (electrode length)Rate of mass transport
Transport-limited current: microdisc electrodeRadiusRate of mass transport
Transport-limited current: RDERotation speedRate of mass transport
† Assumes that the potential is high enough to reach the limiting current, otherwise the rate of electrochemical kinetics is also significant. Steady-state methods

For steady-state experiments, the limiting current and half-wave potential (see below) are two convenient fixed points which allow comparison with theory. It is less common to analyse steady-state data by fitting entire voltammograms.

The species resulting from the homogeneous chemical reactions may be electroactive and thus augment the current measured at the electrode. By varying the rate of mass transport (controlling how much of the product species reach the electrode) and observing the variation in current augmentation, the rate of the chemical reaction may be measured relative to the rate of mass transport. The commonly used measure of current augmentation is Neff, the effective number of electrons transferred. It is defined as the ratio of the limiting current to that of a simple one electron transfer.

The coupling between a homogeneous process and a reversible electrode reaction can also be exploited to measure the homogeneous rate constant. The homogeneous reaction perturbs the equilibrium between the species in the electrochemical couple. The effect is to shift the voltammogram along the potential axis. The half-wave potential, E1/2, defined as the potential corresponding to half the limiting current, is used to measure this shift. Again the rate of mass transport may be varied to control the time which kinetically unstable species spend near the electrode, which modulates the shift in half-wave potential. This method suffers from a few more experimental complications than the measurement of the limiting current. The potential of the reference electrode must be constant throughout the experiment and have a known value. This may present problems with pseudo-reference electrodes such as silver or platinum44. Transient methods

In chronoamperometric experiments the current is recorded as a function of time by an oscilloscope, digital sampling etc. As the system relaxes, the rate of mass transport decreases (towards its steady-state value in the case of hydrodynamic or microelectrodes). The transient current is often analysed as the ratio to the steady-state value (I/Iss). There are no obvious choices for fixed points along the transient, so it is usual to fit theoretical and experimental curves.

In linear-sweep voltammetry, the peak currents and potentials provide natural 'fixed' points to compare with theory. In cyclic voltammetry, ratios of peak currents may be used to quantify homogeneous processes. The peak separation (along the potential axis) may be used to measure heterogeneous rate constants.

1.5.2 Choice of experimental technique: steady-state vs. transient experiments

The experimental technique and electrode geometry should be selected to match the kinetic time scale of the reaction being studied. For a given electrode geometry, a transient method is able to probe faster kinetics than a steady-state method since the rate of mass transport is higher at early times. In linear sweep or cyclic voltammetry, high scan rates can be used to accelerate the time scale beyond that of the steady-state response. Typically scan rates of 10 - 3000 Vs-1 are used, corresponding to a time scale of (1-10ms), although 'fast-scan' CV may go as fast as 1x106Vs-1, reaching time scales of 10ns45.

The price paid for this is precision. Steady-state methods give highly reproducible results. Often in transient methods, such as chronoamperometry or cyclic voltammetry, the current changes very rapidly with time - a small error in the measurement of the time at which the current is sampled could introduce a large error. This is particularly an issue when using transient methods with microelectrodes to probe very fast time scales - the response of the recording device (oscilloscope, A/D converter etc.) must be significantly faster than the experimental time scale. For this reason, steady-state methods are generally the most popular for quantitative work. Transient methods are also complicated by increased Ohmic drop (due to the larger currents passed) and capacitive currents after large jumps in, or rapidly scanning, the potential. Cyclic voltammetry is a generally popular method especially for qualitative work46. The main focus of this thesis will be on these methods: steady-state voltammetry at a microdisc, rotating disc, channel and wall-jet electrodes, representing increasing sophistication, in addition to microdisc cyclic voltammetry.

The models assume that Ohmic drop is negligible and that the current is purely Faradiac. If the Ohmic resistance can be measured experimentally, then modification of the steady-state models with an Ohm's law term to account for the potential drop would be relatively trivial. Cyclic voltammetry would be less straightforward, given the constantly changing potential with time46b.

1.6 The solution method

The electrode geometry, mechanism and experiment define the system of partial differential equations which needs to be solved in order to simulate the experiment. The final component of the model is therefore, not surprisingly, a method to solve such a system of equations. The difficulty of solving these systems depends on the complexity of the material balance equations, and they are linked to each other by the kinetic terms. The main considerations for the 'complexity' of the system are:

For the simpler cases where a low (usually 1) dimensional linear PDE may be solved in isolation, the system is analytically tractable. For anything more than model problems at most practical electrode geometries, numerical methods are currently the only way in which the equation systems may be solved.

1.6.1 Analytical solutions

For a simple electron-transfer it is possible to solve the time-dependent mass transport equation analytically15 for the transport-limited current at semi-infinite planar electrodes:


and spherical electrodes:

where (1.22)

For the microdisc, cylinder and band, the equation has not been fully solved analytically, but a number of approximate asymptotic solutions have been derived15,11. A detailed summary of those for the microdisc electrode is presented in Chapter 9.

Table 1.3: 'Levich' equations for the steady-state transport-limited current due to a simple electron transfer at hydrodynamic electrodes

For hydrodynamic electrodes, in order to solve the convective-diffusion equation analytically for the steady-state limiting current, it is necessary to use a first order approximation of the convection function(s) (such as the Lévêque approximation52 for the channel). These approximate expressions for the steady-state transport-limited currents are shown in Table 1.3.

For geometries such as the spherical electrode, where the mass-transport is a diffusion function in one dimension, it is possible to solve the material balance equations for a few first-order mechanisms (see Chapter 9). In order to tackle second order kinetics, more complex mechanisms and model other geometries with more complex mass transport, it is necessary either to resort to numerical methods, or make cruder approximations such as the reaction-layer approach, where the kinetically unstable species are assumed to be free from the effects of convection and have a linear concentration gradient. Integral transforms53 (e.g. Laplace, Hanckel etc.)

Integral transformations have been used to solve a range of electrochemical problems54-58. By transforming one (or more) of the variables in the material balance equations via an integral transformation, the problem may simplify to an ordinary differential equation which can be solved analytically. The 'difficult' part is solving the resulting integral equation (which is usually of the singular and nonlinear Volterra type56) and performing the inverse transform for the current . In some cases this may be done analytically, but otherwise the integral equation must either be numerically integrated or an approximation must be made, such as solving in a limit. The main disadvantages are that the transformation may only be applied to linear equations (ruling-out second-order kinetics) and for each system (geometry/mechanism/experiment) a new algebraic treatment is required, hence these methods are not at all suitable for 'general' simulators, unless provided in a library. In some cases there may also be difficulties with singularities when integrating the transformed co-ordinate. The main advantage of numerical integration vs. full numerical simulation is the reduced dimensionality of the transformed problem, requiring less memory and CPU time57.

1.6.2 Numerical simulations

These all rely on the approximation of a continuous quantity (usually a derivative) as an array of discrete values which may be manipulated using the simple (+,-,×,÷) operators available in a computer's instruction set. Finite difference methods

The space/time over which the problem is formulated is covered with a mesh of points (often referred to as 'nodes'). At each point the derivatives in the material balance equation are approximated as differences of the concentrations at the given and surrounding points.

Consider a spatial co-ordinate, x, divided into NGX points at a distance Dx apart. k is used to index a particular point (k=1...NGX). The concentration gradient at point k, may be represented in three ways:

; ; (1.23)

where uk is the concentration at point k. The first, known as upwind differencing, centres the derivative at point k-1/2. The last (downwind differencing) centres the derivative at k+1/2. Only the central difference spanning both points either side results in a derivative centred at point k. If either upwind or downwind differencing is used to represent the 1st derivative, the offset introduces a second order error known as numerical dispersion, which acts as an artificial viscosity or diffusion term.

The upwind and downwind differences may be combined to give an expression for the second derivative which is centred around point k:


The conversion of derivatives to differences leads to a set of linear equations (based on a 5-point stencil in two dimensions - each node is related to its 4 nearest neighbours) which can be solved to give the solution to the PDE. The methods are well-suited to simulations in rectangular regions, a situation which is often compatible with an electrochemical cell. These are by far the most popular methods for electrochemical simulations and are used throughout this thesis. Methods of weighted residuals (MWR)

Rather than approximating the differential operator as in the finite difference method, the differential equation is converted into a variational form (an integral equation) in which a stationary point (maximum or minimum) corresponds to the solution of the differential equation. A set of 'trial' (interpolation) functions is used to approximate the concentration distribution. The next two methods fall into this category:

Orthogonal collocation (OC) (Villadsen and Stewart59, 1967). This is an efficient technique based on polynomial curve fitting. Points are 'selected' along the co-ordinate where the polynomial fits the solution exactly. In order to minimise errors between the points (from approximation theory60) the points picked should be the zeros of an orthogonal set of polynomials. By choosing different orthogonal sets, it is possible to modify the distribution of points along the co-ordinate. Derivatives in the material balance equation are replaced with the derivatives of the polynomial. Using matrix algebra, the coefficients of the polynomials may be eliminated from the equations leaving a linear system of equations for the steady-state case, or a system of Ordinary Differential Equations (ODEs) for a time-dependent simulation. The high efficiency of the method more than offsets the computational effort required to integrate the ODE at each collocation point. See Britz61 or Speiser2 for more details. Yen and Chapman62 used a two-dimensional 'double collocation' method to simulate radial diffusion and migration at a rotating disc.

Finite element methods (FEM). The space is divided into a number of subregions known as elements, which join at points known as nodes. The integral equation then breaks up into a sum of integrals over each element. In an element the concentration distribution is approximated by the weighted sum of a set of trial (interpolating) functions of the nodal values - these functions are chosen to have a value of unity at just one node, and zero at all others. These functions are substituted into the variational form which leads to a set of linear equations for each element (the coefficient matrix is known as the 'stiffness matrix' and the RHS as the 'load vector'). Summing over all elements leads to a set of linear equations (which is generally symmetric) for the unknown (nodal) values.

One major advantage of finite element methods over finite difference methods is the way in which they naturally incorporate non-uniform meshes. They can therefore be applied to problems with a complex geometry - for example elevated and recessed electrodes63 and simulation of rough electrodes64. However, finite element methods are more complex to program, especially when simulating chemical steps, and result in a linear system of equations which is not neatly banded. Prentice and Tobias65 briefly discuss the relative merits of FE and FD for electrochemical simulations. Stephens and Moorhead66,67, Penczek et al.68,69, and Stephens and Fisher70-72 have simulated a number of diffusion-only and hydrodynamic problems using FEM. Cecilia et al.56 used a Galerkin FEM for the simulation of a 1-D diffusion-controlled adsorption process. Aoki et al.73 used a FEM simulation of the charge transfer in a polymer film coating the electrode. Finite analytic method (Chen74, 1984)

The space is divided up into a mesh of (rectangular) elements. In each element it is assumed that the coefficients of the derivatives in the PDE may be treated as constants. If suitable boundary conditions are specified, the PDE can be solved in the local element using a separation of variables. The value in the particular element is related to that of its neighbours though the boundary conditions, leading to a set of linear equations based on a nine-point stencil in two dimensions (i.e. each element is related to its 4 nearest and 4 next-nearest neighbours). For more details on this method see Jin et al.75,76 Network approach

It is possible to represent the entire electrochemical system including the instrumentation (potentiostat etc.) as a single electrical circuit. The solution is usually spatially discretised into a network of resistance/capacitance elements77,78 and then solved using a circuit simulation program such as SPICE or PSPICE79.

1.7 Outline of the work presented in this thesis

Figure 1.9: Components of an electrochemical finite difference simulation

Figure 1.9 shows how the components of the electrochemical simulation fit together. It is clear that building an electrochemical simulation model, and especially solving the resulting equation system, can be a difficult task. For complex mechanisms at microelectrodes, simultaneous, often non-linear, multidimensional partial differential equations must be solved.

In the next chapter, the numerical 'nuts and bolts' (triangular symbols in Figure 1.9) are examined: finite difference schemes for discretisation of material balance equations and the algorithms used to solve the resulting sets of simultaneous equations. These solvers can consume an enormous amount of resources, both in memory and especially in CPU time, because in some cases over a million points are necessary to describe the concentration distribution accurately80. In Chapter 3, ways in which the efficiency of the finite difference simulations may be improved are investigated.

Two state-of-the-art solver methods are introduced in Chapter 4. The multigrid method, arguably the most efficient method for solving linear systems of finite difference equations, is applied to electrochemical simulations. Preconditioned Krylov subspace methods are introduced which allow the finite difference representation of an electrochemical experiment to be treated as a general sparse linear system. This approach enables complex reaction mechanisms to be simulated without resorting to an explicit approximation. It also allows the simulation of multidimensional mass transport in conformal space and would, in principle, allow simulation of the potential distribution leading to migration effects. The merits of both methods are realised in a new linear solver - a multigrid method based on PKS methods. Finally, the frontal solution of pseudo-2D problems such as the ChE and WJE is addressed. The PKS methods provide an alternative to Rudolph's FIFD method for fully implicit simulation of the coupling kinetics in these systems. A parallel algorithm for time-dependent frontal simulation is also developed.

One of the criteria for the selection of the linear solvers introduced in this work is that they are available as subroutines in common mathematical libraries. This enables electrochemists to use them without having to learn how to write, debug and optimise them. However, the experimental electrochemist is still faced with an ominous task if a model must be custom-built for every experiment. There is clearly a need for a general electrochemical simulator which can handle a range of mechanisms, geometries and experiments and therefore can be applied to a wide range of experimental systems. The FIFD method achieves this to a limited extent in the commercial package, DigiSim(tm). Although this package is restricted to simulating cyclic voltammetry and one-dimensional diffusion, it does allow a general electrochemical mechanism to be simulated. In Chapter 5 a new method is presented which can simulate a general mechanism at any geometry for a range of electrochemical experiments. Speiser2 raised concern about the danger of experimentalists misapplying a 'black box' which they do not fully understand. For example DigiSim(tm) cannot simulate two-dimensional diffusion to a microdisc electrode and experimentalists may attempt to approximate the microdisc behaviour using a planar or hemispherical model which could lead to inaccurate results81 (see Chapter 9 for more about this). The answer to this particular case is simple - a well-written piece of simulation software should include an assessment of its limitations and check the input parameters make physical 'sense' (for example, attempting to simulate a 1mm 'planar' disc electrode should at least provoke a warning!). It is counterproductive to force each experimentalist to acquire detailed knowledge of the numerical methods used to solve the PDEs. However there are parameters such as the number of nodes, which need to be adjusted to ensure the results obtained are accurate. These must be drawn to the attention of the user in the operation and documentation of the software so that they are set correctly. Ideally software should 'check' such parameters are set correctly (e.g. it would be possible to perform the simulation on the user-defined mesh and one slightly coarser/finer and note the difference).

Chapters 6-9 apply the simulation methods developed in this thesis to investigate the behaviour of microdisc and hydrodynamic electrodes: channel (microband), wall-jet (or tube) and rotating disc electrodes. Particular attention is devoted to the assessment of approximate analytical equations which are much more easily applied than numerical simulations. This includes the approximation of the response of a microdisc electrode with that of an 'equivalent' spherical electrode.

The electrochemical simulation can be viewed as a 'black-box' which, given a set of experimental parameters and a proposed mechanism, predicts a value of the experimental observable (e.g. the current). Electrochemical experiments usually involve recording a 'data set' comprising of a number of such 'points' across a range of time scales (by varying potential, flow rate etc.). At the very least, the simulation must be run once for each of these points. Often the object of the experiment is to measure a rate constant and/or discriminate between mechanisms by minimising the error between experimental and simulated data sets, so the simulation must be run many times for each current in each data set. As later chapters will illustrate, even using very sophisticated linear solvers, the amount of CPU time required to simulate each point can make the entire procedure laborious or prohibitively expensive. There is an alternative approach: simulated results can be stored for a library of mechanisms in the form of working curves or surfaces from which the observable can be interpolated for the experimental conditions. In Chapters 6-9, working curves and surfaces are generated for common mechanisms occurring at microdisc, rotating disc, wall-jet and channel (microband) electrodes. In Chapter 10, working surface interpolation is investigated and a service is established to allow fitting of experimental data using these curves and surfaces via the World Wide Web. The thesis concludes with a quantitative comparison of the kinetic time scales which may be investigated using steady-state voltammetry at the electrode geometries discussed in this work.


1 R.W. Cattrall, Chemical Sensors, Oxford Chemistry Primers, OUP, (1997), p30.
2 B. Speiser, Electroanalytical Chemistry, ed A.J. Bard, Vol 19, Marcel-Dekker, NewYork, (1997).
3 A.C. Testa, W. Reinmuth; Anal. Chem., 33, (1961), 1320.
4 Q. Zhuang, H. Chen, J. Electroanal. Chem., 346, (1993), 29.
5 L.K. Bieniasz, B. Speiser, J. Electroanal. Chem., 441, (1998), 271.
6 R. Meunierprest, E. Laviron, J. Electroanal. Chem., 410, (1996), 133.
7 J.A.V. Butler, Trans. Faraday Soc., 19, (1924), 729 and 734.
8 T. Erdey-Gruz and M.Volmer, Z. Physik. Chem., 150A, (1930), 203.
8b P.H. Reiger, Electrochemistry, 2nd Edition, Chapman and Hall, (1994), p319.
9 F. Prieto, R.D. Webster, J.A. Alden, W.J. Aixill, G.A. Waller, R.G. Compton, J. Electroanal. Chem., 437, (1997), 183.
10 E. Laviron, J. Electroanal. Chem., 391, (1995), 187.
11 C.A. Amatore, B. Fosset, Anal. Chem. 68, (1996), 4377.
12 A.J. Bard, L.R. Faulkner, Electrochemical Methods: Fundamentals and Applications, Wiley, New York, (1980), p145.
13 J. Heyrovsky, Chem Listy, 16, (1922), 256.
14 M. Fleischmann, S. Pons; Ultramicroelectrodes, Datatech, Morganton, NC (1987); J. Wang; Microelectrodes, VCH, New York (1990); I.N. Montenegro; Research in Chemical Kinetics, pg. 1, Elsevier, London, (1994); J. Robinson, Comp. Chem. Kinetics, 29, (1989), 150.
15 R.M. Wightman, D.O Wipf, Electroanalytical Chemistry, ed A.J. Bard, Vol 19, (1987), 267.
16 K. Aoki, Electroanalysis, 5, (1993), 627.
17 Microglass instruments, Greensborough, Victoria, Australia. (email chemk@lure.latrobe.edu.au).
18 W.G. Oldham, Scientific American, Script 77, 110.
19 See references cited in [15], p328.
20 J.A. Alden, J. Booth, R.G. Compton, R.A.W. Dryfe, G.H.W. Sanders, J. Electroanal. Chem., 389, (1995), 45.
21 R.M. Wightman, E. Strope, P.M. Plotsky, R.N. Adams, Nature, 262, (1976), 145; J.B. Zimmerman and R.M. Wightman, Anal.Chem., 63, (1991), 24.
22 J.B. Justice (ed.), Voltammetry in the Neurosciences, Humana, Clifton, New Jersey, (1987).
23 R. Menzel, Spectral sensivity and colour vision in invertebrates (a review) in Handbook of physiology vol. VII-6A, Ed H. Autrum, Springer-Verlag, Berlin, (1979).
24 Y. Saito, Rev. Polarogr. 15, (1968), 177.
24b C.A. Amatore in I Rubinstein, Physical Electrochemistry, (1995), p131.
25 W.J. Albery and M.L. Hitchman, Ring Disc Electrodes, Oxford University Press, Oxford, (1976).
26 J. Wang, J. Electroanal. Chem., 140, (1982), 141.
27 E. Levart, J. Electroanal. Chem., 187, (1985), 247.
28 J. Yamada, H. Matsuda, J. Electroanal. Chem., 44, (1973), 189.
29 J.V. Macpherson, S. Marcar and P.R. Unwin, Anal. Chem., 66, (1994), 2175; J.V. Macpherson, S. Marcar and P.R. Unwin, J. Chem. Soc. Faraday Trans., 91, (1995), 899.
30 R.G Compton, P.R. Unwin, J. Electroanal. Chem., 205, (1986), 1.
31 N. Ibl and U. Dossenbach, Comprehensive Treatise of Electrochemistry, Plenum, New York, 6, (1983), 133.
32 R.G. Compton and B.A. Coles, J. Electroanal. Chem., 144, (1983), 87; W.J. Albery and S. Bruckenstein, J. Electroanal. Chem., 144, (1983), 105.
33 R.E. Meyer, M.C. Banta, P.M. Lantz, F.A. Posey, J. Electroanal. Chem., 30, (1971), 345.
34 J.A .Alden, M.A. Feldman, E. Hill, F. Prieto, M. Oyama, B.A. Coles and R.G. Compton, Anal. Chem., 70, (1998), 1707.
35 H.R. Corti, D.L. Goldfarb, M.S. Ortiz, J.F. Magallanes, Electroanalysis, 7, (1995), 569.
36 N.V Rees, R.A.W Dryfe, J.A Cooper, B.A Coles, R.G Compton, S.G Davies, T.D McCarthy, J. Phys. Chem., 99, (1995), 7096.
37 N.V. Rees, J.A. Alden, R.A.W. Dryfe, R.G. Compton, B.A. Coles, J. Phys. Chem., 99, (1995), 14813.
38 R.A.W. Dryfe, D. Phil Thesis, University of Oxford, (1996), p248.
39 J.A. Cooper and R.G. Compton, Electroanalysis, 10, (1998), 141.
40 P.R. Unwin, R.G. Compton, Comprehensive Chemical Kinetics, 26, (1989), 173.
41 R.G. Compton, A.C. Fisher, G.P. Tyley, J. Appl. Electrochem., 21, (1991), 295.
42 R.G. Compton, P.R. Unwin, J. Electroanal. Chem., 205, (1986), 1.
43 C.G. Phillips, J. Electroanal. Chem., 296, (1990), 255.
44 D.T. Sawyer, A. Sobkowiiak, J.L. Roberts, Electrochemistry for Chemists, Wiley, New York, (1995), p197.
45 D.O. Wipf, R.M. Wightman, Anal. Chem. 60, (1988), 2460.
46 F.A. Armstrong, H.A.O. Hill, N.J. Walton, Acc. Chem. Res., 21, (1988), 407.
46b J. NavarroLaboulais, J Trijueque, J.J. GarciaJareno, F. Vicente, J. Electroanal. Chem., 422, (1997), 91.
47 V.G. Levich, Physicochemical Hydrodynamics, Prentice-Hall, Engelwood Cliffs, New Jersey, (1962).
48 W.J. Albery, S. Bruckenstein, J. Electroanal. Chem., 144, (1983), 105.
49 J.V. Macpherson, S.Marcar, P.R. Unwin, Anal. Chem., 66, (1994), 2175.
50 G. Wranglen, O. Nilsson, Electrochimica Acta, 7, (1962), 121.
51 W.J. Blaedel, C.L. Olsen, L.R. Sharma, Anal. Chem., 35, (1963), 2100.
52 M.A. Lévêque, Ann. Mines. Mem. Ser., 12/13, (1928), 201.
53 J.W. Miles, Integral Transforms in Applied Mathematics, Cambridge, (1971).
54 J.C. Myland, K.B. Oldham, Fundamentals of Electrochemical Science, Academic Press, San Diago, California, (1994), p21.
55 S. Coen, D.K Cope, D.E. Tallman, J. Electroanal. Chem.,215, (1986), 29; D.K. Cope, C.H. Scott, U. Kalapathy, D.E. Tallman, J. Electroanal. Chem., 280, (1990), 27; D.K. Cope and D.E. Tallman, J. Electroanal. Chem., 396, (1995), 265.
56 J. Cecilia, J. Galceran, J. Salvador, A.J. Puy, F. Mas, Int. J. Quant. Chem., 51, (1994), 357; J. Galceran, J. Salvador, J. Puy, J. Cecilia, D.J. Gavaghan, J. Electroanal. Chem., 440, (1997), 1.
57 M.V. Mirkin, A.J. Bard, J. Electroanal. Chem., 323, (1992), 1.
58 W.M. Leslie, J.A. Alden, R.G. Compton, T.Silk, J. Phys. Chem, 100, (1996), 14130.
59 J.V. Villadsen, W.E. Stewart, Chem. Eng. Sci., 22, (1967), 1483.
60 L. Fox, An introduction to Numerical Linear Algebra, OUP, (1964).
61 D. Britz, Digital Simulation in Electrochemistry, 2nd Edition, Springer-Verlag, (1988).
62 S.C. Yen, T.W. Chapman, J. Electrochem. Soc., 143, (1987), 1964.
63 R. Ferrigno, P.F. Brevet, H.H. Girault, Electrochimica Acta., 42, (1997), 1895.
64 J. Fleig, J. Maier, Solid State Ionics, 94, (1997), 199.
65 G.A. Prentice and C.W. Tobias, J. Electrochem. Soc., 129, (1982), 72.
66 M.M Stephens, E.D. Moorhead, J. Electroanal. Chem., 220, (1987), 1.
67 E.D. Moorhead, M.M. Stephens, J. Electroanal. Chem., 282, (1990), 1.
68 M. Penczek, Z. Stojek, J. Osteryoung, J. Electroanal. Chem., 170, (1984), 99.
69 M. Penczek, Z. Stojek, J. Electroanal. Chem., 181, (1985), 83.
70 N.P.C. Stevens, S.J. Hickey, A.C. Fisher, Anales de Quimica, 93, (1997), 225.
71 N.P.C. Stevens, A.C. Fisher, J. Phys. Chem. B, 101, (1997), 41.
72 N.P.C. Stevens, A.C. Fisher, Electroanalysis, 10, (1998), 16.
73 K. Aoki, K. Tokuda, H. Matsuda, M. Oyama, J. Electroanal. Chem., 176, (1983), 41.
74 C.J. Chen, J. Comput. Phys., 53, (1984), 209.
75 B.K. Jin, W.J. Qian, Z.X. Zhang, H.S. Shi, J. Electroanal. Chem., 411, (1996), 19.
76 B.K. Jin, W.J. Qian, Z.X. Zhang, H.S. Shi, J. Electroanal. Chem., 417, (1996), 45.
77 B.A. Coles, R.G. Compton, R.A. Spackman, Electroanalysis, 5, (1993), 41.
78 J. Horno, M.T. García-Hernández, C.F. González-Fernández, J. Electroanal. Chem., 352, (1993), 83.
79 For more information see http://www.eeap.aston.ac.uk/eeap/documents/user-guides/pspice/pspice-ug.html.
80 R.G. Compton, M.B.G. Pilkington, G.M. Stearn, J.Chem.Soc. Faraday Trans. 1, 84, (1988), 2155-2171.
81 J.A. Alden, F. Hutchinson, R.G. Compton, J. Phys. Chem. B, 101, (1997), 949.