Evaluation of pH variation in cathodic protection conditions by FEM simulation

Cathodic protection is an electrochemical technique used to control the corrosion of metals. It does this by sending a proper amount of current equal to the current required by the cathodic processes at the protection potential. One of the chemical effects of cathodic protection is a local pH increase at the surface, which is an important factor in the reduction of corrosion rate, as carbon steel may work in a passive condition. In the present study, a tertiary current distribution finite element method model is introduced. A dynamic boundary condition is considered, by defining a Fe anodic Tafel slope as a pH function, and oxygen reduction reaction limiting current density as a function of available oxygen at the cathode surface. The model has been validated by direct pH measurement during a short‐term polarization test, and subsequently used to predict the pH and the available oxygen at the cathode surface for long‐term condition.


| INTRODUCTION
Cathodic protection (CP) refers to the polarization of local cathodes to the most active anodes. [1,2]In this circumstance, the potential difference of local microcells, the driving force of corrosion cells, would be extremely reduced or eliminated.[5] The first effect, oxygen consumption, is linked to the limiting current density of oxygen reduction reaction (ORR) [2] : 8][9][10] The limiting current density for ORR, which is the ratedetermining step of the cathodic reaction in the reasonable protection potential range, could be considered as the protection current density. [2,7,11,12]15][16] Another possible reaction at the surface is the hydrogen evolution reaction (HER) which is the predominant reaction at more negative potential and could be written in two forms: (2) The behaviors of each mentioned cathodic partial reaction, ORR and HER, as a function of time was modeled with exponential functions. [8,9]The core idea of using the exponential function for time-dependent reactions was calcareous deposit formation and expected changes in the diffusion of the ions.][29] pH variation from one aspect is a result of the chemical balance between hydrogen ions and hydroxyl ions and it is linked to the electrochemical reactions from another aspect.For finite element method (FEM) simulation of CP, and evaluating pH variation at the cathode surface, a comprehensive model is needed.This model should include (1) electrochemical polarization of the cathode, including ORR, HER, and Fe polarization, which links to (2) diffusion of oxygen and other ions under the concentration gradient and (3) mobility of ions in the electrolyte inward-outward of electrodes under the potential and concentration gradient.Moreover, (4) electrolyte neutrality equilibrium is essential.A model containing all four parts can predict a potential and pH stability point.Understanding the pH effect and its stability at the surface is applicable in analyzing CP phenomena in soil and concrete [30,31] and predicting the effect of interference, [29,31] intermittent CP, [27,28] and some proposed mechanism of CP of structure under disbonded coating. [32]he stability condition with specific pH and potential trigger chemical reactions with mutual interactions to form stable and unstable compounds. [33]The formation of these groups of compounds, which are not entirely understood yet, is not included in the present study model.Moreover, a high pH environment, low oxygen concentration, and coating type could lead to different types of corrosion, [34][35][36][37][38][39][40] which is also omitted from the scope of this study.In the present study, as an electrochemical reaction, ORR, HER, and Fe polarization are included in the model, and water dissociation is included as a chemical stability reaction.Moreover, instead of the exponential function for these three electrochemical reactions, [8,9,15,19] more generic assumptions are introduced to the model.This work has two main parts.First, a proposed model is fitted to some experimental data, then long-term data are extracted from the model.For proofing the model at the first step, pH measurement has been done simultaneously.

| EXPERIMENTAL
First of all, the polarization and pH near the surface of the carbon steel sample in the simulated soil solution were recorded.Then, a mathematical model was built, and the gained polarization from the experiment was introduced to the model as a boundary condition.The relationship between pH and anodic Tafel's slope was another variable extracted from a potentiodynamic test series in NaOH solution and introduced as part of the model boundary condition.Other general variables, such as the diffusion coefficient of oxygen and other species, were defined from other references.The polarization steps similar to the experiment were applied to the model, with the pH variation data being extracted from the mathematical model thereafter.After finding the best boundary condition for fitting the model to the experimental result, the model was used to simulate and predict cathode behavior for a more extended period.

| Experimental study
Potentiostatic polarization tests and pH measurements have been done simultaneously.The test solution was 700 mL distilled water containing 0.2 g/L NaCl and 0.2 g/ L Ca 2 SO 4 .2H 2 O.The electrolyte resistivity was about 10 Ω m, and the pH was around 6.8-7.The silver/silver chloride/saturated KCl (SSC) reference electrode is used for all the polarization tests.For the polarization test, an AUTOLAB PGSTAT30, and for pH measurement, a pH meter made by DeltaOHM© (model HD 8705) with a resolution of 0.1 pH is used.The pH electrode, produced by Mettler-Toledo©, was inserted into a cylindrical plastic shaft containing a gel electrolyte and connected to the pH meter by a fixed cable.The sample was carbon steel (Grade C, ASTM A 283/283M). [41]The distance between the tip and the sample is fixed at around 40-100 µm (Figure 1).Lowering the distance may lead to an increased chance of blockage of the pH probe tip with corrosion product and more than this may lead to measuring the electrolyte pH, not the local electrode pH.By changing the direction of the microscope light, the real surface and shadow were separated, so the surface of the sample to the pH probe tip is highlighted with a black line in Figure 1.
The test solution was exposed to the air.Potential steps were 50 mV and start with corrosion potential to 700 mV toward the cathodic polarization with a 200 mV/ h scan rate.In each potential step of the polarization test, data were recorded at minutes 1, 3, 7, and 15.The reported current density in Figure 2b is a value of 15 min, before shifting to the next polarization level.The pH probe tip was a semisphere, and the sample was a small rod.The sample holder body integrated with the pH probe has a small degree to allow the produced bubbles to leave the surface.Moreover, the pH probe tip, sample holder, cell, and solution were transparent, and the distance between the sample and probe tip was visually checked to verify the absence of gas blockage.The distance between the pH probe tip and the reference electrode to the sample was the minimum that is practically suitable and possible, so the local values are measured.In this circumstance, the factor shape of the sample was minimized, and the two-dimensional (2D) model could be successfully used for the simulation.
The polarization graph extracted at the sample close to the pH probe, that is, 40-100 µm, is presented in Figure 2a, which also presents the fitted graph for extracting electrochemical parameters.The current density of the electrode in cathodic polarization versus pH is presented in Figure 2b.
Another variable in the model is the relation between pH and the anodic Tafel slope.For this purpose, another potentiodynamic polarization test, with a 0.6 V/h scan rate, was done for different alkaline solutions containing sodium hydroxide (NaOH).Samples were immersed in the solution for a suitable time according to pH value.Suitable immersion time for pH equal and more than 11 is 3-5 days and 1 h for lower pH values. [30]The exposed surface was 1 cm 2 in the PVC sample holder.Potential was scanned from −0.1 V of free corrosion potential to 0.6-0.8V versus SSC, where a breakdown in polarization happens, and a sharp increase in the current density appears on the polarization graph.The relation between pH and Tafel's slope of the sample anodic branch and the selected value for the modeling are presented in Figure 3.

| Approach
Modeling the CP with all the effective parts is highly nonlinear and complicated.The present model consists of: I.Chemical species in the electrolyte, comprising of H + , OH − , Fe 2+ , and O 2 .It is assumed that H 2 gas is freely moved in the electrolyte and not affected by other variables.II.Polarization of the cathode and anode.Cathode polarization is modified with available oxygen at the electrode surface.III.Fe corrosion, ORR, and HER at the cathode surface are linked to their polarization and current output.In some ranges of potential, oxygen is consumed, and OH − is produced, and in higher polarization potential, H + is consumed to produce hydrogen molecules.IV.H + production from reverse HER at the anode surface.The portion of HER is linked to the potential, as for the H + reaction too.V. Diffusion of oxygen towards the cathode surface from the reference point (constant concentration point) under different concentration driving force with a time variable gradient.Oxygen diffusion generic function, that is, Equation (10), is used in the modeling, not a simplified or specific equation for specific conditions. [13]I. Diffusion of H + , OH − and Fe 2+ under different concentration driving forces between anode and cathode.VII.Under the potential gradient in the solution Fe 2+ , H + , and OH − ions migrate and move based on their charges.VIII.The solution has a defined resistivity.
IX.There is an equilibrium between H + and OH − based on the water dissociation reaction.X.The solution is electroneutrality stable.So, the model has three main parts; electrochemical reaction at the surface, the motion of species, and equilibrium chemical reactions.In the present modeling, a 2D model was built with a dimension of 0.2 × 0.15 m.The surface of the cathode is 0.02 m, and the anode has a 0.13 m length and 2 × 0.13 m surface.The surface of the cathode is a straight line, and the shape factor is not considered in the modeling. [42]The mesh size should be small enough to gain consistent results, especially when diffusion is included in the model.With a very small mesh size, the calculation time increases too much, without improvement in the accuracy of the data.In the present study, the mesh size varies from 0.025 to 7.4 mm.At the electrode surface, five mesh layers with the smallest size are considered (boundary layer).Each simulation time with PC (Core i5, 2.3 GHz process, and 8 GB RAM), and 1 s time step, takes up to 8 h.The model is schematically presented in Figure 4.The centerline, which connects the center of the cathode to the anode, is used for extracting some data.

| Electrochemical reactions
Iron corrosion, ORR, and HER are the three main equations.The equilibrium potential of all three equations (E Fe eq , E O eq 2 , and E H eq 2 ) are considered constant at pH = 7 as presented in Table 1, so: where where E is electrode potential, Tafel slope (β Fe ) is a function of pH, as presented in Figure 3.After producing the ferrous ions at the surface, pH and other available ions incorporate to form some stable and unstable compounds but are not included in the model.For ORR which is presented in Equation From the electrochemical point of view of the ORR, available oxygen will affect the current of the reaction in both the activation and diffusion control regime. [10,16]In diffusion control polarization, the current of the ORR reaction becomes equal to the limiting current density, i L .It is possible to calculate i L : and for HER In the normal range of variables, the HER will be under activation control, and the pH will change the equilibrium potential of the reaction.In the present model, it is assumed that hydrogen gas without any interaction could be removed from the reaction surface.
The polarization behavior of carbon steel is a nonlinear function in the modeling.As primary boundary conditions, the experimental results of Figure 2a are fitted to Equations ( 5), (6), and (8).The polarization presented in Figure 2a is used as the base boundary condition for cathode polarization behavior.Moreover, a combination of Tafel's slope dependency on the pH and varying the limiting current density of ORR as a function of oxygen availability makes polarization behavior dynamic.Polarization is the main part of the model,  which is dynamic and entirely linked to other phenomena such as the diffusion of species and their reaction in the solution and surface of the electrode.

| Motion of species
A well-known equation for describing the motion of species in the electrolyte under CP is the Nernst-Planck equation.This equation has three main parts: Diffusion: It could be explained with ∇ D c i i , this means the diffusion of the species is proportional to the gradient of concentration multiplied by the diffusion coefficient.The driving force is the gradient of concentration.
Mobility: It is defined by where u m i is the mobility and calculated from Nernst-Einstein relation- In this equation, mobility is proportional to the concentration and potential gradient.Here, the driving force is the gradient of the potential.For species like oxygen molecules that do not have an electric charge, mobility is not an effective parameter in fluxes (charge z i = 0).
Convection: It is presented in the form of c μ i , where μ is convection speed.Convection, which is proportional to concentration, is free movement and generally is against the flux of species.A preliminary study shows that considering the convection in the modeling, it could improve the homogeneity of the electrolyte pH, but not at layers close to the surface.To reduce the cost of simulation, the convection was eliminated from the model.Moreover, this concept will be used for the soil environment so that the convection could be eliminated.
Hence, for flux for each species, the used Nernst-Plank equation for diffusion and mobility would be: for oxygen: and for species with electric charge:

| Chemical equilibrium reaction
Different species and intermediate components could be considered in the steps of Fe corrosion.In the present study, the equilibrium equation of water dissociation is added to the model: where k log( ) = −14.

| General reaction rate equation
In any case, the solution should be electrically neutral for species with the charge, that is, Fe 2+ , H + , and OH − : Generally, the proposed mechanism for the presented condition should consider chemical and electrochemical reactions beside the flux of the species, so for every species: Above equation is the general term that includes the flux of all species (Equations 10-13) and their equilibrium reaction (Equation 14).Water dissociation is assumed as an equilibrium reaction, which links H + and OH − concentrations in all parts of the solution domain and electrode surfaces.It is assumed the hydrogen bubbles could leave the surface freely, so hydrogen gas diffusion is not included in the model.
Based on Equation 4, transforming Fe to Fe 2+ or producing Fe 2+ means corrosion, at the rate of production of the ferrous ions in time, could be considered as corrosion rate.The general scheme of the model is presented in Figure 5.

| RESULTS AND DISCUSSION
The experimental results of Figure 2a, which is used to define the model's boundary condition, and Figure 2b, which is representative of the current density versus pH, show a correlation.In the first step of polarization, where the activation polarization is predominant and the pH changes from the natural values, there is a linear correlation between cathodic current density with both potential and pH (Figure 2a,b).In the concentration polarization regime, where the changes in the potential did not change the current density, pH becomes more constant.At high negative potential values, where the HER is predominant, again the current is increasing, more than 1 A/m 2 here, showing a continuous increase in pH.
A potential profile similar to the potentiostatic test, 0.05 V step every 15 min, is applied to the model.In the simulation, boundary conditions could be constant or dynamic.In dynamic boundary conditions, one or more variables of the model change with the progress of some phenomenon.The time step in recalculating the dynamic boundary condition is 1 s, and data are reported every 6 s (0.1 min).Comparing experimental pH-polarization data with the simulations with constant and dynamic boundary conditions is a validating step for the model.
After validating the proposed model, some long-term constant current, 100 days, is applied to the model, and pH variation and oxygen concentration are extracted at the cathode surface.For better understanding, changes in the polarization behavior of the cathode with a 1-day time step in the current sweep, which is affected by dynamic boundary conditions, are presented and discussed too.All the reported data of the cathode are collected at the center point of the electrode.

| Cathodic polarization modeling
In polarization modeling, the potential of the electrode was changed with time.So, the potential time is an input of the model, and the current density at the cathode surface is one of the outputs of the modeling.[45] This comparison shows that the model follows the defined polarization correctly.In the dynamic boundary condition, the polarization behavior is a function of oxygen and pH at the surface.So, imposing local potential and current on the first polarization boundary condition is not useful.
An important parameter in the analysis of dynamic polarization conditions is pH.In the present study pHpotential graph, that is, the Pourbaix diagram, is used for validating data of the modeling in the dynamic polarization behavior (Figures 6-9). [12,26,46]The Pourbaix diagram depicts the thermodynamical state of a metal in an aqueous electrolyte as a function of pH and potential.In the present study, the Pourbaix diagram is used to show the expected state of the cathode surface, even if it has not reached the final stability condition.Moreover, pH versus time is used to evaluate the time effect.49][50]  The boundary conditions used for the four following simulations are presented in Table 1 and Figure 2a, except for the following details; simulation 1 is extracted with the polarization behavior, as presented in Figure 2a, and values of Table 1 as boundary conditions; in simulation 2, ORR cathodic branch is modified during the progress of simulation with Equations ( 6) and (7).The constant values of Equations ( 6) and ( 7) are presented in Table 1, while the remaining part is constant and similar to simulation 1.In simulation 3, the anodic Tafel's slope of Fe is modified as a function of pH, Figure 3, and all other parts, such as ORR, are constant and similar to simulation 1.Finally, in simulation 4, the boundary condition is dynamic in both the ORR cathodic branch and anodic Tafel's slope of Fe.In simulation 4, the ORR reaction will follow Equations ( 6) and ( 7), and the anodic Tafel's slope is extracted from Figure 3 and is time-dependent.
In Figure 6, a fixed polarization is assumed.In fixed polarization, the relation between current and potential is fixed and similar to Figure 2a.pH could vary as hydroxyl ions produce and link to the current.Therefore, the changes in the pH could be linked to potential with a fixed polarization behavior, as presented in Figure 2a.In the range of diffusion polarization in Figure 2a, the current density that passed through the interface is constant, so the rate of the reaction that needs electrons becomes constant.This constant current rate leads to a constant pH, which is like a vertical line in the Pourbaix diagram.Pourbaix diagram, Figure 6a, is used for comparing the theoretical stable condition and the measured pH-potential values; however, the tested conditions are not entirely on "stable condition," as the Pourbaix diagram is.In simulation 1, the stable pH, which is a reasonably vertical line, is around 12. This pH is a little higher than the measured values, 11.
In simulation 2, Figure 7, the ORR cathodic polarization based on Equations ( 6) and ( 7) is modified in time.After a few steps in the potential and moving to the cathodic branch of the model, the oxygen in the vicinity of the electrode is consumed and the assumed ORR limiting current density is proportionally decreased.The stable pH at the second half of the test, after 120 min, reaches 11, which is close to the measured value.
In simulation 3, Figure 8, the anodic branch as a function of pH is added to the model, and the ORR polarization is considered a fixed condition.As illustrated in Figure 8a, at potential around −0.6 V versus SSC, the changes of the pH are gradual, with respect to the previous condition of Figures 6 and 8.This is a result of the modification effect of the anodic Tafel slope.As presented in Figure 8b, pH shows a fast increase in the first minutes, so this effect in the anodic Tafel's slope is significant (Figure 3).The range of potential in polarization that is affected by the anodic Tafel's slope changes is limited in this case and is just close to the corrosion potential.
In the last simulation, simulation 4 in Figure 9, both modification of ORR and the anodic Tafel's slope by time are implemented in the modeling.The pH is increased fast in the first few minutes too.So, based on the defined boundary condition, the anodic Tafel's slope will increase suddenly.After the consumption of oxygen in the vicinity of the cathode, the ORR cathodic branch is modified and shifted to the lower values.Decreasing limiting current density could expand the potential range affected by Tafel slope modification.On the other hand, reducing oxygen in the vicinity of the cathode and increase in the pH has a synergic effect on increasing their effects too.Despite the potential shift between experiments and the simulation, the pH-potential in the complete dynamic boundary condition shows similar behavior.Implementing time in Tafel's slope-pH function makes the model closer to the experiment.The result of the fourth simulation and slower changes in the pH in experiments may be linked to the high-pH film formation at the surface, which is not in the scope of the present study.
In Figures 6a, 7a In the experimental results, pH during the first 60 min, where the cathodic polarization is around 0.2 V, does not change too much.While in all four simulations after the first potential step, the pH increased to values around 11-12.It may have two reasons, non-Faradic reaction at the surface and probe distance to the surface, which needs further analysis.[53] Based on Figure 9 results, the assumed formed film is expected not to allow the pH to change rapidly and simultaneously facilitate the HER reaction and reduce its activation energy.This film needs further investigation.slight increase in the potential.pH, as presented in Figure 10b, continuously increased to the value of 13 after 20 days, and after that, the slope of pH variation is reduced.After a few days, the oxygen concentration comes close to zero, but longer simulating times show that the oxygen could find the opportunity to diffuse to the surface of the cathode.The sign of this diffusion is an increase in oxygen concentration to around 0.5 mol/L in Figure 10c.The potential gradient in this circumstance is constant, so this diffusion and increase in oxygen concentration is under the gradient of concentration, ∇ D c i i .This increase, even after 100 days, is around 20% of the general reduction of oxygen concentration to around zero.Reduction in oxygen concentration in the vicinity of the surface and as a function of cathodic polarization potential was previously reported, which is in general agreement with the simulated results. [49]other long-term extracted parameter is pH at the cathode centerline, presented in Figure 11.pH at the cathode surface increases over time and decreases at the anode surface.Consumption of oxygen at the cathode will lead to the formation of a high pH environment around the cathode, and this area is expanded.Simulation results show that after more than 10 days, the high pH zone could come close to the anode.The hydroxyl ions are produced at the cathode and increase the pH.Because of its charge, it moves toward the anode.In the present model, cell size is small, and the slope of pH changes could be high and nonrealistic in a long-term study.

| Long-term modeling
As presented in Figure 10c, oxygen concentration will be reduced drastically over time in this model, while pH increase needs appropriate time.For a better analysis of the model, cathodic current steps are applied to the model, and potential response is extracted.Local potential versus current density at the center of the cathode is presented in Figure 12.Each current step is applied for 1 day, and potentials are extracted from a short-time and 1-day time span.This means the potentials are extracted at the beginning and the end of each current step.
In dynamic boundary conditions, ORR limiting current density will change by available oxygen, and the anodic Tafel's slope is a function of pH and is increased by increasing cathodic polarization.The polarization of Figure 12a after a short time is similar to the schematic of polarization B in Figure 12b, which is modified with a reduction of limiting current density to 30 mA/m 2 .
Increasing pH to a range of 12 or more could increase the anodic Tafel's slope to more than 0.9 V/decade, Figure 3 and Figure 10b.This effect is schematically presented in Figure 12b  of 0.94 V/decade.A combination of these two effects schematically is presented in Figure 12b, Polarization C, which is similar to the 1-day polarization of Figure 12a.The schematic of Figure 12 shows that if the limiting current density or ORR does not apply in the model, slope modification of the anodic branch is not an effective parameter in the protection potential range.The pH measurement versus potential in the vicinity of the cathode confirms this combined effect on the real polarization of the cathode.

| CONCLUSIONS
A comprehensive model for FEM simulation of the electrode under CP is presented.Complete polarization of the electrode, motion of the species, and water equilibrium are considered in the modeling.Moreover, a dynamic boundary condition is assumed, that links the anodic Tafel's slope to the pH, and the ORR limits current density to oxygen availability.Complex film formation was not included in the model.The study shows: -Changes in the oxygen concentration at the surface are fast, and while the CP system is working, oxygen compensation is much faster than its diffusion rate.-pH increased during the CP current application to more than 11.-Reduction in limiting current density and pH effects of the anodic Tafel's branch could change the real polarization behavior of the electrode.

F I G U R E 1
Distance between sample and pH probe tip (40-100 µm).[Color figure can be viewed at wileyonlinelibrary.com]F I G U R E 2 (a) Potentiostatic polarization, (b) cathodic current density versus pH measured from 40 to 100 µm of carbon steel sample in simulated soil solution.
Anodic Tafel slope of carbon steel in NaOH solution.F I G U R E 4 Anode and cathode arrangement of the model.[Color figure can be viewed at wileyonlinelibrary.com]

T A B L E 1
Coefficient and constant values.

F
I G U R E 5 Graphical presentation of the modeling.[Color figure can be viewed at wileyonlinelibrary.com]F I G U R E 6 Carbon steel in the simulation solution in the potentiostatic test, comparison with simulation 1 with fixed B.C. (a) Pourbaix diagram and (b) pH-time.F I G U R E 7 Carbon steel in the simulation solution in the potentiostatic test, comparison with simulation 2 with oxygen reduction reaction modification effect in B.C. (a) Pourbaix diagram and (b) pH-time.F I G U R E 8 Carbon steel in the simulation solution in the potentiostatic test, comparison with simulation 3 with anodic Tafel-pH function and fixed oxygen reduction reaction modification effect in B.C. (a) Pourbaix diagram and (b) pH-time.

F I G U R E 9
Carbon steel in the simulation solution in the potentiostatic test, comparison with simulation 4 with anodic Tafel-pH function and oxygen reduction reaction modification effect in B.C. (a) Pourbaix diagram and (b) pH-time.
, 8a, and 9a, the Pourbaix diagram is used as a background to show possible stable cathode surface state, even if the test condition is transient, so the stability zones of the diagram should not be considered a definite state of the cathodic surface.The simulated pHpotential profiles versus the measured value from Figures 6a to 9a show a better correlation in dynamic boundary conditions, confirming a strong improvement in the simulation.
With dynamic boundary conditions, ORR limiting current density modification and Tafel's slope relation to the pH, long-term changes in the pH and the oxygen content at the surface of the cathode are studied.For example, the simulation results of the mentioned geometry under 0.31 A/m 2 cathode current density for 100 days are presented in Figure10.The boundary conditions used for Figures10 and 11are similar to the fourth dynamic model, Figure9.The simulation results in Figure10ashow a reduction in the potential on the first day.From the initial days up to 100 days, there is aF I G U R E10 FEM simulation in constant current mode and dynamic boundary condition with I cathode = 0.31 A/m 2 : (a) Potential (b) pH, and (c) oxygen concentration at the center of the cathode.
as a red dashed line with a slope F I G U R E 11 FEM simulation of pH at the centerline of the cathode for I cathode = 0.31 A/m 2 .F I G U R E 12 (a) FEM simulation of current steps after a short time and 1 day, and (b) schematic of the effect of oxygen reduction reaction limiting current density and anodic Tafel slope modification.[Color figure can be viewed at wileyonlinelibrary.com]