Modeling localized corrosion in biofuel storage tanks

1.   Introduction

The global transition from traditional fossil fuels to biofuels marks a significant step toward reducing greenhouse gas emissions and promoting sustainable energy practices. Biofuels, which are derived from renewable biological resources, offer an environmentally friendly alternative to conventional petroleum-based fuels ^[1,2]. However, this shift introduces new challenges, particularly regarding the storage and handling of biofuels. Atmospheric Storage Tanks (ASTs), typically made by carbon steel ^[3], are essential infrastructure for the biofuel industry. However, they are susceptible to corrosion, which is a major concern for their long-term integrity, as it can lead to the release of hazardous materials causing, in turn, serious damages to the environment and living beings.

Corrosion is the deterioration of a material due to a chemical or electrochemical interaction with its environment ^[4]. Localized corrosion is a specific type of corrosion that occurs at discrete sites on a metal surface rather than uniformly. This can result in the formation of pits, crevices, or other defects ^[5]. Pitting corrosion, which is a form of localized corrosion resulting in the formation of small holes ^[6], is especially prevalent in tank bottoms. One of the challenges that green industries face is the increased susceptibility of ASTs to localized corrosion when storing biofuels ^[7]. The bottom of an AST is particularly vulnerable to these corrosion processes due to several factors: water accumulation from product contamination, condensation, or rainwater infiltration; the presence of contaminants within the stored product; and stagnant conditions in certain areas of the tank bottom that speed up the process of localized corrosion ^[8]. In particular, water can dissolve various substances present in the fuel, such as low molecular weight carboxylic acids, carbon dioxide, and chlorides, which are derived from biomass decomposition and other processes ^[9]. Unlike traditional fossil fuels, biofuels introduce unique challenges due to the presence of corrosion-causing byproducts such as organic acids, aqueous carbonic acids, and chlorides (depending on the feedstock). The term "organic acid" generally refers to a wide variety of organic compounds that contain the carboxyl (COOH) functional group, which plays a crucial role in the corrosion mechanism ^[9]. Biofuels, such as Fatty Acid Methyl Esters (FAME) and Palm Oil Mill Effluent (POME), are primarily composed of organic fats ^[10,11,12]. Hydrolysis, oxidation, and microbial activity are the primary mechanisms contributing to the deterioration of FAME and POME. These degradation processes result in the formation of free fatty acids, glycerol, and other corrosive compounds, thus leading to an increase in the total acid number of the biofuel ^[7,12]. Among the major corrodents identified in biofuel streams, there are low molecular weight organic acids, such as formic acid, acetic acid, and propionic acid, as well as carbon dioxide and chlorides. These stem from biomass decomposition, mineral acids used for hydrolysis, and chlorides found in algae grown in salt water. While some of these chemicals can be removed before further processing, others are either generated at high temperatures or remain uneconomical to eliminate. This issue necessitates the use of corrosion-resistant materials in biofuel processing units ^[9]. The presence of these corrosive compounds significantly contributes to the development of localized corrosion within ASTs. However, it is crucial to clarify that carbon steel does not exhibit passivity in neutral or low-pH environments ^[13]. Pitting corrosion in carbon steel ASTs is primarily attributed to material heterogeneities, weld-induced defects, and non-uniform water distribution rather than a classical pitting corrosion mechanism observed in passivated materials ^[14]. The presence of organic acids in biofuels can contribute to the overall corrosion rates by reducing the pH of accumulated water; however, the primary mechanism of localized attack remains associated with material imperfections and environmental factors.

To ensure the integrity of ASTs and prevent potential leaks, it is essential to regularly monitor the thickness of their bottom sections. This typically involves comprehensive inspections that are conducted every 10 years or more ^[15,16]. However, discrete thickness measurements cannot exactly determine the maximum corrosion depth of the bottom of storage tanks, where materials may exhibit localized corrosion in the form of pits. To assure the reliability of a given material under a certain environmental condition, quantitative predictions of the local corrosion rate as a function of exposure time and under various environmental influences should be made, which must include the most significant chemical reactions therein occurring.

This work constitutes an attempt to fill this gap, as it aims to develop a predictive mathematical tool for the estimation of corrosion rates in the bottom of ASTs. To achieve this purpose, the model needs to track the time evolution of the interface between the metallic tank and the electrolyte solution, which is generated from the deposition of water in the bottom of the tank. In building such a framework, it becomes crucial to identify those underlying mechanisms which, once the pit initiation stage has been concluded, rule the spatio-temporal propagation of such corrosion pits.

A consolidated framework, named the so-called phase-field (PF), has been used to account for continuously changing complex interfaces and to track the pit growth and the mutual interaction that arises when more localized corrosion sites form and grow in close proximity ^[17]. This method, commonly known as the "diffuse interface" model, is a mathematical framework that is employed to study interfacial phenomena ^[18]. The PF description has been employed to describe the evolution of an interface by means of an auxiliary order parameter that undergoes a smooth variation between two phases (metal and electrolyte, in our context) and whose dynamics is ruled by an Allen-Cahn balance equation. This approach avoids the complexities associated with the need of explicitly tracking moving boundaries.

Additionally, modeling the corrosion in the bottom of fuel tanks needs to be complemented by a set of Nernst-Planck equations to describe the transport of ionic species balance of metal ions ^[17,19,20] and by a Poisson equation which rules the spatial distribution of the electrostatic potential.

In this paper, a relatively simple PF model is employed to qualitatively capture the most relevant electro-chemical reactions that occur at the metal/electrolyte interface in the bottom of a biofuel tank. Indeed, to keep the complexity of the proposed model at a reasonable level, some additional effects, such as those associated with stress corrosion cracking, mechanically- and microbiologically-assisted corrosion, are disregarded here.

The novelties here are introduced are manifold, as they encompass the possibility to numerically investigate how the electrolyte solution is affected by the presence of a biofuel, which additional chemical reactions can be responsible for the larger corrosion rates typically observed in experiments on biofuels in comparison with traditional fuels, how the pit growth rate can be modeled, which dynamics among multiple interacting pits takes place, which characteristic stages the pit undergoes during its growth and which model parameters govern such dynamics.

The paper is organized as follows: in Section 2, the generalized electro-chemical PF model of pitting corrosion, which includes the additional chemical reaction due to the presence of organic acid in the biofuel, is presented; in Section 3, results of numerical investigations are addressed to characterize the time evolution of the pit depth and the corrosion rate observed in configurations with single and multiple interacting pits; and in Section 4, concluding remarks and some future research directions are provided.

2.   Materials and methods

To investigate the evolution of the metal/electrolyte interface and predict the corrosion rate caused by single and multiple interacting localized corrosion areas on carbon steel, a PF-based numerical framework is presented, which will be named as pits for simplicity.

Unlike sharp interface models, which require tracking a free boundary, the PF method introduces an auxiliary field variable, φ, defined across the entire domain Ω. This variable undergoes smoothly transitions between the two physical states at the interface, mimicking the two different phases (solid and liquid). In the model presented here, φ = 1 corresponds to the solid phase (metal) while φ = 0 corresponds to the liquid phase (NaCl-based electrolyte solution). A smooth variation of φ occurs within the interfacial region, where the pitting corrosion front evolves. Consequently, the free boundary problem is formulated as a partial differential equation of the Allen-Cahn type, which describes the spatiotemporal evolution of the PF variable. This method allows us to effectively manage complex geometries and to handle topological changes ^[18,19].

To describe the localized corrosion process in detail, let us consider a 2D domain composed of a metal phase, such as carbon steel, which is exposed to a NaCl solution simulating water solution that can be found on the AST bottom ^[21]. A schematic of the electrolyte/metal interface, together with some pictures of real biofuel tank bottoms that reveals the presence of circular pits, is shown in Figure 1. Once the dissolution process starts, it leads to the release of metal cations into the electrolyte. The general form of the anodic reaction is represented by the following:

where Fe represents the corroding metal, and e^- refers to electrons released during the process. To reduce the model complexity, it is assumed that the corrosion process solely takes place within the pit area. To this aim, outside the pit area, an inactive layer is introduced to separate the solid phase from the liquid one.

To model the evolution of corrosion waves at the interface between the electrolyte and metal, it is necessary to specifically consider the presence of the following anionic and cationic species: Fe²⁺, FeOH⁺, Na⁺, Cl^-, H⁺, and OH^-. Moreover, additional chemicals need to be included because of the extra reactions associated with the specific nature of biofuels that are mainly composed of methyl esters produced through the esterification of oils. When esters interact with water, they break down into organic acids and alcohol ^[22], thus leading to more organic acid formation. In turn, these acids dissociate in the aqueous environment, thus affecting the acidity of the mixture and lowering the pH.

The above additional chemical reactions will be specified in Section 2.2. In detail, the following subsections tackle different topics. Section 2.1 focuses on the theoretical framework, based upon the PF model, which is used to describe the localized corrosion process. Section 2.2 describes the transport of dilute ionic species, thus providing a detailed analysis of the key chemical reactions involved in the investigated phenomenon of biofuel-induced corrosion. Section 2.3 addresses the spatial distribution of electrostatic potential.

2.1.   A phase field description of pitting corrosion

To describe the localized corrosion that occurs in the bottom of biofuel tanks, it is necessary to generalize the framework used with success in the literature ^[17,19,20]. The time evolution of the PF variable φ is described by an Allen-Cahn equation ^[23], which reads as follows:

where L denotes the interface kinetics coefficient, Ψ represents the total electrochemical energy, and τ is the gradient energy coefficient. Moreover, the coefficient L can be correlated with the corrosion current density i ^[19,20], or the overpotential η ^[17] in the case of activation-controlled corrosion, through the following:

where i₀ represents the exchange current density, and L₀ represents the interface kinetics coefficient when the overpotential η vanishes. This latter parameter is defined as η = φ₀ -φ_l -E_eq, where φ₀ represents the applied (solid) potential, φ_l represents, the electrostatic potential, and E_eq represents the equilibrium potential (assumed to be zero, for simplicity) ^[17].

To estimate the corrosion current density i, the Butler-Volmer equation is used under the assumption that the only significant reaction at the dissolution interface is the anodic one:

in such a way that the corrosion current density i takes the form of a Tafel-type ^[21,24]. In (4), the parameter a represents the anodic charge transfer coefficient, n_M represents the charge number of the metal, F is the Faraday constant, R_g is a gas constant, and T is the temperature. According to the above assumptions, L becomes the following:

The total electrochemical energy Ψ is an integral of the corresponding energy density ψ over a 2D domain Ω, enclosing the electrolyte solution and the metal layers. The above energy gathers the following contributions:

where ψ^(α), ψ^(β), and ψ^(γ) represent the chemical, interfacial, and electrostatic free energy densities, respectively. The chemical free energy density can be considered as ψ^(α) = ψ^{(α, d)}+ψ^{(α, s)}, where ψ^{(α, d)} is the energy for metal dissolution and ψ^{(α, s)} is the energy stored in the dilute solution$.$ In line with the Kim-Kim-Suzuki (KKS) model proposed in ^[25], each material point is considered as a mixture of solid and liquid phases. Given that the two phases coexist with different concentrations, while maintaining equal chemical potentials, the energy for metal dissolution ψ^{(α, d)} can be expressed as follows:

where the degradation function is defined as h(φ) = -2φ³ + 3φ². This function is such that h(φ = 0) = 0 and h(φ = 1) = 1, which identify the fully corroded and the undamaged metal regions, respectively.

The other functions that appear, namely $\psi _S^{\left({\alpha, d} \right)} = A{\left({{c_S} - {c_{{S_e}}}} \right)^2}$ and $\psi _L^{\left({\alpha, d} \right)} = A{\left({{c_L} - {c_{{L_e}}}} \right)^2}$, represent the chemical free energy density terms associated with the concentrations of the solid phase, c_S, and the liquid phase, c_L, respectively. The quantities involved are as follows: ${c_{{S_e}}} = {c_{solid}}/{c_{solid}} = 1$, which represents the normalized equilibrium concentration for the solid phase; ${c_{{L_e}}} = {c_{sat}}/{c_{solid}}$, which represents the normalized equilibrium concentration for the liquid phase; A, which represents the free energy density parameter (assumed to be identical for both solid and liquid phases); ${c_{solid}}$, which represents the concentration of atoms in the metal; and ${c_{sat}}$, which represents the average saturation concentration. Moreover, additional constraints are given in line with the KKS model ^[25]. In detail, $\frac{{\partial \psi _S^{\left({\alpha, d} \right)}}}{{\partial {c_S}}} = \frac{{\partial \psi _L^{\left({\alpha, d} \right)}}}{{\partial {c_L}}}$ and, by introducing the normalized concentration of metal as ${c_M} = {{c_{_M}^0} \mathord{\left/ {\vphantom {{c_{_M}^0} {{c_{solid}}}}} \right. } {{c_{solid}}}}$, being $c_{_M}^0$ the ion metal concentration, and requiring that $c_{_M}^{} = 1$ inside the metal and $c_{_M}^{} = 0$ in the electrolyte solution (far from the evolving interface), it is possible to express the metal concentration as follows:

Finally, by combining all the above assumptions, the energy for metal dissolution can be described by the following equation:

On the other hand, the energy stored in the dilute solution is defined as follows:

where c_i represents the concentration of the i-th ionic species, and $\mu _i^0$ represents the reference chemical potential.

The interfacial energy density, ψ^(β), is given by the sum of double-well potential energy and the energy corresponding to the PF gradient:

where τ is the gradient energy coefficient, and ω is the height of the double-well potential g(φ) defined as g(φ) = φ²(1 − φ)². In the present PF corrosion model, these quantities are correlated to the interface energy per area σ and its thickness l via the following ^[17]:

The electrostatic energy density, ψ^(γ) , is expressed as a function of the charge density as follows:

where n_i is the charge number of the i-th species.

2.2.   Transport of diluted species

The transport of each aqueous chemical species in the dilute electrolytic solution, with a concentration c_i, is ruled by the mass conservation law and consists of three parts ^[21,24,26]: the migration of ions under concentration gradients, the migration under electrostatic potential gradients, and the rate of depletion or production by chemical reaction. It reads as follows:

where N_i and Π_i are the total flux in an electrolyte and the chemical reaction term of the i-th species, respectively. The total flux includes phenomena of diffusion (according to Fick's law) and migration (dictated by the electrolyte potential), whereas convection is neglected here, as usual in the literature ^[17,19,20].

In order to define the Nernst-Planck equations that govern the transport of ionic species, the saline solution that mimics the water collected on the bottom of AST is modeled here as a NaCl-based solution that contains the following six ionic species: c₁ ≡ c_M = [Fe²⁺], c₂ = [FeOH⁺], c₃ = [Na⁺], c₄ = [Cl^-], c₅ = [H⁺], and c₆ = [OH^-]. This approach holds true when a traditional fuel is stored into the AST. In this case, the most relevant chemical reactions therein occurring are as follows ^[17,21,26]:

where k_b and k_f represent the rate constants for the forward and backward reactions, respectively.

Alternatively, when AST is used to store biofuels, it is necessary to account for additional factors related to the occurrence of extra chemical reactions. In detail, biofuels primarily composed of methyl esters ${\text{R - COOR'}}$ produced through the esterification of oils, such as palm oil with methanol, can undergo hydrolysis in the presence of water (Ester Hydrolysis) ^[27]. Here, "R" denotes a straight or branched chain. Ester hydrolysis is a critical chemical reaction to consider when evaluating the stability and corrosive effects of biodiesel, especially in the context of its interaction with metallic materials ^[12]. This reaction breaks down the ester bonds, thus leading to the formation of organic acids ${\text{R - COOH}}$ and alcohols ${\text{R' - OH}}$:

where ${\text{R'}}$ represents an organyl group. This equation represents the reaction of an ester with water to produce an aliphatic (fatty) acid and an alcohol. In turn, these organic acids dissociate in the electrolyte, thus generating additional ions:

The dissociation of organic acids contributes to the overall acidity of the electrolyte by introducing additional H⁺ ions. This increased acidity can accelerate localized corrosion by promoting the dissolution of the metal at the anodic sites ^[22].

The above chemical reaction constitutes the most significant generalization introduced to describe the different corrosive behavior exhibited by biofuels in AST. This leads to include the presence of two additional species in our framework: c₇ = [R-COOH] and c₈ = [R-COO^-].

By including all the above considerations, and accounting for the concentrations of the mentioned eight ions, the overall Nernst-Planck equations (14) specialize into the following:

where D_M and D_i are the diffusion coefficients of the metal cation and of the i-th species, respectively. According to (15), (16), and (18), the reaction terms for each involved species, take the following form:

Notice that sodium Na⁺ and chloride Cl^- ions are assumed not to directly participate in chemical reactions (Π₃ = Π₄ = 0). However, they play a crucial role in facilitating the transfer of current within the pit by contributing to the overall ionic conductivity of the electrolyte ^[24].

2.3.   Electrostatic potential

The distribution of electrostatic potential φ_l can be estimated through a Poisson-type equation. This equation relates the local variations in the electrostatic potential to the distribution of electric charge within a specific region. It reads as follows ^[17]:

where θ is the electric conductivity of the electrolyte.

3.   Results

System (2), (19), (20), and (22) rules that the evolution of localized corrosion in the bottom of biofuel tanks is in the form of a Nernst-Planck-Poisson system ^[28], which is rewritten for the sake of clarity as follows:

where c_M = [Fe²⁺], c₂ = [FeOH⁺], c₃ = [Na⁺], c₄ = [Cl^-], c₅ = [H⁺], c₆ = [OH^-] and, due to the presence of the biofuel, c₇ = [R-COOH] and c₈ = [R-COO^-]. The reactions terms Π_M and Π_i which appear in (23.2) and (23.3) are given in (21).

The complexity of the above system prevents the possibility of addressing further analytical investigations. Therefore, to capture some aspects of qualitative behavior, it is necessary to rely on numerical investigations.

To highlight the main role played by biofuels, a distinction is made between two different frameworks: classical fuel and biofuel. In the former case, organic acids are absent and thus the concentrations of species c₇ = [R-COOH] and c₈ = [R-COO^-] are set to zero. In this case, eight partial differential equations (six balance laws for the species c₁, ..., c₆, together with those for the PF variable φ and the electrostatic potential φ_l) describe the key physical phenomena involved. In the latter one, organic acids are present in the electrolyte solution; therefore, it is necessary to account for the chemical reaction (18) and to include two further balance laws for the species R-COOH and R-COO^-, coupled with the previously-mentioned ones.

The set of PDEs is numerically integrated by means of finite-elements framework of COMSOL Multiphysics® ^[29]. To keep the computational cost at a reasonable level, a narrow region of interest was selected at the interface between a layer mimicking the NaCl electrolyte solution having a thickness H_l = 0.6 μm and another layer representing the top surface of the carbon steel fuel tank (the solid phase), having a thickness H_s = 0.4 μm (Figure 1).

Moreover, in order to allow the corrosion process to solely take place within the pit area, an "inactive" layer, having a thickness H_in = 1 µm, is introduced as a meshless layer separating the solid (metal) phase from the liquid (electrolyte) one (see Figure 1(c)). By means of this mathematical artifice, it is intended to describe a layer which does not affect any chemical reaction, does not undergo degradation and, thus, does not participate in the localized corrosion process.

It is also assumed that the pit initiation stage has already occurred ^[30]; therefore, the processes here solely numerically investigated the pit growth. During this stage, as known, the development of some of the initiated pits occurs, the rate of growth starts to increase and, possibly, the interaction among different pits takes place. In particular, pits with a radius r_pit = 10 µm constitute the initial condition for all the numerical simulations. This choice allows us to improve the convergence of the numerical algorithm and is in line with the literature ^[17,19].

To accurately assess the corrosion rates, the simulations are carried out over a time window of 30,000 sec and considers a computational domain having a width W = 0.6 μm.

Table 1 outlines the initial conditions for various field variables within the system. These conditions represent the state of the system at the beginning of the simulation. The variables include the PF variable (φ), the metal ion concentration c_M, various ionic species c_i, and the electric potential φ_l. The liquid and solid phases are differentiated, with specific values assigned to each phase for each variable. Table 2 details the boundary conditions for the PF variable, the metal ion concentration, and other ionic species, imposed at the top surface of the electrolyte layer (x = 1000 μm) and at the bottom of the metal layer (x = 0 μm). Moreover, zero flux boundary conditions are applied to the lateral boundaries. The additionally parameter values used in all simulations are reported in Table 3.

The numerical values reported in these tables are extracted from the literature ^[17,21,24]. The only exception is constituted by the values of the initial concentrations of R-COOH and R-COO^- reported in Table 1, which are not generally available in the literature. In the following, the procedure to deduce their values is outlined.

Let us introduce the concept of the partition coefficient κ, which quantifies the distribution of a solute between two immiscible phases ^[31]. In our case, κ will quantify the distribution of the organic acid between the aqueous electrolyte (where it dissociates) and the palmitic oil phase:

Moreover, as outlined in ^[17], the equilibrium constants K₁ associated with reaction (15) quantifies as follows:

Similarly, K₂ characterizes the equilibrium of reaction (16) and is determined by the following:

In the specific case of acid dissociation, where a weak acid forms its conjugate base, the equilibrium constant is given by K_a = [R-COO^-][H⁺]/[R-COOH], and is referred as the acid dissociation constant. On the other hand, pK_a is related to the acid dissociation constant through pK_a = -log(K_a) ^[32,33,34].

According to the reaction (18), the dissociation equilibrium constant K₃ ≡ K_a of the organic acid in the electrolyte solution is given by the following:

Therefore, taking the partition coefficient for acetic acid between oil and water given in ^[35] κ = 0.075 into account, and extracting the mass concentration of acetic acid in oil s_o through the knowledge of typical molar concentration and molar mass ^[36],

it is possible to deduce the following:

Then, let us consider a mixture containing 0.15% water and 99.85% oil. On the total amount, the weight percentage of acetic acid in water W_w and the weight percentage of acetic acid in oil W_o are given by W_w = 1.96% and W_o = 98.04%, respectively.

Therefore, the solubility of acetic acid in the mixture s_aa is determined by considering the solubility in both water and oil:

Once the solubility has been determined, considering the molar mass of acetic acid PM^{(acetic_acid)} = 60.052 g/mol, the concentration of the undissociated form [R-COOH] is calculated by the following:

Finally, considering that the pK_a value for acetic acid is 4.7, as reported in ^[37], and that the initial pH of the solution is 4.8 (see Table 1), the concentration of the dissociated form [R-COO^-] is given by the following:

Numerical investigations have been performed by following this criterion. First, the corrosion process is modeled under the simplistic assumption that only a single pit forms and evolves. In such an analysis, an attempt is made to extract a suitable predictive law and highlight the most relevant differences between classical fuels and biofuels. Then, the analysis solely focuses on the biofuel, and the most relevant features that occur when two or three pits interact and evolve on the metal surface in close proximity are described.

3.1.   Single-pit corrosion

It is expected that the most relevant physical quantities that contribute to the evolution of a single pit are as follows: time t and solid potential φ₀. To investigate this dependence in detail, several numerical simulations have been carried out by considering a computational time window t = [0, 3e4] sec, and by varying the solid potential in the range φ₀ = [0.3, 0.42] V. To accelerate observations on pit growths within a computationally feasible timeframe, the chosen solid potential values are intentionally larger than those that might be encountered in real-world applications.

Then, a nonlinear regression analysis on these numerical data is performed, for both classical and biofuels, in order to find the best fit for the function depth(t, φ₀), which describes the time evolution of the bottom of the pit (i.e., the worst situation). This analysis has led to the following expression:

where B₀, B₁, and B₂ are fitting parameters, whereas t^* and φ₀^* are normalization coefficients. In detail, B₀ is a characteristic depth which also determines the pit growth rate, and B₁ and B₂ describe the power laws for time and solid potential, respectively. The normalization coefficients are arbitrarily set in such a way that t^* = 31,536,000 sec, which corresponds to the period of 1 year, whereas φ₀^* = 0.42 Volt is the largest value considered for solid potential.

For classical fuel, the best-fit values of these parameters are as follows:

whereas, for biofuel,

Among the possible functional dependencies that can be considered to fit a real function of two-variables, it should be mentioned that this law appears to be attracting since it encloses just three fitting parameters and, above all, the power laws have been already used with success in the literature, especially to describe the time evolution of pit growth ^[30,38].

The numerical fits reported in Eq (33) are depicted in Figure 2 through surfaces for classical fuel (panel (a)) and biofuel (panel (b)). In the same figures, the results of numerical simulations are represented by solid lines. A regression analysis for both classical fuel and biofuel yielded an adjusted R-square value of 0.997, thus indicating a very strong fit of the model to the data, as it also appears from the close agreement reported in the figures. From the close agreement between the numerical and analytical data, it can be argued that the model in Eq (33) appears to be quite suitable to predict both short-time and long-time behavior (i.e., from the conclusion of the pit initiation stage to the achievement of quasi-steady-state).

Additionally, the results depicted in this figure clearly indicate that the pit depth observed in biofuels is always larger than that reported for classical fuels, which is in qualitative agreement with literature ^[22].

Moreover, from the inspection of the above fitting data, it is worth noticing that the parameter providing the most relevant contribution to the higher corrosion rate exhibited by biofuel is the prefactor B₀, which exhibits an increase of about 12% with respect to that of a classical fuel. Interestingly, all the other parameters exhibit negligible differences (smaller than 2%), so that the power laws describing time and potential dependencies may be ruled by the same function. Noteworthy, this enables the possibility to characterize the typology of a fuel solely by means of a single parameter, that is, the role played by additional chemical reactions due to the presence of organic acids in the electrolyte solution may be enclosed in the sole coefficient B₀.

In our view, the results discussed here represent one of the most relevant achievements of this work as the possibility to manage a best fit function allows us to address some quantitative predictions on the time required to achieve a given pit depth once the solid potential is known. This constitutes a more than helpful tool for green industries to optimize the protocol adopted to test the integrity of bottom tanks.

For instance, it is possible to address a quantitative comparison with experimental data reported by several groups ^{[12,22,27,39,40,41]}, where a quite wide range of corrosion rates due to biofuel in contact with carbon steel was reported, from 1 μm/year to about 2 mm/year. For instance, by considering the value of pit depth equal to 17 μm reported after 8 hours of laboratory observations (that roughly corresponds to our computational time window of 30.000 sec) (see Figure 28 in ^[40]), it is possible to get an indirect measure of the potential φ₀ = 0.3125 V by means of Eq (33). Then, by considering this value of applied potential, the same Eq (33) allows us to estimate the final pit depth after 1 year, which is 1.66 mm, which is again in line with the above literature.

In Figure 3(a), the time evolution of a single pit for the two fuels is compared for different values of solid potentials φ₀. From this figure, it can be seen that the larger corrosivity of biofuel increases with the rise in solid potential. Apart from the time evolution of the pit depth, it is interesting to evaluate its time derivative, which is representative of the corrosion rate and provides significant information on the different stages occurring during the corrosion process. Figure 3(b) shows the behavior obtained for a specific applied potential (φ₀ = 0.4 V), where the behavior associated to the classical fuel is given as the solid line while the one exhibited by biofuel is given via the dashed line.

Two key pieces of information might be extracted from this figure. First, the higher corrosivity of the biofuel manifests itself from the very beginning. Second, the time evolution of the rate of corrosion, in the case of a single pit, is a two-steps process characterized by an initial sudden increase up to a given maximum, followed by a slow decrease over time. This pattern, which appears to be independent of the fuel type, suggests that it might be a general and common feature of the time evolution of localized corrosion. It is also interesting to notice that this two-steps process was also observed in the literature of corrosion of metal in biodiesel ^[41] and even in other contexts (e.g., marine corrosion) ^[30]. The initial rapid growth is typically associated with favorable conditions (i.e., large supply of chemical agents) that support a high rate of corrosion. However, this rate reduces to a quasi-steady state situation that is supposed to continue indefinitely and whose precise underlying mechanisms are still unclear.

3.2.   Multiple-pit corrosion

From now on, the investigations will solely focus on the behavior of the biofuel. Additionally, for simplicity but without loss of generality, the analysis is limited to a solid potential of φ₀ = 0.4V. The behavior associated with other values of the potential is qualitatively unchanged.

In Figure 4, the time evolutions of the pit depth (panel (a)) and corrosion rate (panel (b)) are depicted for frameworks where a single, double, and triple pits are initially formed. In simulations with multiple pits, the initial condition consists of pits which are equal-in-size (radius r_pit = 10 μm) and located at a distance d = 20 µm from each other.

In the single pit scenario (dash-dotted line) already discussed in the previous figure, the maximum depth achieved at the final simulation time is approximately 180 µm. Notably, the presence of multiple pits significantly impacts the corrosion process. Indeed, in the double pit configuration (dashed line), the final depth is around 200 µm and reaches about 220 µm in the case of triple pits (solid line) (i.e., an increase of about 20% (according to this parameter set)). These observations strongly suggest that the interaction between adjacent pits exacerbates the corrosion process, thus leading to a more rapid increase in the pit depth. The analysis of the corrosion rate (panel (b)) revealed that the maximum corrosion rate obtained in the multiple pit configuration (double and triple), about 13.5 nm/s, is about 40% larger than that reported in the single pit case, 8.5 nm/s. Two key differences may be pointed out between the two multiple-pits configurations: in the triple one, the peak is slightly delayed and, after it is achieved, the decreasing trend of corrosion rate appears to be slower than that observed in the double pit.

A superficial analysis of Figure 4 might erroneously lead to the conclusion that the corrosion process in the presence of multiple pits is a two-stage process similar to the one observed in Figure 3(b). However, deeper investigations reveal that it is not. To prove it, let us consider a double-pit configuration with an initial spacing d = 20 μm and track the time evolution of the corrosion rate as well as the spatial profiles of the PF variable φ at each time instant, with a particular emphasis at those at which the corrosion rate exhibits a change in slope. The results of this analysis are reported in Figure 5.

In detail, the simulation starts with the initial (t = 0 s) configuration with two equal pits, see point P₁ in panel (b). As the simulation time progresses, at P₂ (t = 750 s), the corrosion rate reaches the value of about 10 nm/s and the curve exhibits a kink. At this point, the rate of corrosion has indeed risen. The reason for that is that a change of morphology is taking place: the boundaries of the two pit areas are now in contact, as shown in Figure 5(c). During the subsequent time instants, from P₂ to P₃, the initially distinct pits progressively merge together and form a unique pit with a larger volume, as can be appreciated in Figure 5(d), in line with both experiments ^[42] and other numerical simulations ^[19,43]. Here, the maximum corrosion rate of about 13.5 nm/s is obtained. Starting from P₃, the evolution is the one typical of a single pit.

To gain further insight into the process of corrosion due to multiple pits, let us now vary the initial spacing between them. In the case of double pits, d is considered in the range (20,160) μm, while for triple pits, d is applied in the range (10, 70) μm. These values are chosen to be compatible with the considered width of each layer and the simulation time window. Figure 6 reports the overall results (panel (a) for double, panel (b) for triple). In the same figure, small squares denote the time instant at which the boundary of the adjacent pit areas become in touch with each other. The results reveal that the closest pits get in touch earlier, as expected, and thus the corresponding pit depth exhibits a faster increase. However, these two combined observations don't imply that the pit depth achieves the largest value at the final simulation time. Indeed, since the corrosion rate starts monotonically decreasing once the maximum is reached, it happens that, after a substantial long time (in our case, the upper bound of the simulated time window), the pit depth observed for d = 20 μm appears to be lower than those observed for larger spacings. This behavior holds for double and triple pits.

The above considerations may be further supported and enriched by addressing an overall comparison among all the time evolutions of the pit depth for single, double, and triple pit configurations, for different values of initial spacings between the pits, which is in line with Figure 6. The results of this analysis are reported in Figure 7.

Here, the dashed line depicts the pit depth evolution for a single pit and serves as a baseline for the comparative analysis, as it represents the lower limit for the corrosion process. Dashed and solid lines denote the behavior for double and triple pits, respectively. The evidence that the increase of interacting sources enhances the pit depth is not the only conclusion suggested by this figure. In fact, these results reveal that the long-time behavior associated to a given number of interacting pits, independently of the considered initial spacings between pits, forms a "band", which appears disjointed with respect to the band formed by a different number of interacting pits. In other words, the long-time dynamics results to be mostly determined by the number of interacting pits rather than the initial spacings among pits. In our view, this is another intriguing feature to be considered in the design of the maintenance system of fuel tanks.

4.   Conclusions

AST bottoms are critical components susceptible to corrosion due to the accumulation of sediments and exposure to biofuels. This study employed a PF model to investigate the influence of biofuels on pitting corrosion at the bottom of ASTs made by carbon steel. To emphasize the role played by biofuel, the model presented here also accounted for the additional chemical reaction which described the dissociation of organic acids which are present in the electrolyte solution. Thus, the overall system included an equation for each of the eight ionic species, one for the PF variable and one for the electrostatic potential. The obtained nonlinear and coupled system was numerically integrated to describe the spatio-temporal evolution observed during the stage of pit growth. In particular, the dynamics of single and multiple interacting pits were inspected, and results agree with some experimental observations.

The most relevant achievements contained in this manuscript can be summarized as follows:

● A best-fit law describing the dependence of the pit depth as a function of time and solid potential was deduced. This function, apart from being characterized by just three free parameters (B₀, B₁, B₂), yielded the most significant advantage of highlighting that the presence of organic acids in the electrolyte solution may be enclosed in the sole coefficient B0. Thus, this fitting parameter can be used to identify and characterize the typology of biofuel.

● It is confirmed that biofuels exhibit a larger corrosivity in comparison with traditional fuels.

● As far as single, isolated, pits are concerned, the pit growth rate is a two-steps process.

● When multiple pits are formed in close proximity, their interaction yields a threefold contribution:

○ it enhances the pit growth;

○ the pit growth consists of three distinct stages; and

○ the time evolution of pit growth may be classified in terms of separated "bands", within which the curves representative of different initial spacings among pits fall.

These latter observations served to stress the importance of considering pit interactions during inspections and risk assessments for ASTs. At the same time, they suggested that while the proximity of multiple pits becomes crucial for the corrosion of the bottom of the tank, as it may potentially led to enhanced localized breaches and leaks, the number of such interacting pits plays a relevant role rather than their initial spacings.

It is believed that the findings described here may contribute to acquiring a better understanding on the influence of biofuels and pit interactions on corrosion behavior at the tank bottom. Additionally, they may help to develop more robust design strategies for AST bottoms, to optimize the inspection procedures to target areas susceptible to multiple interacting pits, and to select appropriate materials to ensure the long-term integrity and reliability of these critical storage systems.

Future research directions could include the following: modeling additional contributions arising from mechanically- and microbiologically-assisted corrosion and the inclusion of stochastic effects describing the initial formation of several, randomly dislocated, pits.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Authors gratefully acknowledge fruitful discussions with Prof. C. Cui and E. Martinez-Paneda. This research was funded by MUR (Italian Ministry of University and Research) PNRR - Missione 4, Componente 2, Investimento 1.1 – Bando Prin 2022, Decreto Direttoriale no. 104 del 02-02-2022, through PRIN2022 Project PRIN_202248TY47 ''Modelling complex biOlogical systeMs for biofuEl productioN and sTorAge: mathematics meets green industry (MOMENTA)" CUP J53D23003580006; by INdAM-GNFM and by INAIL within the BRIC/2021 ID = 3 project DRIVERS.

Conflict of interest

The authors declare there is no conflict of interest.