Development of a Model for Experimental Data Treatment of Diffusion and Activation Limited Polarization Curves for Magnesium and Steel Alloys

Potentiodynamic polarization curves are often performed to characterize the bulk corrosion properties of a material in terms of its corrosion potential, rate, and Tafel kinetics in a given environment. This requires approaching them as one of two extreme cases, with purely activation or diffusion controlling the measured currents. For the intermediate behavior displayed by many real samples, making this assumption is a signiﬁcant source of error. In this work, a ﬁnite element model has been developed which considers both of these contributions simultaneously for experimental data with a broad material scope. This has been used to examine both extremely corrosion vulnerable samples in two magnesium alloys and extremely corrosion resistant samples in two steels. The effectiveness of several goodness of ﬁt statistics in comparing the similarity between experimental and simulated data has also been assessed.

Models of corrosion are ubiquitous in the literature. Mathematical formulations of corrosion processes have been developed to describe pitting, 1-4 galvanic, [5][6][7][8] and uniform corrosion [9][10][11] which in turn have been used to study factors influencing the rate of corrosion in specific systems including but not limited to chloride concentration in pitting initiation, [12][13][14] oxygen diffusion in kinetic limitations, 9,[14][15][16] and sample geometry in crevice corrosion. [16][17][18] These systems have been solved both analytically 4,8 and numerically. 4,5,7,15,17,19,20 Finding an analytical solution requires making simplifying assumptions such as a 1D-geometry, uniform surface reactivity, and negligible change in solution composition as a result of faradaic reactions. 8 This severely limits the scope of applicable systems. Numerical methods are by far the option of choice for any system with an added degree of complexity, whether that be a non-uniformly reacting surface, multi-dimensional geometry, or solution undergoing significant pH or other compositional changes which in turn affects the corrosion rate.
Numerical models of corrosion take one of two forms: (1) A model-first approach, 13,14,18 which begins with a theoretical system and looks to predict material evolution, or (2) An experiment-first approach, 9,17,21 which begins with a physical system and looks to explain observed behavior. Both approaches begin with model inputs for the initial concentrations and diffusion coefficients of species in solution, metal and solution conductivity, and Tafel slopes. These Tafel slopes (β) are extracted from potentiodynamic polarization curves (PDPs), which can also be used to measure the corrosion potential (E corr ) and corrosion current density (j corr ) of the system. In practice, extraction of these parameters from experimental data is problematic: the selection of the linear region is crucial for reproducibility and validity of analysis. If the linear region is selected based on a target R 2 value, the fitted region will be inconsistent between samples. Determining the potential at which to begin fitting is also problematic: Tait recommends linear fitting begin 50 mV from E corr and end after a one decade increase in current, 22 while Poorqasemi et al. suggest fitting can begin as far as 100 mV from E corr over the same change in current. 23 Many textbooks on the subject do not specify this range at all 24,25 as the linear region can change from metal to metal.
Significant discrepancies in the analysis method between materials exist: previous works have opted not to perform fitting, 26 perform a partial fit to a single branch, [27][28][29] or perform a full fit to both branches in both ideal and non-ideal cases. 23,30,31 If the two branches of the PDP are analyzed separately, they often yield different values for the corrosion current density. 32 Many use the value extracted exclusively from the cathodic branch as extensive surface changes and corrosion product formation occur in the anodic branch; however, the same issue may occur in the cathodic branch in acidic solutions or at potentials sufficiently negative to reduce the native passive film. 23 Fitting to the Tafel equation further assumes purely activation-controlled currents. If the reactions at the surface are under diffusion control, the fitting protocols discussed earlier cannot be used. Some work has been done to quantify systems of this nature involving iron 33 and copper, 34 but these have focused on a limited material scope or decoupled mass transport effects by introducing convection to the system. The agreement between corrosion current densities determined from both branches is also strongly dependent on experimental conditions used, such as the electrolyte and the scan rate. 32 All of these factors are significant sources of error for the PDP technique, which leads to poor agreement with other techniques such as mass loss and electrochemical impedance spectroscopy.
Herein, we present a material independent, finite element model that provides a systematic method for the more accurate data treatment of polarization curves. This is done through the incorporation of mass transport limiting current density to traditional Tafel kinetics. Experimental polarization curves are then fit through an iterative refinement process, with the experimental and simulated results statistically compared to validate model accuracy. This method is applicable to a broad material scope -specifically, extremely corrosion active magnesium alloys and extremely corrosion resistant steels have been presented here. As discussed previously, taking a numerical approach allows this model to be extended through the incorporation of complex, multi-component materials with asymmetrical geometries. Our model is therefore a significant improvement on those currently available in the literature, as ours is the first real example that offers reliable and reproducible comparison of the corrosion properties across such a broad range of materials regardless of the experimental approach.

Experimental
Materials.-Magnesium rods (99.9% purity, 5 mm diameter) were obtained from Sigma Aldrich. Steel A516 rods (11.3 mm diameter) and AM60 Mg alloys (1 cm × 1 cm) were obtained from General Motors Canada. Stainless steel 444 samples were obtained from Hy-droQuébec and machined into cylinders with 1 cm 2 surface area. All samples were fixed in cold mounting epoxy (Epofix, Streurs Canada) and the metal face was exposed using a TegraPol-25 polishing wheel (Streurs, Canada) through a series of successive grinding steps with increasingly finer grits of SiC foils (800, 1200, 4000). This was followed by polishing with a 1 μm diamond suspension and a 0.05 μm alumina paste using an MD-chem polishing pad to obtain a mirror finish. Solutions were prepared with reagent grade NaCl (ACP) and nanopure water (18.2 M -cm, Millipore).
Instrumentation.-Electrochemical measurements were made using an ELProscan 1 system (HEKA, Germany) with Potmaster software (version v2 × 66) and a VSP-300 system (Biologic, France) with EC-Lab (version 11.01). A 1M Ag/AgCl reference electrode in a Luggin capillary was prepared following a literature procedure. 25 A commercial SCE reference electrode was used for the measurements on stainless steel 444. A platinum mesh was used as the counter electrode in all cases. The finite element models were built using COMSOL Multiphysics version 5.2a equipped with the Corrosion Module.
Data collection and procedure for Tafel fits.-Polarization curves were recorded in 0.10 M NaCl over a range of −300 mV to +300 mV relative to open-circuit potential. Initial values of E corr , j corr , β a , and β c for use in the parametric sweep were extracted from experimental data using MATLAB 2015b. A linear fit was performed for each branch using a traditional Tafel fitting procedure where fitting began 50 mV from E corr and ceased when the current density had increased one decade from the start point. 22 These values were used as the initial guess of E corr , j corr , β a , and β c , which were then refined by simulating a parametric sweep and selecting the best statistical match. In cases where mass transfer limitations appeared to be present, a fifth parameter (j lim ) was included in the parametric sweep. The initial guess was then updated and the process was repeated using finer increments in the parametric sweep until further refinement yielded negligible improvement in the quality of the fit.

Results and Discussion
Factors influencing choice of model geometry.-The geometry of a corroding system is ultimately three-dimensional: uneven oxide formation dependent on the crystallographic orientation, heterogeneous corrosion rates and microgalvanic corrosion associated with the underlying microstructure, and localized processes such as pitting and crevice corrosion all alter surface topography. In the context of a laboratory measurement, metal samples are polished to a smooth finish prior to testing and the topography changes minimally relative to the overall sample dimensions when testing under mild conditions.
The standard electrochemical cell during experimental corrosion measurements employs a three electrode setup where the metallic sample of interest serves as the working electrode (WE). In a simulation, a combined reference and counter electrode (CE) can be used. An isolated region of the WE with well-defined surface area is exposed to the electrolyte through encasement of the sample in epoxy or the use of an O-ring. This geometry is depicted in Figure 1. This cell possesses an axis of symmetry through the center of the sample. This enables the original 3D geometry to be simplified to its 2Daxisymmetric equivalent ( Figure 1A). This reduces the computational power needed to solve the mesh and aids in convergence, which is necessary when additional partial differential equations (PDEs) with a spatial dependence are incorporated into the system -in this case, mass transport of additional species in solution. space. The dimensions were chosen such that the working electrode (WE) has an exposed surface area of 1 cm 2 , with m r = 0.564 cm, m h = 0.2 cm, c r = 0.564 cm, c h = 0.2 cm, ep r = 2 cm, ep h = 0.4 cm, el h = 2 cm. Tafel kinetics at the electrode/electrolyte interfaces in the WE and CE govern the currents measured at the WE. A boundary condition of electrolyte neutrality is applied to the walls of the container (dashed line). In versions of the model where the electrolyte composition is specified directly, a pH equilibrium is also specified within the electrolyte.
Simulations involving both 3D and 2D-axisymmetric geometries have been performed and are discussed in this work. For a description of the meshing technique employed for both geometries, refer to Figure S2. In general, the most appropriate geometry is determined by the goal of the simulation and computational power available to solve it. A 3D model is advantageous for asymmetric systems or localized corrosion processes, while a 2D-axisymmetric model is useful for a more detailed study of the individual reactions occurring during corrosion.

Decoupling activation and diffusion-controlled contributions to current.-Extraction of initial model inputs E
full description of the underlying equations used can be found in the Supplemental Information. In a finite element formulation, the total current at the corroding surface is calculated according to: Where i total is the overall current at the electrode surface, i loc,m is the faradaic current associated with reaction m, and i dl is the nonfaradaic current associated with double layer capacitance. The current distribution within the electrolyte is described by the following: i l = −σ l ∇φ l [3] Where i l is the current, Q l the charge, σ l is the conductivity, and φ l the electric potential within the electrolyte. At high overpotentials, the Butler-Volmer equation simplifies to a linear Tafel expression. Where η is the overpotential, R the universal gas constant, T the temperature, α the transfer coefficient, F Faraday's constant, and i corr , the corrosion current. More commonly, the corrosion current density (jcorr) is reported rather than the cTafel slope is reported rather than the transfer coefficient: Where β a and β c are the anodic and cathodic Tafel slopes respectively, and n is the number of electrons transferred during reaction. The overpotential may be calculated as: Where φ s , φ l , and E corr are the electrode, electrolyte, and corrosion potentials respectively. Numerical analysis of a standard PDP takes advantage of these kinetics to extract four corrosion-relevant parameters from the Tafel region: these are β a , β c , E corr , and j corr (where j corr is i corr /area). In systems under activation control, these four parameters are sufficient for describing the shape of a polarization curve. In systems under diffusion control, the reaction rate is limited by mass transport of reacting species to and from the electrode and so the current density plateaus at high overpotentials. Diffusion-limited currents are often observed while carrying out polarization measurements on a corroding sample where significant gas evolution occurs or extremely passivating oxide layers are present, 22 as these phenomena reduce the availability of electroactive species at the metal/electrolyte interface. Under these conditions, an additional parameter is needed to adequately describe the rate of reaction: j lim , the limiting current density. Depending on the underlying cause of this diffusion control, this may only apply to the anodic branch (j lim,a ) or cathodic (j lim,c ). The faradaic current at the electrode surface may therefore be calculated as: For a given reaction m, i loc is the total faradaic current, i lim the mass transport contribution to current, and i Tafel the activation contribution to current.
Optimization of kinetic and corrosion parameters.-A given PDP simulation was performed by linearly polarizing the WE from −300 mV to +300 mV relative to E corr at 0.167 mV/s. A multi-parametric sweep was then performed which took the initial values of β a , β c , E corr , and j corr extracted from 0 and examined the effect of changing these in coarse increments. All combinations were examined statistically and the best set of parameters was selected as the new input for a subsequent multi-parametric sweep using finer increments. In cases where diffusion control appeared to be present (initial values of |β| > 240 mV/decade, which would correspond to a process involving >4 electrons if it were purely kinetic controlled), traditional Tafel fitting procedures provided poor estimates of the corrosion kinetics. Therefore, a five-parameter sweep over a broader range was performed. The kinetics at the CE were assigned values corresponding to that of pure platinum: these are j corr = j exchange = 0.794 mA/cm 2 , E corr = +1.188 V vs. NHE, and idealized kinetics of |β| = 118 mV/decade, 24,35 To aid in convergence, a boundary condition is applied to the walls of the container (the region of bulk solution) of E electrolyte = 0 V. The simulated PDPs for four different metallic systems and the corresponding experimental data to which they were fit (Figure 2) demonstrate the broad material scope of the model, with excellent agreement found in all cases. The low E corr and high j corr of the two magnesium-based samples (Figures 2A and 2B) is consistent with their poor corrosion resistance. The j corr determined in this work is lower than those previously found in a Tafel analysis (80 μA/cm 2 in this work compared to 120 μA/cm 2 ): notably, this value was extracted using only the cathodic branch and under the assumption of purely activation controlled kinetics, 28 an analytical procedure commonly used for Mg alloys including AM60. 36 Similar β c values were found for the Mg alloys of ∼170-230 mV/decade; while a value for β a was not given for direct comparison, the steep slope in this region and sharp contrast in the shapes of the anodic and cathodic branches are consistent with polarization curves reported on Mg and zinc-containing Mg alloys. 37,38 Previous numerical analysis has been limited to the cathodic branch based on the analysis that since the hydroxide film formed on Mg and Mg alloys is only partially protective, there is no single anodic or cathodic reaction occurring in each branch and so corrosion parameters extracted from Tafel analyses can be problematic. 39 This highlights the challenge of correctly interpreting the meaning of these values after they are extracted. As is, the currents measured are a convolution of multiple reactions occurring at the metal-solution interface. Their respective contributions can be elucidated using the concentration profile of each species in solution ( Figure 4). This is particularly relevant for understanding the mass transport limited current densities that appear to be present in the cathodic branches for these samples, where hydrogen evolution is the dominant reaction. The aggressive formation of H 2 gas bubbles blocks the surface, effectively decreasing the surface area and resulting in a mass transport limited current density. No such limitation exists in the anodic branch where the dominant reaction is Mg (s) dissolution and so activation controlled reactions are observed.
In contrast to the Mg alloys, the dominant cathodic reaction for the steels tested ( Figures 2C and 2D) is oxygen reduction and no limiting current density is observed in the cathodic branch for these samples. It has previously been shown that oxygen reduction on passivated stainless steel is activation-controlled within the Tafel region. 40 Higher j corr values have been reported for S444 41,42 as only the anodic branch was extrapolated. The extracted parameters in this work are consistent with the difference in corrosion susceptibility between carbon and stainless steels. This is further evidenced in the appearance of a mass transport limited current density in the anodic branch of the stainless steel. This is characteristic of a diffusion-controlled process where transport of reacting species to the metal surface is slowed from surface blocking due to gas bubble formation or corrosion product restricting the electroactive species access to the metal surface. Such situations are typical explanations for diffusion-controlled processes occurring during the cathodic sweep of a PDP. As before, this deviation from an activation-controlled reaction rate poses a problem for traditional Tafel extraction that cannot decouple the two contributions.
In the anodic branch, however, the dominant reaction is metal dissolution. Therefore the limiting current density in the anodic branch for Figure 2D arises due to the extremely passivating nature of the surface oxide layer that controls the rate of dissolution. Activationcontrolled behavior in the anodic branch of the carbon steel alloy shows that this material does not successfully passivate under the presented conditions and the metal dissolution is kinetically limited.

Evaluation of statistical metrics for goodness of fit.-Examin-
ing the accuracy of a simulation with respect to theory or experimental data in a quantitative manner requires the selection of an appropriate goodness of fit statistic. Previous simulations of electrochemical experiments have considered a variety of metrics toward this end, including percent difference for chonoamperometry 43 and cyclic voltammetry, 44,45 mean deviation for microelectrode approach curves, 46 and least-squares residuals analysis for impedance. 47 Given the non-linear nature of a polarization curve, there are a number of fit statistics which might be considered as candidates for determining the optimal fit. Of those mentioned earlier, percent difference and residuals are metrics that have previously been used in this context. Deviation would be most applicable in a system where a model is being compared to a series of measurements rather than individual data sets as examined here. In addition to percent difference and residuals, the chi-squared goodness of fit test has been examined as a possible metric as it accounts for the number of parameters being fit.
In all tests, the largest discrepancy between experimental and simulated curves occurs in a narrow region around the corrosion potential where the model shows currents that are orders of magnitude smaller than experiment. This is due to the combination of two things. Firstly, when performing a linear sweep a potentiostat increases the potential in a discrete rather than continuous manner. This means that while a number of points immediately surrounding E corr will be sampled, it is unlikely that a data point will be collected at exactly E corr . It may be argued that measuring the initial open-circuit potential (OCP) of the sample serves this purpose, but this has two drawbacks: the OCP of an immersed sample often fluctuates even after an extended equilibrium time, and E corr is measured over an actively corroding surface undergoing surface changes. Secondly, the assumption that electrode kinetics are Tafel in nature is most valid at high overpotentials. For small overpotentials, the current rather than the log current is proportional to potential; this is the regime of polarization resistance. 29 For this reason, the model is not expected to be a good fit in this region and so goodness of fit statistics should account for this. To this end, a weighting function has been developed ( Figure 3A) to place higher value on statistical comparison in the Tafel region as opposed to the polarization resistance region. A narrow region of 10 mV has been chosen to investigate the difference between weighted and unweighted values while retaining as much information as possible over the entire curve.
The potential dependence of the three chosen statistics for a PDP performed on pure Mg (Figure 3) communicate the same information about regions where the fit is better or worse. In practice, using these metrics to compare two sets of parameters and draw overall conclusions about goodness of fit requires summarizing the values for each point of comparison as an aggregate (residuals, chi-squared) or average (% difference) as summarized in Table II. At this point, the differences between metrics become evident. The sum of residuals provides higher sensitivity for selecting between two possible sets of parameters. Percent difference is an intuitive metric which is independent of the order of magnitude of currents observed and thus generalizable across samples. The chi-squared sum retains the proportionality of percent difference, though it loses the sensitivity to sign with respect to the simulation over-or underestimating measured currents. When the chi-squared goodness of fit test was used to examine goodness of fit, fits that appeared both good and poor passed the test, suggesting this is not a suitable metric for this system. Residuals and percent difference are better characterizations of goodness of fit, the former for comparing two fits and the latter for comparing the final fit to experimental results.
As seen in Table III, the fit for pure Mg samples scanned at 0.167 mV/s ( Figure 2A) and 10 mV/s ( Figure 3) is an excellent match statistically (<3% difference) though the fits are very different. The lower E corr and higher β a obtained for faster scan rates are consistent with previous methods studies conducted by Zhang et al., 32 which explained this difference in terms of a disturbance in the charging current. If scanned too quickly, the double layer formed at the metal surface will not reach equilibrium before the measurement is taken and so the non-faradaic contributions to the overall current will be larger. The larger this effect (the faster the scan rate), the larger the change in E corr . This highlights one of the limitations of the model, which is reliance upon trusted experimental data. The extracted parameters are specific to a given data set, and so drawing conclusions about overall sample behavior still requires a consideration of the experimental conditions and associated studies on reproducibility.
Analyzing diffusion-limited current densities in terms of solution mass transport.-A full description of the underlying equations used can be found in the Supplemental Information. The basic version of the model described earlier is suitable for data treatment, with good comparison between Tafel parameters determined using traditional methods and those optimized via simulation for activation-controlled systems. Furthermore, the model provides access to these parameters in diffusion-controlled systems where traditional methods fail. To understand these diffusion-limited current densities requires a detailed description of the chemical reactions involved and the environment both prior to and during corrosion, with the concentration, charge, and mobility of all electroactive species and the supporting electrolyte explicitly specified (Table S2 in the supporting information.) The current distribution within the electrolyte is described by the following: Where i l is the current within the electrolyte, F is Faraday's constant, z i and R i are the charge and flux of species i respectively, and Q l is the charge within the electrolyte.
The flux of species 10 is calculated as the sum of contributions from mass transfer 11 and chemical Reaction 12 as described by: Where c i is the concentration, D i the diffusion coefficient, and u i the mobility of species i. With respect to reaction m, R i,m is the flux of species i, v i,m the stoichiometry coefficient, i loc,m the faradaic current, and n m the number of electrons transferred during reaction. The total current at the electrode surface is calculated in the same way as discussed previously.
To reduce the computational power needed to solve this system, a 2D-axisymmetric geometry was used to reduce the dimensionality of the model. The polarization curves simulated previously were rerun in the new geometry to confirm neither the geometry nor the physics change affected the measured currents.
Accurate mass transport within solution was validated through comparison of calculated transport numbers to literature values available for electrolytes of equivalent concentration. These transport numbers were determined from a simulated chronoamperogram performed at E corr + 300 mV. For a given ion in solution, this was calculated as the current associated with diffusive and migratory flux of this ion within the electrolyte domain relative to the total current within this domain. For an equivalent electrolyte concentration of 0.1 M NaCl, the calculated transport numbers are 0.3963 and 0.6037 for Na + and Cl − respectively compared to 0.3854 and 0.6146 in the literature; 24 these agree within 3%. The change in this values over the course of the experiment due to the formation of new charged species is negligible ( Figure S2).
The kinetics of reactions at the WE and CE are described using the same Tafel relationships as in 0 where the limiting current densities previously determined can now be analyzed in terms of the reactions occurring at the surface. All standard potentials given are versus NHE. 24 For the Mg-based samples, the reactions considered were: For the Fe-based samples, the reactions considered were: The cathodic reaction is determined by the range of potentials considered: Mg begins corroding at potentials where the rate of hydrogen evolution is much faster than that of oxygen reduction, whereas the opposite is true for Fe. Further reactions of the metal cations to form solid oxides/hydroxides have not been included at this stage of model development, but are planned for future works. To understand the source of anodic and cathodic limiting current densities, a simulation was performed using the kinetic parameters previously extracted for pure Mg and S444 (Table I). The concentration profiles of protons, metal ions, and dissolved gases were examined as a function of distance from the corroding sample in each case.  The concentration profiles calculated for a PDP on pure Mg (Figure 4) for E < E corr are important for examining the j lim,c observed in 0. When mass transport limited currents are present in the cathodic branch, possible causes include low availability of H 2 O to react at the surface or blocking behavior of the produced H 2 (g) . Negligible Mg 2+ has formed at these potentials and so transport of the metal ion through the oxide film is not a factor in this phenomenon. As the concentration of water as a solvent far exceeds everything else in solution and is essentially constant near the electrode, the most likely scenario is the produced H 2 (g) participating in surface blocking. This is consistent with the high activity of magnesium alloys toward hydrogen evolution. 48 From the concentration profile of H 2 (g) observed in Figure 4B, the region of high H 2 (g) concentration extends quite far into solution (mm scale) for a PDP performed at the slow scan rates suggested by standardized methods. 49 The size of this region will also depend on the scan rate and potential range employed: faster scans over a narrower potential range should see this effect to a lesser degree due to reduced hydrogen depletion near the electrode at the same overpotential. As the potential increases to give E > E corr , the trend is consistent with previous time points. This is inconsistent with the negative difference effect (NDE) observed on Mg alloys, where the rate of hydrogen evolution increases with current density during anodic dissolution of the metal. 48 This is due to the nature of the applied kinetics as Tafel in nature: a model which encompasses the NDE would require either an additional or altered kinetic description of the material behavior. Some first principles work to define this has been done recently, 50 but this requires an understanding of the mechanism involved which has been debated going back to the 1950's. 51 As mentioned previously, the equations considered in this simulation do not currently include the formation of magnesium hydroxide: [17] This assumption should only effect the concentration profiles observed for where Mg 2+ begins to be formed.
The concentration profiles calculated for a PDP on S444 ( Figure  5) for E > E corr are of interest for examining the j lim,a observed in 0. The mass transport limiting currents observed in this branch can be related to Equation 16 where Fe dissolution is the dominant reaction. Stainless steels are well known for their corrosion resistance in aqueous environments due to the high amounts of chromium present in the alloy. 52 Chromium forms a thin passivating oxide film over the metal surface. Fe dissolution is still the dominant reaction but the rate slows down due to hindered transport of species through this film. When a current plateau is observed in the anodic branch, this film is stable: current increases associated with pitting initiation and passive film breakdown are not observed until much higher overpotentials are applied than those simulated here.
The question is therefore whether it is the hindered transport of oxygen or of dissolved Fe 2+ through the film which has the largest impact on the currents measured. Figure 5 considers the concentration profiles of both these species in the absence of a passive film. Under Table II. Summarized statistics for the data fits shown in Figure 2 and Figure 5. U refers to unweighted values, W refers to values determined using a weighting function which reduces the weight of values determined for |η| < 10 mV.

Alloy
Scan    these conditions, the region of high Fe 2+ concentration is limited to a much smaller distance from the sample than that of low O 2 (g) . This region will be compressed if diffusion away from the surface is hindered. This has the potential to produce a steep concentration gradient, reducing the rate of further Fe dissolution.

Conclusions
The model developed in this work is able to treat polarization curves in a systematic fashion in order to extract information about the underlying corrosion properties. Traditional analytical procedures tend to treat the system as either activation-controlled and extract Tafel slopes or as diffusion controlled and extract a limiting current density; the ability to consider both contributions to current is useful for systems with intermediate behavior which do not behave as either one of these two idealized cases. This allows for the analysis of a broader material scope. In this work, both extremely active magnesium-based alloys and extremely passivated stainless steels, ranging over 1V difference in E corr and two orders of magnitude difference in j corr . The usefulness of common goodness of fit statistics has been assessed, with the sum of residuals and percent difference providing the strongest insight into the accuracy of the fits determined, with excellent agreement between simulation and experiment within 3% difference. The fits determined agree both with statistics and literature where comparison values are available, with similar Tafel slopes observed to previous studies on magnesium and similar corrosion currents observed to previous studies on stainless steel.
As presented, the model neglects the formation of solid species. Future work will build upon the model developed here in order to include the formation of oxides and hydroxides as a corrosion product, accounting for the corresponding effect on surface topography. The model may also be extended to examine localized corrosion processes such as pitting or micro galvanic corrosion.