Electrochemical tomography as a nondestructive technique to study localized corrosion of metals

Abstract

We present an approach, termed electrochemical tomography (ECT), for the in-situ study of corrosion phenomena in general, and for the quantification of the instantaneous rate of localized corrosion in particular. Traditional electrochemical techniques have limited accuracy in determining the corrosion rate when applied to localized corrosion, especially for metals embedded in opaque, porous media. One major limitation is the generally unknown anodic surface area. ECT overcomes these limitations by combining a numerical forward model, describing the electrical potential field in the porous medium, with electrochemical measurements taken at the surface, and using a stochastic inverse method to determine the corrosion rate, and the location and size of the anodic site. Additionally, ECT yields insight into parameters such as the exchange current densities, and it enables the quantification of the uncertainty of the obtained solution. We illustrate the application of ECT for the example of localized corrosion of steel in concrete.

Introduction

Localized corrosion is a common phenomenon in a wide range of materials and environments, including corrosion of underground metallic structures (e.g., pipes, storage tanks, borehole casings, geotechnical, etc.)^1,2,3,4, reinforced and prestressed steel in concrete^5,6,7,8,9, and medical implants such as vascular stents in the human body¹⁰. Since corrosion potentially impairs the functionality and serviceability of these different metal objects, and as their potential failure presents serious risks to human safety, health, and the environment, there is a need for methods, that quantify in-situ the local corrosion rate. For the non-destructive measurement of instantaneous corrosion rates on metal objects in their actual service environments (as opposed to laboratory samples and conditions), generally, only electrochemical techniques are feasible^{11,12,13,14,15,16}.

A particular problem relevant to all examples of corrosion in engineering mentioned above, is that the metal is surrounded by a solid, opaque, porous medium. This complicates the non-invasive measurement of the corrosion rate for various reasons. First, any measurement probes (such as reference and counter electrodes in electrochemical techniques) have to be positioned at a certain minimum distance from the metal. This generally dilutes the acquired signals, and thus limits the accuracy and reliability of the measurement. Second, since the corroding metal cannot be visually inspected, the anodically behaving area of the galvanic corrosion cell remains unknown. Therefore, it is difficult to relate the measured corrosion current to an area, and subsequently convert it to a corrosion current density, or corrosion rate in terms of the thickness or mass loss per unit time. These aspects have received considerable attention in the field of chloride-induced corrosion of steel in concrete since the 1980s^{13,17,18,19,20,21}. A number of approaches have been proposed to overcome the problem of the unknown anodic surface area, which is recognized to potentially induce an error in the obtained corrosion rate by up to several orders of magnitude. Such approaches include attempts to confine the polarization current applied during electrochemical corrosion rate measurements^16,17,22,23. However, the feasibility and success of these approaches have been repeatedly questioned, both on the grounds of experimental validation studies^16,22 and theoretical reasoning^24,25,26,27.

The unknown anodic surface area in localized corrosion of metals embedded in porous media remains the major limitation for non-destructively quantifying local corrosion rates. This fundamental problem is relevant for virtually all different electrochemical techniques. It universally limits the interpretation of results from techniques based on linear polarization resistance measurements, galvanostatic pulse measurements, electrochemical impedance spectroscopy, etc.

Here, we present an approach, for the in situ study of the instantaneous state of corroding systems (see Fig. 1a), termed electrochemical tomography (ECT). The approach is based on combining electrochemical measurements with a numerical model, which then serves in the solution of a Bayesian inverse problem. The electrochemical measurements are performed non-invasively, at the surface of a porous electrolytic domain, in which the locally corroding metal is embedded. We consider it a particular advantage that, while quantifying the instantaneous rate of localized corrosion, the approach not only yields the corrosion current, but also the location and size of the anodic area. Hence, ECT enables the determination of the actual corrosion rate, in terms of thickness or mass loss per unit time, without the need for assumptions related to the anode area. A further advantage of this approach is that it is capable of yielding insight into the fundamental electrochemistry of the system under study, namely through determining parameters such as the anodic and cathodic exchange current densities and Tafel slopes. Additionally, the use of a stochastic inverse method allows us to take into account the inherent uncertainty within the various model parameters.

Fig. 1: The methodology of electrochemical tomography (ECT) and the importance of external polarization.

a Overview of ECT. b The geometry of the 2D axisymmetric forward model, including model parameters W_an and z_an describing the size and location of the anode and the position of the reference electrodes (RE) and counter electrodes (CE) at the surface of the concrete. The marginal probability density function (pdf) estimates of the MCMC results for simulations with and without external polarization of −0.2 and +0.08 mA and a spacing of Δz = 5 cm between the RE, for model parameters c W_an, d z_an, e log10(i_0,Fe), f ${\rm{log}}10({i}_{0,{{\rm{O}}}_{2}})$. The vertical dashed lines represent the means of the corresponding MCMC results, the black lines the known (input) model parameters.

Inverse methods are widely applied in engineering and form the basis of many non-destructive testing methods^28,29. In the field of localized corrosion, however, the inverse approaches are still in their infancy. Several authors performed pioneering work and showed the potential of these methods for localized corrosion in reinforced concrete^30,31,32,33. The most recent study by Marinier and Isgor³³ was the first to use a numerical model based on the Butler−Volmer kinetics, describing the polarization behaviour of steel. They applied a deterministic inverse method on potentials measured under open-circuit conditions. They demonstrated that this approach is able to localize the anode. However, they encountered problems in the stability of their method, when a smaller, but more practical number of measurements were used. This is a problem that our proposed approach overcomes; firstly by including the response of the electrical field to external polarization and secondly, by applying a stochastic inverse method, instead of using regularization methods, that are commonly used in deterministic approaches to overcome instability problems.

In the current work, we illustrate the application of ECT for corrosion rate measurements by applying the method to the example of localized corrosion of steel in concrete. A first reason for selecting the example of steel corrosion in concrete, is that this particular case of localized corrosion has so far received considerable attention in the literature on corrosion rate measurements. An additional reason is the urgent need for techniques to quantify corrosion processes in situ in civil engineering structures. Various studies have revealed the staggering socio-economic relevance and huge economic impact of corroding concrete structures^34,35, and the potential risks for human safety are well known. The need for reliable inspection methods will only grow, because of the ageing of concrete structures. In the next 20 years, for example, the fraction of bridges with an age greater than 50 years will rise from currently 40−50% to 70−80% in all industrialized countries³⁶. The consequence will be a wave in the need for the replacement of this vital infrastructure. As a result, reliably diagnosing the condition of structures, in order to prioritize their repair or maintenance, will become increasingly important.

Results

ECT applied for the example of localized corrosion in reinforced concrete

A typical case of chloride-induced reinforcing steel corrosion in concrete was considered. A locally corroding steel bar segment, embedded in a concrete cylinder, was assumed (see Fig. 1b). For the electrochemical measurements, an array of reference electrodes (RE), as well as two counter electrodes (CE), were positioned on the concrete surface. The counter electrodes enable a temporary polarization of the corroding system, and the reference electrodes are used to measure the resulting system response.

ECT—described in detail in the ‘Methods’ section below—seeks to find the location and the size of the anodic part of the galvanic cell by means of a deterministic forward model and a stochastic inverse solver. The forward model gives the electrical potential field in a certain conducting environment (the electrolyte) as a function of the location, size, and electrochemical properties of an actively corroding metal surface, and under the effect of external polarization. The inverse solver seeks to invert for certain parameters encoding the state of the corroding steel, using the electrical potential field measured at the surface of the concrete (${\overrightarrow{\phi }}^{{\rm{obs}}}$), including the potential fields that were obtained after polarizing the corroding system by help of the counter electrodes. In this example, the parameters used in the forward model to simulate these system responses at the electrolyte surface, hereafter referred to as ’model parameters’, include, besides the size (W_an) and location (z_an) of the actively corroding metal surface, also the exchange current densities for iron oxidation (i_0,Fe) and oxygen reduction (${i}_{0,{{\rm{O}}}_{2}}$) (see ‘Methods’ section for more details).

The Bayesian inverse approach applied here, as opposed to deterministic methods, does not only supply us with the ’best values’ for the sought after model parameters, but also gives us the posterior probability: the probability distribution of the model parameters given the observed electrical potentials (${\overrightarrow{\phi }}^{{\rm{obs}}}$) and a priori knowledge about the model parameters (the prior probability densities). The marginal posterior distributions of the model parameters are subsequently sampled using a Markov chain Monte Carlo (MCMC) algorithm.

To demonstrate and optimize the approach for the example of reinforced concrete, synthetic electrochemical measurements (${\overrightarrow{\phi }}^{{\rm{obs}}}$) were generated using a set of known model parameters. These synthetic data were then used as input for the inverse model (see ‘Methods’ for more details).

Increasing the size of the measured data vector

Figure 1c−f shows how ECT succeeds in finding the location and size of the anode as well as the exchange current densities. These specific results also show that the temporary external polarization of the system is important in enhancing the inverse reconstruction of the corroding system. Simulations show that without external polarization, one has to go to an impractically small stepsize (<5 cm) between the measurements at the concrete surface, in order to see improvement in the MCMC results (see Supplementary Figs. 2 and 3). This finding is in general agreement with the earlier study of Marinier and Isgor³³, where the authors had to use a large amount of closely spaced data points to reach stable results in their inversion. The results of the current work show that external polarization is much more effective in reducing the ill-posedness of the problem. Decreasing the stepsize from 5 to 1.25 cm, to increase the size of the measured data vector, has considerably less effect on the resulting MCMC results, than tripling the size of the data vector by externally polarizing the steel with just one negative and one positive current. Using external polarization, a spacing of 5 cm seems to already result in relatively accurate estimations of the model parameters.

Information gain

To assess how successful the inversion was in finding the model parameters, we compute the information gain. The information gain, in our case, is determined to be the difference between the sampled posterior (the MCMC inversion results) and the prior probability densities, which describe the a priori belief about the model parameters. The information gain tells us how much our knowledge about these parameters has improved, relative to the assumed prior information. The information gain is computed using the Kullback−Leibler divergence ³⁷. The Kullback−Leibler divergence is a special case of Shannon’s measure of information content, and is given for discrete probability distributions by

$${D}_{{\rm{KL}}}(P| | {P}^{{\rm{ref}}})=\mathop{\sum}\limits_{x\in {\mathbb{M}}}P(x){\rm{log}}10\left(\frac{P(x)}{{P}^{{\rm{ref}}}(x)}\right)$$

(1)

where P is the sampled posterior, P^ref the prior probability density, ${\mathbb{M}}$ the model space, and D_KL the Kullback−Leibler divergence in units of dits.

Figure 2 shows the information gain for varying quantities M, that is the number of considered cases of different external currents I_CE applied at the counter electrodes, varying from M = 3 to M = 11. As the improved results shown in Fig. 1c−f were obtained after the application of just one positive and one negative polarization current (M = 3), one could expect the MCMC results to improve further if a larger set of discrete currents were to be applied. However, Fig. 2 shows only a gradual increase in the information gain. This improvement is relatively small, considering that an increasing number of I_CE, leads to a substantial increase of required computational power to perform the inversion. The computation time of simulations with M = 11 is three times greater, than that with M = 3. It seems that the extreme values of the applied range of I_CE weigh the heaviest, and applying additional currents within this range has only a limited effect. These extreme values of −0.2 mA and +0.08 mA are the maximum external currents, for which our forward model is still valid. Taking into account the computation time, we use M = 9, instead of performing the inversion for larger lengths of I_CE.

Fig. 2: Information gain upon increasing the number of different cases of external polarization.

The information gain (Kullback−Leibler divergence) D_KL, for the four model parameters, for an increasing number of different cases of external polarization (M). For M = 3, the applied current I_CE = [−0.2, 0, 0.08] mA. For M = 7, I_CE = [−0.2, −0.08, −0.02, 0, 0.02, 0.04, 0.08] mA. For M = 9, I_CE = [−0.2, −0.08, −0.04, −0.02, 0, 0.02, 0.03, 0.04, 0.08] mA. And finally for M = 11, I_CE = [−0.2, −0.08, −0.04, −0.02, −0.008, 0, 0.008, 0.02, 0.03, 0.04, 0.08] mA.

Figure 3 shows the probability density function (pdf) estimate of the MCMC results for M = 9, together with the assumed prior probability densities. These plots allow us to visualize the information gain for each model parameter as given in Fig. 2. The inverse solver is able to recover all model parameters and for most parameters the information gain is substantial. Especially concerning the width of the anode, W_an, our a priori knowledge was very limited, as it was assumed that the anode width is equally likely anywhere between 0 and 22 cm (a flat pdf). The results show a large information gain and allow us to confidently estimate this width between 4 and 10 cm. Concerning the location of the centre of the anode, z_an, the MCMC results show a narrow posterior. The location can be estimated within an error range of just a couple of cm.

Fig. 3: Markov chain Monte Carlo (MCMC) sampling results for simulations with M = 9.

The blue histograms show the probability density (pdf) function estimates of the MCMC results. The red curves are the corresponding prior probability densities for a W_an, b z_an, c log10(i_0,Fe), d log 10(i₀,_O2). The vertical lines represent the known ground truth value, the dashed lines the mean of the MCMC results.

We are also able to accurately estimate the exchange current densities log10(i_0,Fe) and ${\rm{log}}10({i}_{0,{{\rm{O}}}_{2}})$. The information gain for log10(i_0,Fe) is smaller compared to the other parameters (see also Fig. 2). An additional simulation, where a shifted prior probability density function was applied, shows, however, that the posterior seems to be mostly independent of the assumed prior (Fig. 4). This gives us confidence that the MCMC did not simply sample the prior distribution for log10(i_0,Fe), but reflects the posterior probability density. Moreover, the results in both Figs. 3 and 4 show us that even if the underlying assumptions (prior knowledge) are very imprecise, the inverse solver is still able to obtain good results.

Fig. 4: Successful solving even with imprecise prior knowledge.

MCMC results for simulations with M = 9 and a shifted prior probability distribution for log10(i_0,Fe). The blue histograms show the probability density function (pdf) estimate of the MCMC results for parameters a W_an, b z_an, c log10(i_0,Fe), d ${\rm{log}}10({i}_{0,{{\rm{O}}}_{2}})$ with the corresponding prior probability densities (in red). The vertical lines represent the known ground truth values, the dashed lines the mean of the MCMC results.

Convergence diagnostics

Finally, we want to stress the necessity to assess convergence when applying MCMC sampling to estimate the posterior probability density. One has to evaluate if sufficient samples were collected to create a histogram that well approximates the posterior probability density. Especially in this example where we are inverting for four different parameters, it is important to evaluate if we have achieved sufficient mixing using the number of samples that were collected. As an a priori analytical convergence rate is impossible to compute, we perform a statistical analysis on the output sample data, to test if convergence has been reached.

We apply a method originally described by Gelman et al.³⁸ and later generalized by Brooks and Gelman³⁹. It compares Markov-chains with different starting positions. Convergence is assumed to be reached when the variability of the individual chains is comparable with the variability of the pooled samples. Figure 5a presents part of the convergence diagnostics set by Brooks and Gelman³⁹. It shows the evolution of the estimation of the ’scale reduction factor’, $\hat{R}$, during the sampling. The ’scale reduction factor’ is the ratio between the variance estimate of the target distribution and the average of the variances of the different Markov chains. The chains are likely to have converged, when $\hat{R}$ has approached 1. The complete convergence diagnostics and details of the computation of $\hat{R}$ can be found in the supplementary note.

Fig. 5: Convergence assessment of the Markov-chains.

a The evolution of $\hat{R}(k)$ over the samples for the four different model parameters. b The evolution of the information gain (Kullback−Leibler divergence) D_KL over the samples, normalized versus the divergence at N = 60,000 samples. In around 30,000 samples $\hat{R}(k)$ has approached 1 and the information gain is stable, suggesting the Markov chains have converged and the posterior is sufficiently sampled.

Figure 5a shows that in fewer than 20,000 iterations, $\hat{R}(k)$ is very close to 1 for all model parameters. This means that we are then likely sampling from the target stationary distribution: the posterior. To evaluate when we have sufficiently sampled this posterior, we compute the information gain as a function of the amount of samples taken (see Fig. 5b). At 30,000 samples the information gain is stable for all model parameters, meaning that further sampling does not change the information content. Therefore we conclude that we have reached convergence and sufficiently sampled the posterior in around 30,000 samples.

Discussion

The described approach, electrochemical tomography (ECT), is a promising method for the instantaneous determination of local corrosion rates of metals embedded in porous media. As opposed to traditional electrochemical methods, such as linear polarization resistance (LPR), its accuracy is not limited by the lack of information about the surface area of the anodic part of the corroding steel.

To illustrate the accuracy and precision with which ECT is able to estimate corrosion rates, Fig. 6b shows the normalized distribution of the corrosion current densities, obtained for every set of model parameters, which were sampled during the inversion. These corrosion current densities were computed by averaging over the anodic steel area.

Fig. 6: Comparison of ECT to the traditional LPR method.

a Overview of the setup of the modelled LPR method, with depicted distances in [cm]. Counter electrodes (CE) are placed adjacent to the reference electrode (RE). Corrosion rates are obtained for different positions of the RE-CE assembly over the concrete surface. b The blue histogram shows the probability density function estimate of the sampled corrosion current densities, corresponding to every set of model parameters sampled during the MCMC. The black line represents the corrosion rate corresponding to the known ground truth model parameters, the dashed black line the mean of the sampled corrosion rates, as obtained from ECT. The red curves represent the corrosion rates obtained by the modelled LPR technique, as a function of position of the RE-CE assembly and for different assumed corroding steel areas [% of the total steel area], with and without IR compensation.

For comparison, in the same figure, corrosion rates, which can be detected with the traditional LPR method, are indicated. These corrosion rates were obtained with a numerical model similar to the forward model used for ECT, described in the ‘Methods’ section below. In the model, an LPR measurement was simulated by applying a polarization current of +/− 5 μA at two counter electrode rings with a width of 3.5 cm, positioned directly left and right of one reference electrode (RE), and recording the related change in potential at this RE (see Fig. 6a). This change in potential is used to compute the polarization resistance, R_p which is defined as the ratio between the measured voltage and the applied current, $\frac{{{\Delta }}\phi }{{{\Delta }}I}$. Subsequently, the corrosion current can be computed:

$${I}_{{\rm{corr}}}=\frac{1}{{R}_{p}}* B=\frac{1}{{R}_{p}}* \frac{{\beta }_{{\rm{an}}}{\beta }_{{\rm{cath}}}}{2.303* ({\beta }_{{\rm{an}}}+{\beta }_{{\rm{cath}}})}$$

(2)

where β_an and β_cath are the anodic and cathodic Tafel slopes, here considered to be the set values for the anodic Tafel slope of iron and the cathodic Tafel slope of oxygen as defined for the ECT application (see the section on ‘The example of localized corrosion in reinforced concrete’ under ‘Methods’).

The LPR measurements were simulated for different positions of the RE-CE assembly at the concrete surface (see Fig. 6a). To account for ohmic resistance of the concrete, the LPR measurements were corrected for the IR drop. To this aim, the IR drop was obtained by modelling the primary current distribution in the actual three electrode setup, namely for the different RE-CE positions with respect to the steel bar (working electrode). Finally, for the conversion of corrosion current to corrosion current density, different assumptions for the anode area were made in agreement with literature. For the application of LPR on reinforced concrete, an anode area of 10% of the total steel area is commonly suggested as a conservative limit⁴⁰. Note that the ’true’ anode surface area, corresponding to the ground truth of W_an, is 7.5%.

Figure 6b shows that it is possible to estimate an accurate corrosion rate using both the traditional LPR method and ECT. For LPR, however, this requires that the RE-CE assembly is positioned directly over the anode (zero-offset), that the IR drop is compensated, and that the correct assumption for the anode surface area is made (between 5 and 10%). Figure 6b shows that the LPR method is heavily sensitive to both the distance to the anodic surface, as well as the assumed anodic surface area. Using a recommended surface area of 10%, here results in an underestimation of the actual corrosion rate, even if the RE-CE assembly is positioned directly above the anode. The corrosion rates sampled by ECT on the other side, are neither dependent on the setup of the measurements, nor on estimations of the anodic steel area. Finally, not only is ECT able to more accurately estimate the corrosion rate compared to the LPR technique, it also gives insight in the precision of this estimation, which we consider a major advantage.

Besides the accurate estimation of the corrosion current density for localized corrosion, a major advantage of ECT is that it also supplies us with additional information. First of all, it is able to determine the location and the size of the anodic surface. Not only is this information crucial for the correct interpretation of other electrochemical techniques, it could also be applied to monitor the spread of the anodic surface over time, giving us more insight in the corrosion behaviour of steel in concrete.

Secondly, we are able to also determine more fundamental parameters, describing the polarization behaviour of the steel. The demonstration of ECT for the example of reinforced concrete has shown us that we are able to accurately invert for the exchange current densities of iron oxidation and oxygen reduction. These parameters have theoretically constant values for a certain metal in a set environment. However, especially for the case of steel in porous media, these values are not accessible from in-situ experiments, partly because these parameters also depend strongly on the state of the steel surface, consistency of the electrolyte, and other factors that are generally unknown for a system under study⁴¹, and partly because the experimental methodology induces errors in the measurement of these parameters^42,43. Thus, exchange current densities are in the literature generally defined as ranges, instead of constants^41,44,45.

On a more technical note, as the exchange current densities have a substantial influence in the modelling of localized corrosion⁴⁶, it is a major advantage that we are able to define these parameters as prior probability distributions, instead of as constants. An erroneously assumed exchange current density could result in a false sense of security, for a wrong corrosion rate. In the example presented here, this was evident for the exchange current density of iron (i_0,Fe). During the MCMC sampling a clear trade-off between the width of the anode W_an, and this exchange current density, was observed (see Fig. 7a).

Fig. 7: 2D histograms of the sampled model parameters.

a The width of the anode versus log10(i_0,Fe) versus the amount of samples. b The width of the anode versus ${\rm{log}}10({i}_{0,{{\rm{O}}}_{2}})$ versus the amount of samples. One sees a clear dependency between the size of the anode and the exchange current density for iron, while this dependency is not visible for the exchange current density of oxygen for example.

The example described here only considers the exchange current densities in the inversion. However, one could also think about taking into account the Tafel slopes for example. Although an increase in the number of parameters would increase the computation time, as well as increase the ill-posedness of the inverse problem, this approach is beneficial to identify which parameters may induce significant errors if defined inadequately.

The need for non-destructive testing methods to detect and quantify localized corrosion for structures where visual inspection is difficult, will keep increasing over the coming years, as engineering structures are ageing. ECT could prove to be an important addition to existing methods. At this moment, the only widely recognized and standardized non-destructive method is potential mapping, or half-cell mapping. Here, the electrical potential field at the surface of the porous medium is measured, to identify possible corrosion sites⁴⁷. Using the same kind of measurements, ECT would need little adjustment of existing techniques and apparatuses to be applied in the field.

While the Bayesian approach has various advantages, as discussed previously in this paper, the main disadvantage is its higher computational cost, compared to deterministic approaches, such as used in earlier studies^30,31,32,33. The Bayesian inverse method applied, requires sampling using Monte Carlo methods, resulting in larger computation times. Thus, the instantaneous determination of corrosion locations and rates in the field may currently be a too ambitious goal for the method. However, we are confident that steadily increasing computational power, as well as improvements of the algorithm and programming model, will enable on-site measurements and model analysis in the future. Our convergence diagnostics indicate that the number of samples could be reduced by 50%. Furthermore, a change of the programming model from COMSOL and Matlab to a compiled language that is more suitable for high-performance computing and parallelization (e.g., C++ and Cuda), is expected to substantially reduce the computing time from several days on a single core.

In summary, as a proof of concept, the application of ECT for the example of localized corrosion in reinforced concrete, has shown that even with limited prior knowledge, ECT is able to reliably identify the corrosion state of the reinforced steel. ECT allows us to not only accurately estimate the corrosion current in the case of localized corrosion, but also to identify the location and size of the anodic surface. Therefore, the direct quantification of the corrosion current density (the corrosion rate) is possible. We consider this an important advantage over traditional electrochemical methods, such as the linear polarization resistance method, where a major limitation is the need for assumptions to be made related to the anodic surface area. Furthermore, the technique also yields insight into fundamental corrosion science parameters describing the polarization behaviour of the reinforced steel, such as the anodic and cathodic exchange current densities. Lastly, thanks to the stochastic approach, ECT allows us to quantify the uncertainty inherent to the obtained solution. This is considered another advantage over traditional electrochemical methods, where no information on the reliability of the actual measurement is available.

Future research should address the optimization of the method concerning the measurement setup, including the position of the reference and counter electrodes, and the sensitivity related to different values for the parameters describing the polarization behaviour of the metal. Finally, validation against experiments and in the field will be essential.

Methods

In the following section, we describe the general approach, valid for any corroding metal in a conducting porous medium. The forward and inverse problems are described, as well as the applied MCMC method, used to sample the solution of the inverse problem. Finally, we show the application for the example of localized corrosion in reinforced concrete.

Forward problem

Macro-cell, or localized, corrosion occurs due to the galvanic coupling between a small active surface (anode) and relatively large passive surfaces (cathodes) on scales up to several metres. Numerical modelling of macro-cell corrosion is well established in different fields, such as corrosion of reinforced concrete and buried pipelines, and is able to accurately reproduce experimentally obtained data⁴⁸. These models have proven to be a useful means to investigate and predict the influence of, for example, the geometry, electrolyte resistivity, and Tafel slopes on the corrosion behaviour of the metal in certain environments^7,46,49,50.

Macro-cell modelling of corrosion involves the solution of Laplace’s equation for the electrical potential, applying Butler−Volmer kinetics⁵¹ to derive boundary conditions at the steel surface. Consider, for instance, a buried metal pipe, with a distinct anodic surface, surrounded by cathodic surfaces (Fig. 8). Assuming an isotropic and homogeneous conductivity, the distribution of the electrical potential, ϕ [V vs SHE], in the conductive medium (the electrolyte) is given by:

$${\nabla }^{2}\phi =0$$

(3)

Fig. 8: Macro-cell corrosion and the current lines for the example of a buried pipeline.

In red, the actively corroding (anodic) surface.

At all medium boundaries, except the metal boundaries, free surface boundary conditions apply ($\frac{{\rm{d}}\phi }{{\rm{d}}\vec{n}}=0$). At the metal surface, the boundary condition is described by Ohm’s law:

$$i=-\frac{1}{\rho }\frac{{\rm{d}}\phi }{{\rm{d}}\vec{n}}$$

(4)

where $\vec{n}$ is the normal to the metal surface and ρ [Ωm] is the electrical resistivity of the electrolyte (the inverse of the electrical conductivity). At the metal bar, the current density, i [A/m²], is proportional to the rate of the chemical reactions taking place at the metal surface, and is the sum of the current densities of all reduction and oxidation reactions.

Though many reactions are simultaneously taking place, we assume that at the anodic surface, the dominant oxidation reaction of the metal is:

$${\rm{Me}}({\rm{s}})\to M{e}^{n+}({\rm{aq}})+n{e}^{-}$$

(5)

This gives an outward current as indicated by the red arrow in Fig. 8. At the cathode, we assume the dominant reaction to be the reduction of oxygen, resulting in an inward current:

$${{\rm{O}}}_{2}({\rm{g}})+2{{\rm{H}}}_{2}{\rm{O}}+4{e}^{-}\to 4{\rm{O}}{{\rm{H}}}^{-}({\rm{aq}})$$

(6)

The current densities, or reaction rates, belonging to these redox reactions are described by the Butler−Volmer reaction. The oxidation of the metal is assumed to be activation-controlled. The relation between electrical potential and current density is therefore directly derived from the Butler−Volmer equation:

$$\phi ={\phi }_{{\rm{Me}}}^{{\rm{rev}}}+{\beta }_{{\rm{an}},{\rm{Me}}}* {\rm{log}}\left(\frac{i}{{i}_{0,{\rm{Me}}}}\right)$$

(7)

where i_0,Me [A/m²] is the exchange current density of the metal, ${\phi }_{{\rm{Me}}}^{{\rm{rev}}}$ [V vs SHE] the reversible potential, and β_an,Me [V/dec] the Tafel slope for the metal oxidation (see Eq. (21)).

The reduction of oxygen is assumed to be both activation- and concentration-controlled. At higher current densities, the diffusion of oxygen to the cathodic surface may limit the reaction. Using the Butler−Volmer equation and ’Ficks first law’⁵¹, the relation between electrical potential and current density can be described by

$$\phi ={\phi }_{{{\rm{O}}}_{2}}^{{\rm{rev}}}-{\beta }_{{\rm{cath}},{{\rm{O}}}_{2}}* {\rm{log}}\left(\frac{i}{{i}_{0,{{\rm{O}}}_{2}}}\right)-\frac{2.303RT}{nF}{\rm{log}}\left(\frac{{i}_{L}-i}{{i}_{L}}\right)$$

(8)

where ${i}_{0,{{\rm{O}}}_{2}}$ [A/m²] is the exchange current density of oxygen, ${\phi }_{{{\rm{O}}}_{2}}^{{\rm{rev}}}$ [V vs SHE] the reversible potential, and ${\beta }_{{\rm{cath}},{{\rm{O}}}_{2}}$ [V/dec] the Tafel slope for oxygen reduction (see Eq. (6)). Further constants R, T, n, and F are the gas constant, temperature, number of electrons transferred, and Faradays constant, respectively. Finally, i_L [A/m²] is the limiting current density, which is dependent on the oxygen concentration in the medium.

In our approach, we also include counter electrodes (CE), within our forward model formulation. These electrodes are used to apply external currents, I_CE, on the system, thus polarizing the metal. This polarization of the metal, using varying quantities of I_CE increases the quantity of observable data, thereby reducing the ill-posedness of our inverse problem (see section 2).

The complete strong-form formulations of the forward model are now obtained by substituting Eqs. (7) and (8) in Ohms law (Eq. (4)) for the anodic and cathodic surfaces, and using Ohms law for the CE surfaces:

$${\nabla }^{2}\phi =0\,\,\,{{\mbox{in the electrolyte}}}\,$$

(9)

with boundary conditions:

$$\frac{{\rm{d}}\phi }{{\rm{d}}\vec{n}}=-\rho * {i}_{0,{\rm{Me}}}* \exp \frac{2.303(\phi -{\phi }_{{\rm{Me}}}^{{\rm{rev}}})}{{\beta }_{{\rm{an}},{\rm{Me}}}}\,\,\,{{\mbox{at the anode}}}\,$$

(10)

$$\frac{{\rm{d}}\phi }{{\rm{d}}\vec{n}}=\rho * {\left[\frac{{\left({i}_{0,{{\rm{O}}}_{2}}\exp \frac{2.303({\phi }_{{{\rm{O}}}_{2}}^{{\rm{rev}}}-\phi )}{{\beta }_{{\rm{cath}},{{\rm{O}}}_{2}}}\right)}^{\gamma }}{1+{\left(\frac{{i}_{0,{{\rm{O}}}_{2}}}{{i}_{L}}\exp \frac{2.303({\phi }_{{{\rm{O}}}_{2}}^{{\rm{rev}}}-\phi )}{{\beta }_{{\rm{cath}},{{\rm{O}}}_{2}}}\right)}^{\gamma }}\right]}^{\frac{1}{\gamma }}\quad \,{{\mbox{at the cathode}}}\,$$

(11)

$$\frac{{\rm{d}}\phi }{{\rm{d}}\vec{n}}=-\rho * {i}_{{\rm{CE}}}\,\,\,{{\mbox{at the CE}}}\,$$

(12)

$$\frac{{\rm{d}}\phi }{{\rm{d}}\vec{n}}=0\,\,\,{{\mbox{at all other boundaries}}}\,$$

(13)

where i_CE is the current density at the counter electrodes (I_CE/Surface area CE). Note that Eq. (11) is an approximation of the inverse relation of Eq. (8), as described by Kim et al.⁵². Here, γ is a curvature-defining constant, controlling the sharpness of the curve in the transition from the activation-controlled to the concentration-controlled current. Kim et al.⁵² suggested γ = 3.

Theoretically, all parameters in Eqs. (9)−(13) are constants, depending on the environment and type of metal. As long as the size and location of the anodic surface(s) are known, the corrosion rate at the metal surface can be computed.

Inverse problem

Our inverse solver seeks to find the ’model parameters’, describing the size and location of the anode, and therefore the rate of the corrosion of the metal, using ’observations’ of the electrical potential at certain locations in the electrolyte. We take a Bayesian approach to the inverse problem, described in more detail by Earls^53,54.

We assume that the measurement of these observations is imperfect, which means there is an additive Gaussian random noise, η_i, in these measurements. Therefore, the observed electrical potentials may be written as

$${\phi }_{i}^{{\rm{obs}}}={\phi }_{i}^{{\rm{est}}}(\overrightarrow{\theta })+{\eta }_{i}$$

(14)

where ${\phi }_{i}^{{\rm{est}}}(\overrightarrow{\theta })$ are the modelled, deterministic potentials using model parameters $\overrightarrow{\theta }$. Where $\overrightarrow{\theta }=({\theta }_{1},{\theta }_{2},..)$, includes, in a simplified case, model parameters describing the position and the size of a geometrical shape approximating the surface of the anode.

As opposed to a deterministic approach, a stochastic approach expresses the data, and therefore the model parameters, as random variables whose probabilities explicitly encode their uncertainties. Instead of finding the ’best values’, one seeks to obtain the posterior probability, $p(\overrightarrow{\theta }| {\overrightarrow{\phi }}^{{\rm{obs}}})$: the probability distribution of the model parameters, $\overrightarrow{\theta }$, given the observed electrical potentials, ${\overrightarrow{\phi }}^{{\rm{obs}}}$. Following Bayes’ theorem, the posterior can be computed with the likelihood distribution and prior information concerning the model parameters:

$$p(\overrightarrow{\theta }| {\overrightarrow{\phi }}^{{\rm{obs}}})=\frac{p({\overrightarrow{\phi }}^{{\rm{obs}}}| \overrightarrow{\theta })p(\overrightarrow{\theta })}{p({\overrightarrow{\phi }}^{{\rm{obs}}})}$$

(15)

The likelihood, $p({\overrightarrow{\phi }}^{{\rm{obs}}}| \overrightarrow{\theta })$, describes the probability of having observed certain electrical potentials, ${\overrightarrow{\phi }}^{{\rm{obs}}}$, in light of the model parameters, $\overrightarrow{\theta }$. The prior probability, $p(\overrightarrow{\theta })$, encodes our a priori belief regarding the model parameters. $p({\overrightarrow{\phi }}^{{\rm{obs}}})$ is the ’model evidence’: basically a normalizing constant within Bayes’ theorem, to ensure the posterior is a proper probability density function. It is given by the integral $\int p({\overrightarrow{\phi }}^{{\rm{obs}}}| \overrightarrow{\theta })p(\overrightarrow{\theta }){\rm{d}}\overrightarrow{\theta }$, and therefore very time-consuming and difficult to compute in high dimensional parameter spaces. However, its computation is of lesser importance here, as its value does not affect the relative likelihoods of the different models.

Equation (14) implies that $p({\overrightarrow{\phi }}^{{\rm{obs}}}| \overrightarrow{\theta })=p({\phi }_{i}^{{\rm{obs}}}-{\phi }_{i}^{{\rm{est}}}(\overrightarrow{\theta }))=p({\eta }_{i})$. Approximating the noise distribution by a Gaussian with variance σ², we express the likelihood as:

$$p({\overrightarrow{\phi }}^{{\rm{obs}}}| \overrightarrow{\theta })=\frac{1}{\sqrt{2\pi {\sigma }^{2}}}{{\rm{e}}}^{-(1/2{\sigma }^{2}){\left({\overrightarrow{\phi }}^{{\rm{obs}}}-{\overrightarrow{\phi }}^{{\rm{est}}}(\overrightarrow{\theta })\right)}^{2}}$$

(16)

The prior probability densities allow us to include a priori knowledge or beliefs about the model parameters. In the least ideal case, we have very limited information, and can only define the prior with a uniform distribution with bounded support. For example, we can constrain the size of the anode, but we can hardly provide any additional information about the relative probabilities within these constraints. One could increase the a priori knowledge between the bounds by using, for instance, experimental data.

MCMC sampling

In this problem, we cannot analytically solve for the posterior probability. It needs to be evaluated point-wise, solving the forward model for each $\overrightarrow{\theta }$. For efficient evaluation of the marginal posterior distributions of each model parameter, we sample the posterior probability density by applying a MCMC algorithm. These methods construct a Markov chain containing sequences of the model parameters which are dependent samples of the required posterior distribution. Here we apply the Metropolis−Hastings algorithm^55,56. From the prior distributions, mean values are computed as starting values $\overrightarrow{\theta }(k)$, with k the iteration number being initially 0. The algorithm computes a candidate $\overrightarrow{\theta }(k+1)$ for each model parameter (j) using a uniform distribution, with a certain maximum stepsize (α_j):

$${\theta }_{j}(k+1)={\theta }_{j}(k)+{\alpha }_{j}* {\rm{random}}(0,1)$$

(17)

We then evaluate the ratio of the posteriors of θ(k + 1) and $\overrightarrow{\theta }(k)$. Evaluation of this ratio eliminates the ’model evidence’, $p({\overrightarrow{\phi }}^{{\rm{obs}}})$, and therefore the computation becomes relatively straightforward:

$$r=\frac{p(\overrightarrow{\theta }(k+1)| {\overrightarrow{\phi }}^{{\rm{obs}}})}{p(\overrightarrow{\theta }(k)| {\overrightarrow{\phi }}^{{\rm{obs}}})}=\frac{p({\overrightarrow{\phi }}^{{\rm{obs}}}| \overrightarrow{\theta }(k+1))}{p({\overrightarrow{\phi }}^{{\rm{obs}}}| \overrightarrow{\theta }(k))}\frac{p(\overrightarrow{\theta }(k+1))}{p(\overrightarrow{\theta }(k))}$$

(18)

Generating a random number u ∈ [0, 1], the candidate is ’accepted’, when r ≥ u. When r < u the candidate is ’rejected’, and $\overrightarrow{\theta }(k+1)=\overrightarrow{\theta }(k)$. Per iteration of the Markov chain, this computation is repeated in turn, for every single model parameter: sometimes referred to as Metropolis within Gibbs sampling. A “burn-in" phase is applied before the actual sampling phase. This is generally used to account for any initially bad “guesses" of the model parameters (i.e., highly improbable). However, in our case, we use it primarily to also ’tune’ the generation of the candidate values, by optimizing α(j) (Eq. (17)) to achieve a desired target acceptance probability. The approach adopted herein was described by Link and Barker⁵⁷. When the candidate is accepted, α is updated using:

$${\alpha }_{j}={\alpha }_{j}* (A-(A-1)\frac{k}{B})$$

(19)

and after the candidate is rejected:

$${\alpha }_{j}={\alpha }_{j}/\frac{1-{p}_{{\rm{accept}}}}{1-{p}_{{\rm{accept}}}(A-(A-1)\frac{k}{B})}$$

(20)

where A = 1.01 and p_accept = 0.42, to obtain an acceptance rate of 42%, as recommended in Link and Barker⁵⁷. The term $(A-1)\frac{k}{B}$ accounts for ’phasing out’ of the tuning process. Here k is the sample number and B the total amount of samples in the burn-in phase.

The example of localized corrosion in reinforced concrete

To show the application of the described approach we take the example of localized corrosion in reinforced concrete. This corrosion commonly occurs due to local breakdown of the passive layer, commonly due to chloride ingress⁵⁸. Galvanic coupling occurs between the relatively small active surface (anode) and the larger passive surfaces (cathodes), on scales up to several metres. It is especially dangerous due to the relatively small anode to cathode ratio, which causes a fast local decrease of the steel diameter at the anode. The quantification of local losses of steel, that is, the rate of localized corrosion, is thus important in the inspection of structures.

We consider a steel rebar in a concrete cylinder (Fig. 9). We assume that the main oxidation at the anode is that of iron:

$${\rm{Fe}}({\rm{s}})\to Fe^{2+}({\rm{aq}})+2{e}^{-}$$

(21)

Fig. 9: Schematic overview of the 2D axisymmetric representation.

The system is symmetric around the z-axis. In red the anodic metal surface, surrounded by cathodic surfaces (black).

Our 2D axisymmetric finite-element model, symmetric around the z-axis, is set up using the COMSOL Multiphysics^® software. The mesh consists of tetrahedral elements with varying sizes from small elements at the metal surface (<0.1 cm) to larger elements at the surface of the concrete (<0.5 cm). These sizes were obtained from a mesh convergence study, to ensure proper spatial discretization. On the boundary representing the steel, we assume there is one actively corroding surface (the anode) surrounded by passive steel (the cathode). Included are 2 counter electrodes (CE) at the surface of the concrete. The corresponding strong formulations of the forward problem are given by Eqs. (9)−(13).

It is difficult to find specific parameters for the boundary conditions at the anode and cathode in the literature. Especially for steel in concrete, these are difficult to determine experimentally, due to, for instance, pH dependence of the Tafel slopes⁵⁹. Values used within models of different researchers exhibit large scatter⁴⁶. Therefore, assuming all the parameters as constants is difficult to justify on the basis of literature data available. The exchange current densities are known to be able to vary over several orders of magnitude between different experimental measurements; thus embodying the largest source of error. Therefore, we include two extra model parameters. Besides the two parameters describing the location of the anode (i.e., the width W_an and the coordinate of the middle of the anode z_an (Fig. 9)), we also take into account the exchange current density of iron oxidation, i_0,Fe, and oxygen reduction, ${i}_{0,{{\rm{O}}}_{2}}$. All other parameters are set to constant values in this example, which can be found in Table 1.

Table 1 Parameters describing the polarization behaviour of the steel (for the anodic and cathodic surface) and model and medium properties of the forward model (see Fig. 9).

Under the assumption that we only have one anodic surface, and including the exchange current densities in the inversion, we have four model parameters $\overrightarrow{\theta }=({\theta }_{1},{\theta }_{2},{\theta }_{3},{\theta }_{4})$. These are, respectively, the width of the anode (W_an), the z-coordinate describing the position of the middle of the anode (z_an), and the exchange current densities of iron oxidation (i_0,Fe) and oxygen reduction (${i}_{0,{{\rm{O}}}_{2}}$).

The observed data vector ${\overrightarrow{\phi }}^{{\rm{obs}}}(x)$ consists of potential measurements at discrete positions at the concrete surface, x_i+kN = (i − 1)Δx; i = 1, . . . . , N, for different external polarization currents ${i}_{{\rm{CE}}}^{k};\,k=0,....,M-1$, where N is the cardinality of the data points at the concrete surface (L/Δx), and M the total number of different currents applied at the counter electrodes.

The prior information regarding the width and the location of the anodic surface was gleaned using experimentally measured potentials. Figure 10 shows a typical potential distribution measured at the concrete surface. We assume to have no information about the width of the anode, W_an, beyond that, for this setup, it will not be larger than the width of the curve at the median of the potential values (see Fig. 10). Therefore we determine a flat (i.e., uniform) prior probability with an upper boundary. The location of the anode is most likely around the minimum potential value measured. There is a very small probability that this location is outside of the range determined by max(W_an). Therefore, regarding the location of the anode, z_an, we assume a normally distributed prior, with the minimum potential as mean and the maximum width of the anode as four times the standard deviation.

Fig. 10: Estimation of the prior information for the anode based on potential measurements at the surface.

The sketch illustrates the determination of the prior information for the location of the centre of the anode (z_an) and the width of the anode (W_an), based on the measurable potential distribution at the concrete surface. W_an is described by a uniform probability density with 0 < W_an < max(W_an), where max(W_an) is determined as shown. z_an can be described with a normal probability function with the mean as shown here and a standard deviation of $\frac{1}{4}{\rm{max}}({W}_{{\rm{an}}})$.

Prior information about the exchange current densities is preferably experimentally determined for the specific situation. Because the exchange current densities are known to vary over several orders of magnitude, we consider the logarithm of the exchange current density. For this example we use the following log10 normal prior distributions: for log10(i_0,Fe) a mean of −3.25 and standard deviation of 0.88; for ${\rm{log}}10({i}_{0,{{\rm{O}}}_{2}})$ a mean of −6.98 and standard deviation of 0.96. These are based on experimental data measured in our laboratory in a Borate Buffer solution of pH 7.5.

To investigate and optimize ECT for the example of localized corrosion in reinforced concrete, we performed a series of synthetic tests: Synthetic observations are generated using the deterministic forward model instantiated using a set of known model parameters. These deterministic predictions were then contaminated with Gaussian white noise consistent with representative sensors. These synthetic data, ${\overrightarrow{\phi }}^{{\rm{obs}}}$, were then used as input data for the inverse model, which then seeks to uncover the input model parameters. Although experimentally obtained data will be preferred ultimately, these synthetic tests are useful to illustrate the approach. It affords us a good control in the initial assessment of the proposed method, without the need to deal with the difficulties associated with physical experiments.

All MCMC sampling chains contained 60,000 samples, obtained after 5,000 burn-in samples. The synthetic data (potentials) were obtained for input model parameters: W_an = 0.06 m, z_an = 0.40 m, log10(i_0,Fe) = −3.11 A/m² and ${\rm{log}}10({i}_{0,{{\rm{O}}}_{2}})=-7.5$ A/m². A normally distributed error with a standard deviation of 10 mV was added to the synthetic data, to account for an unavoidable measurement error during potential measurements.