Numerical simulation of carbon steel atmospheric corrosion under varying electrolyte-film thickness and corrosion product porosity

A finite element model is developed to study dynamics of atmospheric corrosion of carbon steel, focusing on the influence of thin electrolyte film thickness under varying corrosion product porosity. Calculations have been done to evaluate the impact of electrolyte film thickness and corrosion product porosity on oxygen diffusion path, and the hindrance effect of corrosion products on the metal surface activity. The time evolution of corrosion current density and controlling steps in the corrosion process are explored. When the corrosion products are loose, oxygen diffusion is the dominant controlling step, and the thicker the electrolyte film, the lower the corrosion rate. When they are dense, the corrosion process is controlled by the mixture of oxygen diffusion and the surface discharge. The oxygen diffusion path is determined only by the corrosion product porosity, and therefore the corrosion rate is not affected by the electrolyte film thickness.


INTRODUCTION
Atmospheric corrosion is generally considered to be slow and mild. But it occurs widely in metal materials exposed to the atmosphere, which generates enormous economic losses yearly and may even cause catastrophic accidents 1 . Thus it is of great significance to explore the dynamics of atmospheric corrosion and establish models for precise prediction of its behavior.
According to the electrochemical theory of metal corrosion 2 , atmospheric corrosion is essentially an electrochemical process under a thin electrolyte film, including the dissolution of the metal and oxygen reduction. The thickness of the film affects the oxygen diffusion rate and further impacts the oxygen reduction reaction rate on the metal/electrolyte interface. The porous corrosion products not only partly cut down the electrode surface activity of metals but also increases the paths for oxygen diffusion 3 . And they can change significantly with the composition of the materials 4 and the serving environment 5 . Therefore, it is particularly crucial to study the influence of electrolyte film thickness and corrosion product porosity on atmospheric corrosion.
Concerning the influence of electrolyte film thickness, Atsushi Nishikata 6-9 estimated the corrosion rates of metals under different film thicknesses by measuring the cathodic polarization curves and electrochemical impedance spectroscopy (EIS). The conclusion can be made that the corrosion rate is independent of the film thickness when the film thickness is greater than 500-1000 μm. And when the film thickness is less than 500 μm, the corrosion rate negatively correlates to the thickness. While the corrosion rate either remains constant or decreases with decreasing film thickness when the film thickness is below a specific low limit. Dolgikh, O. 10 summarized the low limit of the measurement varies in different publications, ranging from 1 to 50 μm, which is related to the properties of the metal, electrolyte, test accuracy, and corrosion products. However, these electrochemical tests did not consider the effects of corrosion products. In fact, the dense rust layer can inhibit the anodic process and hinder the diffusion of dissolved oxygen. The denser and more stable the rust layer, the lower the corrosion rate. These conclusions have been confirmed by various experimental methods, including outdoor atmospheric exposure tests [11][12][13] , indoor accelerated tests [14][15][16] , and electrochemical and surface analysis combined reserches [17][18][19][20] . Therefore, to study the effect of electrolyte film thickness on atmospheric corrosion, it is necessary to consider the roles of corrosion products simultaneously. In addition, traditional experimental methods are challenging to quantify the influence of electrolyte film thickness and corrosion product porosity on the corrosion process owing to the evaporation too easily of the thin film. And it is difficult to monitor the control steps in the corrosion process and the evolution of ion distribution in the electrolyte domain.
In recent years, numerical simulation has been proposed as a method to link corrosion measurement results and corrosion process quantification, and then models for precise prediction of the atmospheric corrosion process have been established [21][22][23][24][25] . Venkatraman et al. 26 published a one-dimensional mathematical model for the natural corrosion of the bare metal surface under a thin electrolyte layer based on the mixed potential theory, indicating that the corrosion current and potential are functions of the electrolyte thickness. Thébault et al. 27 presented a twodimensional zinc/iron galvanic atmospheric corrosion model to illustrate the influence of electrolyte layer thickness and area ratio, taking electrolyte conductivity and oxygen diffusion. Furthermore, Hans Simillion 28 developed a mechanistic multi-ion transport and reaction model (MITReM), which considered the influence of electrode geometry and electrolyte film thickness, predicting the local concentration of ions and electrolyte potential. However, most of these models were for galvanic corrosion and didn't consider the influence of corrosion product deposition on the subsequent corrosion process. For the roles of corrosion products, San-Qiang Shi 29 established a model to demonstrate the pitting corrosion model under the deposited corrosion products by using the phase field method. In our series of the previous publications [30][31][32][33][34] , we step by step optimized a model of the micro galvanic corrosion dynamics triggered by the intermetallic particles embedded in aluminum alloy to investigate the coverage effect of corrosion products on electrode surface activity, the chemical effect on interfacial reaction due to the composition and concentration variations during corrosion process, and the volume hindrance effect of corrosion products of different porosities on the surface and also in the liquid phase. In addition, the repassivation due to the deposition of chromium hydroxide during pitting corrosion of stainless steel induced by MnS was explored 35 . But these models were established for immersed corrosion and didn't cover the corrosion process of carbon steel under a thin electrolyte film. Compared with the immersion corrosion model, the atmospheric corrosion model is very different, including the mass transport of species at the airelectrolyte interface, the limiting diffusion current density, etc.
Besides, atmospheric corrosion is a highly complex process. It is affected by many climatic factors, including temperature, relative humidity, rain, pollutant gas, and salt deposition [36][37][38] . The formation of corrosion products during the evaporation and condensation of the electrolyte film is very complicated 3,39 . Thus, it is extremely challenging to establish an accurate atmospheric corrosion prediction model. As the first step, we build a model describing the atmospheric corrosion of carbon steel under the thin electrolyte film. The model uses the measured polarization curves as boundary conditions, considering the mass transport of oxygen and other species, dynamic deposition of corrosion products, and corroded surface movement. In this paper, we focus on the modeling work to investigate the influences of electrolyte film thickness and the corrosion product porosity on atmospheric corrosion by quantitatively calculating the impact on oxygen diffusion and electrode surface reactions. Meanwhile, the time evolution of corrosion current density and control steps in the corrosion process is explored. This model can help develop accurate models for atmospheric corrosion prediction and provide insights into the kinetics of atmospheric corrosion under thin electrolyte films.

RESULTS AND DISCUSSION Model development
As shown in Fig. 1, the computational domain is a thin electrolyte layer above the metal surface. The thickness of the electrolyte layer (δ) ranking from low to high is 5 μm, 10 μm, 25 μm, 50 μm, and 100 μm, respectively. The mesh is only divided laterally due to the uniformity of the electrode surface to reduce calculation. As for the bottom boundary in contact with the substrate, the electrode reactions of iron dissolution and oxygen reduction, and for the top boundary, the electrolyte/air interface exchange of oxygen is considered. The liquid volume of the entire electrolyte domain remains unchanged during the corrosion process, thus movement of the corroded moved surface and the occupancy of the corrosion products can influence the instant position of the electrolyte/air interface. In this model, the initial electrolyte concentration of NaCl is set to 1 mol m −3 , with an initial pH value of 7. The initial oxygen concentration is the oxygen saturation value of 1 mol m −3 NaCl, about 0.18 mol m −3 .
In addition, the corrosion product porosity (ε) is introduced to quantitatively describe the compactness of the finally formed corrosion products. The porosity(ε) can be calculated according to Eq. (1): Here, C cp is the molar concentration of the corrosion products, C compact is the molar concentration of the compact solid FeðOHÞ 2 , the ρ FeðOHÞ 2 and M FeðOHÞ 2 are the density and molar mass of the compact solid FeðOHÞ 2 , respectively. Thus, the denser the corrosion products, the smaller the porosity. The critical minimum porosity (ε c ) indicates the most compact state the corrosion product can reach. In other words, the porosity (ε) of the corrosion products in the model is continuously changing, and when a critical minimum porosity (ε c ) is reached, the minimum compactness no longer changes. In this model, ε c is preset to be 0.01 and 0.4, representing compact and loose corrosion products, respectively. The accumulation of corrosion products can reduce the active sites on the electrode surface, which is described by the surface coverage(θ) and calculated by 1-ε in the first layer of mesh near the electrode surface. The oxygen reduction limiting diffusion current density is critical in atmospheric corrosion and is usually calculated as where F is the Faraday constant, D O2 is the coefficient of oxygen, c O2sat is the oxygen saturation concentration, and δ is the electrolyte film thickness. This calculation is reasonable without corrosion products. But when corrosion products are formed, the limiting diffusion current density value is a time-dependent transient value due to the inhibition of oxygen diffusion by corrosion products, calculated by Eq. (3).
where ðJ lim Þ z¼h represents the limiting diffusion flux through the first layer of mesh near the electrode surface.
Corrosion product porosity distribution Figure 2 shows the time evolution of corrosion product porosity of carbon steel along the longitudinal distribution in the electrolyte with different thicknesses. The critical porosities of the corrosion products ε c are equal to 0.01(a) and 0.4(b), respectively. The x-axis and y-axis represent the time and longitudinal distribution of the electrolyte, respectively. The color scale exhibits the porosity of the corrosion products. In both compact and loose cases, corrosion products deposition occurs in the vicinity of the electrode surface, gradually accumulating over time and reaching the densest state. In both ε c cases, the thicker the electrolyte film, the longer it takes for the corrosion products deposition to reach the densest state. This indicates that the thicker the film, the slower the initial deposition rate of corrosion products. In addition, it can be distinctly seen from Fig. 2 that the amount of the deposited corrosion products when ε c = 0.01 is less than that when ε c = 0.4. In other words, the proportion of positions containing loose depositions in the liquid phase is higher than that of the compact ones. Especially when the electrolyte thickness is relatively small, the loose corrosion products gradually cover the entire liquid phase over time, while the compact corrosion products are still mainly concentrated near the electrode surface. This difference leads to an apparent difference in the steric hindrance effect of two corrosion products on mass transport. Figure 3 presents corrosion product porosity distribution above the carbon steel immersed in different thickness electrolytes of 5, 10, 25, 50, and 100 μm for 12 h. The critical porosities of the corrosion products ε c are also 0.01 and 0.4, respectively. It is seen that, after 12 h, the electrode surface moved downward, accompanied by the accumulation of corrosion products. By comparing the two ε c cases, it is found that the corrosion is more serious with loose corrosion products, and the deposited corrosion product layer is thicker, which is consistent with the result shown in Fig. 2. When the corrosion products are compact (ε c = 0.01), the corrosion depth and the cumulative amount of corrosion products are basically the same at different electrolyte film thicknesses. As the electrochemical activity of corrosion products is ignored, it reflects a situation close to passivation. But unfortunately, the corrosion products of carbon steel are usually quite loose, for an instance ε c = 0.4. Under this situation, a porosity gradient in corrosion products is gradually generated along the longitudinal direction, with the most compact corrosion products located on the surface of the carbon steel, and a looser layer in the electrolyte domain. As the thickness of the electrolyte increases, the corroded volume and corrosion products deposition decrease, which implies that the electrode reaction rate reduces with the increase of electrolyte thickness.

O 2 concentration and pH distribution
Oxygen diffusion plays a significant role in the atmospheric corrosion of carbon steel. In order to quantitatively describe the influence of electrolyte thickness on oxygen concentration under different corrosion product porosities, Fig. 4 shows the time evolution of O 2 concentration along the longitudinal distribution under electrolyte film with different thicknesses of 5, 10, 25, 50, and 100 μm. The ε c of corrosion products is 0.01 (a) and 0.4 (b), respectively. It is evident that, with the deposition and gradual accumulation of corrosion products, the distribution of O 2 concentration varies significantly with the two corrosion product porosities. In all cases, there is an O 2 concentration gradient in the electrolyte domain at the initial stage of corrosion. The thicker the electrolyte film, the lower the O 2 concentration near the electrode surface, and the wider range of the O 2 concentration in the electrolyte domain. With the corrosion products covering the surface, with ε c = 0.01, the O 2 concentration in the electrolyte increases gradually, and the gradient disappears and reaches stable. In the other case of ε c = 0.4, when the electrolyte film thickness is 100 μm, the electrolyte/air interface has the same oxygen concentration as the initial state, while the O 2 concentration on the electrode surface stabilizes at zero, suggesting a fully O 2 -diffusion-controlled electrode kinetics. When the thickness is 5, 10, 25, and 50 μm, there is still a decrease in electrolyte/air interfacial oxygen concentration, and the gradient exists in the whole electrolyte film. However, it is noted that all the O 2 concentrations near the electrode surface under these electrolyte films are not zero, whose value declines as the increase of electrolyte film thickness. In the next section, the control steps in the atmospheric corrosion process of carbon steel are discussed in detail in combination with the current density and coverage of the electrode surface.
The pH value distribution is another chemical environment of concern in the atmospheric corrosion process, greatly influencing the formation of corrosion products. The pH value of the electrolyte domain is mainly affected by the OH − ions generated by oxygen reduction reaction and the H + ions generated by hydrolysis of metal ions. In our calculation case, the pH value in the electrolyte domain is kept constant after a rapid increase from the initial value of 7 to 10.9. The pH value is slightly higher near the electrode surface, but the gradient of pH longitudinal distribution is very small. The difference between the minimum and maximum values is only 0.02. This is because the electrolyte film is thin (5-100 μm) and has no ion exchange with the outside environment in this model. A large number of OH − ions generated by the oxygen reduction reaction on the surface of the electrode immediately fill the entire electrolyte domain except for combining with metal ions to form hydroxides. So the overall pH value can still be raised and gradually reach equilibrium with the hydroxides. In addition, as the electrolyte film thickness increases, the rate of pH rise decreases. But this difference is negligible for pH changes throughout the corrosion process because this rises too fast. Similarly, different porosity affects the final stable pH, but the effect of corrosion products hindering mass transport is insignificant due to too much OH − generated in the initial process. A detailed introduction is in Supplementary Fig. 2. The calculation results show that the pH value of the electrolyte film in the neutral atmospheric corrosion of carbon steel is maintained in an alkaline environment independent of the thickness and the porosity of the corrosion products. O. Dolgikh 10 measured the pH value distribution in an electrolyte with a thickness of 100 μm, showing that the pH near the working electrode is 11.5, which is coordinated with our simulation results.
Some researchers have inferred the evolution process of carbon steel atmospheric corrosion products through outdoor atmospheric exposure tests and X-ray diffraction or Raman testing methods. That is, Fe(OH) 2 is first transformed into lepidocrocite (γ-FeOOH), and then evolved from γ-FeOOH to goethite (α-FeOOH) or hematite (Fe 2 O 3 ), and converted into magnetite (Fe 3 O 4 ) in the inner layer 3,4,20,39,40 . The simulation results support the evolution process of carbon steel atmospheric corrosion products obtained by the experiments from the perspective of the oxygen and pH concentration distribution in the electrolyte film. After the deposition of Fe(OH) 2 predicted in our model, the high pH and low oxygen concentration in the initial stage of corrosion make it easier for Fe(OH) 2 to be converted to γ-FeOOH. This is because γ-FeOOH is an amorphous iron oxyhydroxide in the early stage, and its thermodynamics is the most unstable. As the corrosion progresses, the oxygen concentration increases and the outer layer corrosion products are progressively converted to hematite Fe 2 O 3 andα-FeOOH. In contrast, the inner layer corrosion products  are more likely to generate Fe 3 O 4 because of the lower oxygen concentration on the electrode surface. In addition, the looser the Fe(OH) 2 is formed, the less likely it is to be converted to a more stable corrosion product. The above hypothesis will be included in our future model to verify its rationality.
The kinetic process of thin electrolyte film corrosion of carbon steel To analyze the kinetic process of carbon steel corrosion under the thin electrolyte film, Fig. 5a-c present the time evolution of electrode surface O 2 concentration (a-1, b-1, c-1), coverage (a-2, b-  Fig. 5a, in the situation of ε c = 1, the surface coverage of corrosion products is always zero, the entire electrode surface keeps active, and the current density for oxygen reduction and metal corrosion is only related to the oxygen concentration. Oxygen is consumed by the oxygen reduction reaction and supplemented through oxygen transport from the air-electrolyte interface, as the result, the O 2 concentration drops sharply from the saturation value to a steady-state diffusion value. As the thickness of the electrolyte film increases, the stable value of the oxygen concentration on the electrode surface gradually decreases due to the rise in the oxygen diffusion distance. Similarly, the current density also declines, which is 0.67, 0.58, 0.38, 0.23, and 0.12 mA cm −2 under 5, 10, 25, 50, and 100 μm, respectively. This implies the corrosion rate is negatively correlated with the thickness of the electrolyte film. As can be seen from Fig. 5b, c, when corrosion products are deposited on the electrode surface, in all cases, the electrode surface coverage increases from 0 until it reaches the critical minimum porosity (0.4 or 0.01 Fig. in Fig. 5b, c, respectively). In comparison to the case of thinner electrolyte film, it takes a longer time for the thicker one to reach the most compact state, i.e., ε = ε c on the surface. This is because the initial current density decreases with increasing electrolyte thickness, which further leads to a decrease in corrosion product accumulation rate. On the other hand, with the same electrolyte thickness, more compact corrosion products results in longer time to reach ε c .
For case ε c = 0.01 (Fig. 5c), the current density gradually decreases until it reaches a minimum during the deposition of corrosion products. The thicker the electrolyte film, the slower the rate of decline. Particularly, when the electrolyte film thickness is 50 and 100 μm, a plateau appears at the initial stage. Correspondingly, consumption of O 2 decreases, and its concentration increases to a maximum. However, the oxygen diffusion path continues to be limited as the continuous accumulation of corrosion products. Thus the O 2 concentration on the electrode surface turns to decline again slowly. For case ε c = 0.4 (Fig. 5b), when the electrolyte film thickness is 5-25 μm, the current density declines with the coverage rising, and the O 2 concentration also increases and then decreases. The difference is that the current density continuously remains relatively steady throughout the corrosion process when the electrolyte film thickness is 50-100 μm. Correspondingly, the O 2 concentration at the electrode surface is always maintained very low.
In addition, the final stabilized current density is also negatively correlated with the electrolyte film thickness when the critical porosity is 0.4, which is 0.29, 0.28, 0.24, 0.18, and 0.11 mA cm −2 under 5, 10, 25, 50, and 100 μm, respectively. However, the final current densities are equal under different electrolyte film thicknesses at the critical porosity is 0.01, which is 0.08 mA cm −2 , indicating that the electrochemical corrosion process is only influenced by the surface protective corrosion products. Of course, please note, this situation is an ideal condition without taking the electrochemical activity of the corrosion products and its effect on the subsequent corrosion process of the matrix.
The control steps of atmospheric corrosion are mainly associated with the diffusion step and/or the surface discharge step. The electrode surface coverage by corrosion products can impede the surface discharge process by reducing the effective active area on the surface, while the occupation of corrosion products in the electrolyte can hinder either the leaving of products or the supply of reactant. In order to further explore the control steps in the corrosion process with different compactness of corrosion products and under different thickness of electrolyte film, combined with the prediction of oxygen concentration in evolution Figs. 5, 6 depicts the schematic diagram of the difference in control step, in which, the ε c of corrosion products is 1.0, 0.4, and 0.01, respectively. The pie chart represents the influence weight of oxygen diffusion and electrode surface discharge in each case, and the size of circle qualitatively represents the duration. The evolution of the weight of the control step can be obtained by calculating the ratio of the oxygen concentration to the saturation concentration on the metal surface over time. Among them, when the oxygen concentration on the metal surface is 0, the corrosion process is completely controlled by the oxygen diffusion process. When it is equal to the saturation concentration, the corrosion process is completely controlled by the electrode surface discharge process. Between the two represents controlled by a mixture of the two processes. The control steps and their weight of influences remain unchanged throughout the corrosion process in the absence of corrosion products. When the electrolyte film thickness is 5 μm, 10 μm, 25 μm, 50 μm, and 100 μm, the weight of oxygen diffusion control is respectively 74%, 85%, 97%, 100%, and 100%, with the rest proportion being surface discharge control. Those are also the control ratios for ε c = 0.4 and 0.01 cases prior to the deposition of corrosion products. When the porosity of corrosion products is 0.4, the proportion of electrode surface discharge increases with time, owing to the gradual deposition of corrosion products. The thinner the electrolyte film, the earlier the corrosion products depositing, and the greater rise in proportion of electrode surface discharge. When the thickness is 5 μm, 10 μm, 25 μm, and 50 μm, the influence weights of electrode surface discharge increase to 37%, 30%, 15%, and 5%. However, when the thickness is 100 μm, it is still completely controlled by oxygen diffusion. Interestingly, under the situation of the ε c = 0.01, because of the more substantial blocking effect from the corrosion products, the corrosion enters a brief period of 95% dominated by the electrode surface discharge. After that, due to the weakening of oxygen consumption, it gradually changes to a mixed control of the electrode surface discharge and oxygen diffusion, with the weight of 55 and 45%, respectively. It is worth noting that the final control state is independent of the electrolyte film thickness under dense corrosion products. A thicker electrolyte film only prolongs the time before the corrosion products is deposited, and then shortens the intermediate process to reach the steady state.
As we know, the current density is solved by the following equation 41 .
where f 1 ðE corr Þ ε and f 2 ðE corr Þ ε represent the electrode kinetics of anodic and cathodic reactions, respectively. I L is the oxygen reduction limit diffusion current density, determined by the corrosion product porosity affecting the diffusion coefficient of species in the corresponding locations, and the thickness of the electrolyte film. The thicker the electrolyte film, the denser the corrosion products and the longer the oxygen diffusion path, resulting in a lower limiting current density. Therefore, electrolyte film thickness and corrosion product porosity affect the final corrosion current density. Combining Figs. 5 and 6, as the corrosion product porosity decreases, the steric hindrance effect of corrosion products on oxygen diffusion becomes stronger, weakening the influence of electrolyte film thickness on oxygen diffusion. This influence decreases to 0 when the porosity is 0.01, resulting in the oxygen diffusion being independent of the thickness. This explains that although oxygen diffusion is still one of the controlling steps when densifying corrosion products, the corrosion current density is independent of film thickness. If the corrosion product is loose and the electrolyte film is thick, the effect of the loose corrosion product on the path of increased oxygen diffusion is much less than that of the electrolyte film thickness itself. This interprets why the current density can remain relatively constant under thicker electrolyte films at a critical corrosion product porosity of 0.4. Figure 7a further displays the average corrosion rate within 12 h, v, with different ε c and under the electrolyte film with various thicknesses, δ, where the blue ball represents the corrosion rate value. Figure 7b is the corresponding interpolation surface graph. The dashed line on the XZ side represents the regular pattern of average corrosion rate with the electrolyte thickness, i.e., v versus δ relationship, under different ε c conditions. It is shown that the average corrosion rate within 12 h decreases with the increasing film thickness for both ε c = 1 and ε c = 0.4 cases. Smaller critical porosity leads to a gentler decreasing trend. The average corrosion rate is independent of film thickness when the ε c is 0.01. Meanwhile, the solid line on the YZ side represents the rule of average corrosion rate with ε c , i.e., v versus ε c relationship, at different thickness of electrolyte. It is seen that at a certain electrolyte thickness, the average corrosion rate increases as the critical porosity of the corrosion products increases, and the thinner the electrolyte film, the more v depends on ε c .
Since the model simulates the atmospheric corrosion process under a constant salt film, it simplifies the transformation process of carbon steel atmospheric corrosion products, which may differ from the actual annual corrosion rate. However, the relationship between the electrolyte film thickness, the corrosion product porosity and the corrosion rate demonstrated and predicted by the simulation results is still valid and consistent with the reported experimental phenomena. For instance, it is believed that the corrosion rate decreases with time, which is assigned to the formation of corrosion products, but this decline is slow and complicated 13,42,43 . Through dry-wet cycle corrosion experiments on weathering steel, Junhua Dong 20,44 found that the corrosion rate drops significantly only when the formed corrosion products are compact. Similarly, our simulation results suggest that the denser the corrosion products, the more pronounced the downward trend of the corrosion rate with time. And the calculation results show that even with the same critical porosity, the looser corrosion products distribute in the electrolyte domain, while the denser ones are located next to the electrode surface, which is consistent with the experimental phenomena 3 . As mentioned in literature reviews [6][7][8][9][10] , the atmospheric corrosion rate increases with the decline of electrolyte film thickness when the film thickness is less than 500 μm and reaches a maximum value under 1-50 μm. The calculation results of the model likewise indicate that the maximum value is respectively below 5 μm or above 100 μm when the corrosion products are loose or compact. This implies that the role of corrosion products should be considered when exploring the electrolyte film thickness corresponding to the maximum rate of atmospheric corrosion. Moreover, the model uses critical minimum porosity (ε c ) to characterize the compactness of different rusts, and quantitative predictions can be made on the protective ability of various rusts. For instance, in the pure atmosphere (no pollutants) or in the atmosphere with low Cl − deposition, the main corrosion products formed on steel are lepidocrocite, goethite, and magnetite. And in the Cl − rich atmosphere, it is easy to form akaganeite. Akaganeite has more pores thus providing less protection to the matrix than lepidocrocite and goethite [45][46][47] .
Compared with experimental research, the advantage of numerical simulation lies in its capability of extracting the influence of any single factor individually in a coupled and complex system, as well as quantifying the influence weight of more factors on the corrosion kinetics under various situations. Of course, the atmospheric corrosion of carbon steel is a very complex system. The current modeling work is only our first step to simulating atmospheric corrosion basing our previous experiences on immersion corrosions 33,34 . The most important contribution of this work is to introduce the deposition and subsequent impact of corrosion products with different porosities into the modeling of corrosion under thin electrolyte film. Moreover, also for the first time, the limiting diffusion current density depending on situations, a characteristic related to atmospheric corrosion, is introduced to the model. In our next work, more important processes will be melted in the model, such as the electrochemical activity of the corrosion products, etc.
In summary, we have developed a FEM model for simulation of atmospheric corrosion of carbon steel under thin electrolyte film, considering dynamic deposition process of corrosion products, especially the effect of electrolyte film thickness and the porosity of corrosion products. Based on the simulation results, the following conclusions can be drawn: The transient corrosion rate gradually decreases with the deposition of corrosion products over time, and the thicker the electrolyte film, the slower the decreasing trend. Especially when the ε c is 0.4 and the electrolyte film thickness is 50-100 μm, the corrosion rate remains relatively stable throughout the corrosion process.
The average corrosion rate over 12 h is positively correlated to the critical porosity (ε c ) of the corrosion products. It decreases with increasing electrolyte film thickness when ε c is high (>0.4), but becomes independent of electrolyte film thickness when ε c is 0.01.
The change of corrosion rate is closely related to the controlling step in the corrosion process, namely oxygen diffusion and metal surface discharge. Corrosion products not only hinder mass transport but also impede the surface discharge of the metal. For thick electrolyte film or loose corrosion products oxygen diffusion control is dominant, whereas for dense corrosion products, surface discharge control is dominant.
The corrosion products have a gradient structure: i.e., denser near the electrode surface and looser in the electrolyte domain, with a pH of about 10.9 in all cases. O 2 in the entire electrolyte film is close to saturation when the corrosion products are dense, whereas it is depleted near the electrode surface when the corrosion products are loose.

Governing equations
The governing equation of the electrolyte domain must satisfy three conditions, i.e., the conservation of mass, principle of electrical neutrality and Laplace's equation. The conservation of mass is described by the Nernst-Planck equation, which mainly includes mass transport and chemical reactions.
The Nernst-Planck equation is modified through a Bruggeman effective coefficient of porosity (ε) as follows: where t is time, F is the Faraday constant. D i , Z i , C i , the diffusion coefficient, charge number and concentration for species i. ∇ϕ l is the potential gradient in the electrolyte and m i is the electric mobility for species i. R i expresses the reaction rate of the chemical reaction of ions in the electrolyte. The specific parameters can be found in our previous publications 34 .
In addition, the principle of electrical neutrality is expressed as the number of charges in the electrolyte is zero.
The effective production rate of species (R i;eff ) is introduced as shown in Eq. (3).
(8) where k forward and k backward are forward and backward reaction rate constant, and ζ and ζ are the corresponding reaction order. r, θ and z are position coordinates, and t is the time. The homogeneous reactions involved in the model and reaction rate are listed in Table 1.

Boundary conditions
The top boundary condition describes the exchange flux of the mass at the air/electrolyte interface. In the model, only the transport flux of oxygen at the interface is considered. The oxygen flux entering the electrolyte is related to the ratio of the local concentration and the saturated concentration, as shown in the following equation: where F O2 is the oxygen flux at the air/electrolyte interface, F maxO2 is the limiting value for this flux and C b O2 is the oxygen solubility in water at a temperature of 25°C. In Table 1 are given the constant values.
The bottom boundary is the electrode surface where local electrochemical reactions occur. The cathodic reaction is the reduction of oxygen and the anodic reaction is the dissolution of iron. The polarization curve data expressing the relationship between the electrochemical reaction potential and current is collected from the carbon steel rotating disk electrode (RDE). The polarization curve test adopts a three-electrode test system, in which the reference electrode is a saturated calomel electrode (SCE), the counter electrode is a platinum sheet, and the working electrode is carbon steel. The test range is −0.5-0.5 v relative to the open circuit potential, the scanning rate is 0.33 V/s, and the rotation speed is 2000 rpm. The final electrochemical boundary conditions are obtained by fitting the polarization curve based on the Tafel relationship. It is worth mentioning that the RDE test provides a uniform transport condition, from which can derive the key electrochemical kinetic parameters of the thin electrolyte layer. The experimentally tested polarization curves and fitted curves shown are in Supplementary Fig. 1.
The electrochemical reactions include the following equations, -Cathodic reaction: -Electrochemical dissolution of the iron: Fe ! Fe 2þ þ 2e À These electrochemical reactions satisfy exponential relations between the current density and the electrode potential, as seen in Eq. (12) and Eq. (13), respectively: where i 0;Fe and i O2 are the exchange current density, E eq;Fe and E eq;O2 are the Equilibrium potential, A Fe and A O2 are the Tafel slope. These parameters are fitted by testing the polarization curve of the RDE, as shown in Table 2.

Calculation platform
The partial differential transient equation of the control model is solved by the FEM method in Comsol Multiphysics R . The movement of the boundary in the model is tracked by an Arbitrary Lagrangian Eulerian method. The mesh size used in the computational domain is small enough to ensure that further reduction of the mesh size does not influence the calculation results.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.