Unravelling the interplay between steel rebar corrosion rate and corrosion-induced cracking of reinforced concrete

Accelerated impressed current testing is the most common experimental method for assessing the susceptibility to corrosion-induced cracking, the most prominent challenge to the durability of reinforced concrete structures. Although it is well known that accelerated impressed current tests lead to slower propagation of cracks (with respect to corrosion penetration) than in natural conditions, which results in overestimations of the delamination/spalling time, the origins of this phenomenon have puzzled researchers for more than a quarter of a century. In view of recent experimental findings, it is postulated that the phenomenon can be attributed to the variability of rust composition and density, specifically to the variable ratio of the mass fractions of iron oxide and iron hydroxide-oxide, which is affected by the magnitude of the applied corrosion current density. Based on this hypothesis, a corrosion-induced cracking model for virtual impressed-current testing is presented. The simulation results obtained with the proposed model are validated against experimental data, showing good agreement. Importantly, the model can predict corrosion-induced cracking under natural conditions and thus allows for the calculation of a newly proposed crack width slope correction factor, which extrapolates the surface crack width measured during accelerated impressed current tests to corrosion in natural conditions.


Introduction
Corrosion is responsible for premature deterioration of 70-90% of all reinforced concrete structures [1,2].One of the most common degradation mechanisms is corrosion-induced cracking, which takes years or decades under natural conditions.For this reason, researchers have strived for decades to accurately simulate and predict corrosion-induced cracking [3].
Currently, there are two methods employed to experimentally simulate corrosion-induced cracking in laboratory conditions [4].The first option is to expose the concrete specimen to an exaggerated climate, i.e., to various intensities of contact with aggressive species such as chlorides or carbon dioxide, and possibly also to variable humidity (wetting-drying cycles).This accelerates the corrosion process to a certain extent.However, accelerated testing can still require months or years to achieve significant corrosion-induced cracking.Also, corrosion current density, steel mass loss and anodic area are evolving dynamically with the penetration of aggressive species.These are not known in advance and need to be measured in other ways, which often requires many samples to be tested at different stages of the corrosion process.
The other method to simulate corrosion-induced cracking experimentally is the impressed current technique.In this case, an electric potential is applied on the steel rebar which serves as the anode, with typically another stainless steel rebar or mesh serving as the cathode.This method results in nearly uniform corrosion and allows to accurately control corrosion current density and mass loss.At the same time, it allows for applying high corrosion current density in hundreds or thousands of µA/cm 2 , which drastically shortens the time of the experiment to days or hours.For all their benefits and simplicity, impressed current techniques have become very popular [3].However, a major downside is that while corrosion current densities from 100 to 2000 µA/cm 2 are commonly applied in impressed current tests [3], natural corrosion current densities are much smaller, typically about 1 µA/cm 2 [5][6][7][8][9].This raises the question of how accurately can the impressed current technique reproduce the corrosion process and corrosion-induced cracking under natural conditions.
To answer this question, Alonso et al. [10] measured the surface crack width on the impressed current test with varying applied corrosion current density and found out that lower current density resulted in faster crack propagation with respect to corrosion penetration (i.e. the thickness of the corroded steel layer).The measured differences were significant: Alonso et al. [10] reported that the slope of linearly fitted crack width/corrosion penetration curve for the corrosion rate of 10 µA/cm 2 was even six times larger than for the corrosion rate of 100 µA/cm 2 .The decreasing trend of the slope of linearly fitted crack width/corrosion penetration curves with increasing applied corrosion current density was confirmed by the comprehensive experimental study of Pedrosa and Andrade [11], who found out that the slope is proportional to a negative power of the corrosion current density.Slower propagation of cracks under an increased applied corrosion current density was also observed by Vu et al. [12], Saifullah and Clark [13] and Mullard and Stewart [14] up to approximately 200 µA/cm 2 .By fitting the experimental results, Vu et al. [12] proposed a loading correction factor dependent on the applied corrosion current density to allow for the extrapolation of the time to cracking measured within accelerated tests to corrosion in natural conditions.It has to be mentioned that for very high currents over 200 µA/cm 2 , some authors reported an increased speed of crack propagation compared to current densities lower than 200 µA/cm 2 [13,15,16].El Maaddawy and Soudki [15] suggested that for a higher corrosion rate, corrosion products have less time to be transported into the concrete pore space and a larger portion of them is thus trapped at the steel-concrete interface, which increases the induced pressure.
Even though the decrease of crack propagation with increasing applied corrosion current up to at least 200 µA/cm 2 is well experimentally documented, there has been no agreement on the origins of this phenomenon.Alonso et al. [10] and Pedrosa and Andrade [11] proposed an alternative explanation, suggesting that the loading rate-dependent material properties of concrete may be responsible.The idea is that the increased strain rate would affect the fracture-related properties such that the tensile strength would increase and delay the cracking process.If this was true, given the range of rates relevant to the corrosion process, the rate-dependency of material properties would have to be caused by creep.Although creep is likely to affect the results of long-term tests conducted under natural-like corrosion current densities, the study of ?] found out that creep alone does not explain the dependence of the speed of the crack propagation (with respect to corrosion penetration) on the corrosion rate.
Recently, a new plausible explanation emerged from the study of Zhang et al. [17] who performed X-ray diffraction (XRD) measurements on the series of rust samples from impressed current tests conducted with different applied corrosion current densities.Zhang et al. [17] found out that the composition of rust, specifically the mass fractions of two main components of rust-iron oxides (mainly wüstite FeO, maghemite (γ − Fe 2 O 3 ) and magnetite (Fe 3 O 4 )) and iron hydroxyoxides (mainly goethite α-FeOOH, akageneite β-FeOOH, lepidocrocite γ-FeOOH)-changed with the applied magnitude of corrosion current.This indicates that as the corrosion current density increases, the mass fraction of iron hydroxy-oxides decreases.Because iron hydroxy-oxides have a smaller density than iron oxides, the decrease in hydroxy-oxide content results in lower corrosion-induced pressures and thus slower crack evolution.Although the study of Liu et al. [18] identified high hydroxy-oxide content even when the corrosion current density was raised to 100 µA/cm 2 , the conclusions of Zhang et al. [17] are supported by a number of other experimental studies.Yuan et al. [4] observed different colours (indicating different chemical compositions) of rust produced from accelerated impressed current tests than those corroded in an artificial climate environment closely simulating the natural corrosion process.Chitty et al. [19] investigated rust from ancient Gallo-Roman reinforced concrete artefacts naturally corroding for roughly two millennia and identified mostly hydroxy-oxide goethite with marblings of iron oxides maghemite and magnetite.The analyses of naturally carbonated corrosion products in carbonated reinforced concrete facade panels by Köliö et al. [20] identified mostly hydroxy-oxides goethite, feroxyhite (δ − FeOOH) and lepidocrocite with some rare findings of iron oxides magnetite and maghemite.
A very high content (more than 56% mass fraction) of iron oxides in four impressed current tests accelerated to 100 µA/cm 2 was measured by Zhang et al. [21].
These recent findings are adopted to build our hypothesis and develop a corrosion-induced cracking model capable of resolving (i) the reactive transport of dissolved iron species released from the steel surface and their precipitation into the dense rust layer adjacent to the corroding steel surface and the concrete pore space, (ii) the change of the rust chemical composition, specifically the ratio of mass fractions of iron oxides and iron hydroxy-oxides, and thus of the density of rust with the applied magnitude of corrosion current density, (iii) the corrosion-induced pressure of corrosion precipitates accumulating in the dense rust layer and the concrete pore space of concrete, and (iv) the corrosion-induced fracture of concrete, as simulated with a quasi-brittle phase-field fracture model.For the first time, the analytically derived estimates are presented: for (I) the corrosion-induced pressure of the dense layer of corrosion products, which takes into account their compressibility, and (II) the Fe 2+ ions flux reduction factor determining the ratio of corrosion products precipitating in the dense rust layer and those transported further into the concrete pore space.Thus, a comprehensive model capable of simulating impressed current tests is provided.Also, the model allows to predict a newly proposed crack width slope correction factor with which the surface crack obtained from impressed current tests can be extrapolated to corrosion in natural conditions.

Theory and computational framework
The formulation of the corrosion-induced cracking model comprises several parts.Firstly, the model for reactive transport and precipitation of dissolved iron released from the corroding steel surface is described in Section 2.1, allowing for the prediction of the distribution of iron oxide and iron hydroxy-oxide rust.This is followed by the discussion of the phase-field fracture model in Section 2.2.Corrosion-induced mechanical stress arises from the accumulation of rust in the form of a dense rust layer (Section 2.3) and in the pore space of concrete (Section 2.4).The damage-dependent diffusivity tensor allowing for the consideration of the role of crack in the reactive transport of dissolved iron species is introduced in Section 2.5.The formulation of the model is concluded in Section 2.6 with the summary of governing equations and related boundary conditions.
Notation.Scalar quantities are represented by light-faced italic letters, e.g.φ, Cartesian vectors by upright bold letters, e.g.u, and Cartesian second-and higher-order tensors by bold italic letters, e.g.σ.The symbol 1 denotes the second-order identity tensor and I represents the fourthorder identity tensor.Inner products are denoted by a number of vertically stacked dots, where the number of dots corresponds to the number of indices over which summation takes place, e.g.
Finally, ∇ and ∇• respectively denote the gradient and divergence operators.

Solution domain, representative volume element (RVE) and primary unknown variables
The problem is solved on a domain Ω = Ω c ∪ Ω s , which comprises the concrete domain Ω c and the embedded steel domain Ω s (see Fig. 1).The outer boundary of Ω and Ω c is Γ c and the inner boundary of Ω c between the steel and concrete domains is Γ s .The symbol n denotes the outer normal vector to the boundary Γ c ∪ Γ s of the concrete domain.Concrete is a porous material and within the reactive transport model, its internal structure is described by the representative volume element (RVE) depicted in Fig. 2. The concrete RVE of volume V is assumed to comprise the solid concrete matrix of volume V s and the pore space of volume V − V s , characterized by capillary porosity p 0 = (V − V s )/V , with all these quantities remaining unchanged in time.The pore space is simplified to be fully saturated with water, which is a reasonable assumption in the vicinity of the steel rebar, especially for corrosion induced by impressed current or chloride ingress.There, the pore space is gradually filled with precipitated rust, and the volume of liquid pore solution V l decreases to the benefit of the increasing volume of precipitated rust V r .Volume fractions θ l = V l /V , θ r = V r /V and saturation ratios S l = θ l /p 0 , S r = θ r /p 0 of the liquid pore solution and the precipitated rust provide useful indicators of their current volume content.The primary unknown variables of the reactive transport problem are c II and c III , in mols per cubic meter of the liquid pore solution.

Solid concrete matrix
Liquid pore solution

Distribution and chemical composition of rust in concrete
Microscopy investigations of corroded reinforced concrete samples [19,22,23] show that a portion of the rust accumulates in the immediate vicinity of the corroded steel rebar as a dense rust layer while the remaining rust is located in the concrete pore space up to a certain distance from the rebar (see Fig. 2).
In samples corroded under both natural and accelerated conditions, the main corrosion products found with X-ray diffraction (XRD) measurements were iron oxides (wüstite FeO, hematite α − Fe 2 O 3 , magnetite Fe 3 O 4 ) and iron hydroxy-oxides (goethite α−FeO(OH), akaganeite β−FeO(OH), lepidocrocite γ − FeO(OH), feroxyhyte δ − FeOOH) [17,19,24].It should be noted that the XRD method can determine only the presence of crystalline materials, but some corrosion products are known to be amorphous [24].One of the critical inputs of the model is the dependence of the mass Fig. 3: The dependence of the mass fraction of hydroxy-oxides rust r h on corrosion current density.r h decreases dramatically when accelerated from the typical range of natural corrosion current densities about 1 µA/cm 2 .Data adopted from Zhang et al. [17].
fraction of iron hydroxy-oxides in rust on applied corrosion current density because this determines the density of produced rust and thus the induced pressure.This curve is obtained by fitting the experimentally measured data of Zhang et al. [17] with a power function k 1 (i a /i a,ref ) k 2 (see Fig. 3), where i a,ref = 1 µA/cm 2 is the reference corrosion current density.The first measurement of Zhang et al. [17] came from a naturally corroded sample for which the value of corrosion current density was not documented and thus it was estimated as 1 µA/cm 2 , which is a typical value measured during the natural corrosion process in reinforced concrete [5][6][7][8][9].Even though experimental measurements of the mass ratio of particular corrosion products are very scarce in the current literature, the results of Zhang et al. [17] agree with the studies on samples corroded under natural conditions [19,20] and with accelerated impressed current tests [21].In Fig. 3one can see that r h drops dramatically from 0.9 for 1 µA/cm 2 to roughly 0.5 for 50 µA/cm 2 , while a subsequent decrease is significantly more moderate.

Iron ion transport and precipitation
An alkaline concrete pore solution allows for the formation of a nanometre-thick protective semiconductive layer on the steel rebar surface, which significantly impedes the corrosion process.
However, this layer can be disrupted, typically by chlorides or due to the carbonation of concrete.
If water and oxygen are present, corrosion can initiate.A corrosion current of density i a starts to flow between anodic and cathodic regions and Fe 2+ ions are released from the steel surface into the pore solution by reaction (1a).The dissolved ferrous Fe 2+ ions then undergo a series of complex chemical reactions forming a number of intermediary products, see for instance [25,26] for more details.For the purpose of mechanical modelling, the complexity of the chemical reaction system is reduced to three arguably critical reactions (see Fig. 4).The first considered option for Fe 2+ ions is to be oxidised to Fe 3+ ions (reaction (1b)), which then further precipitate into iron hydroxyoxides (FeO(OH)); see reaction (1c).The alternative option to this process is direct precipitation of Fe 2+ ions into Fe(OH) 2 , which is assumed to be further transformed into iron oxides, especially into mixed Fe 2+ -Fe 3+ oxide Fe 3 O 4 and Fe 2 O 3 [27], though this transformation is not explicitly Experimental studies indicate that iron oxides Fe 2 O 3 and Fe 3 O 4 together with iron hydroxy-oxides FeO(OH) constitute the majority of produced rust for both accelerated impressed current tests and samples corroding in natural conditions [17-20, 24, 28].However, critical for this study, based on the experimental results of Zhang et al. [17], it is assumed that the ratio of mass fractions of iron oxides and hydroxy-oxides changes with the corrosion current density (see Fig. 3).The reactive Fig. 4: Scheme of considered chemical reactions transport of Fe 2+ and Fe 3+ in Ω c is described by mass-conserving diffusion equations previously derived in [29].Small deformations are assumed, the velocity of the solid concrete matrix is neglected, and the flux term is scaled with liquid volume fraction following Marchand et al. [30].
In Eqs. the corrosion current density and is estimated such that the mass ratio documented in Fig. 3 is achieved for a closed system with a given initial concentration of Fe 2+ ions.It can be shown that, in the limit of t → ∞, the ratio M o and M h refer to the molar masses of oxide rust and hydroxy-oxide rust respectively.These estimates are theoretically accurate only for large times, but the half-time of Fe 2+ ion transformation is only about 2900 s for i a = 1 µA/cm 2 and steadily decreases to approximately 560 s for i a = 500 µA/cm 2 .This is very short compared to the advance of corrosion penetration-the shortest further discussed case studies lasted over 26 days or 6 hours for the cases of i a = 1 or i a = 500 µA/cm 2 , respectively.For all calculations, the constant oxygen concentration c ox = 0.28 mol m −3 is assumed as in the study of Zhang and Angst [33].
The volume fractions θ o and θ h of precipitated iron oxides and hydroxy-oxides are respectively governed by in which ρ o and ρ h are densities of oxide rust and hydroxy-oxide rust respectively, and M h and M o denote the molar masses of hydroxy-oxide rust and oxide rust, respectively.Immobile precipitated rust gradually fills the concrete pore space and reduces the liquid volume fraction Reactive transport equations ( 2)-( 3) require boundary conditions.On the outer concrete boundary, zero flux is assumed such that n on Γ c .On the corroding steel rebar surface Γ s , the inward influx of ferrous ions follows Faraday's law where i a is the corrosion current density, F is Faraday's constant and z a = 2 represents the number of electrons exchanged in anodic reaction (1a) per one atom of iron.Some of the released Fe 2+ ions precipitate into a dense rust layer instead of being released further into concrete pore space.This effect is taken into consideration with the flux reduction coefficient k f ∈ 0, 1 , which is discussed in the following section.

Flux reduction coefficient k f
During corrosion, ferrous ions are released from the steel surface by electrochemical reaction (1a) and their flux follows Faraday's law J II,F ar = i a /(z a F ).A portion of the released Fe 2+ ions and emerging Fe 3+ ions precipitate to form a dense rust layer in the vicinity of the rebar while the remainder is transported further into the pore space to eventually precipitate there (see Fig. 5).
The thickness of the rust layer is typically much smaller compared to the rebar radius.For this reason, the rust layer is not explicitly geometrically considered during domain discretization for numerical calculation, as this would strongly exacerbate the computational requirements.Instead, the steel domain is handled as unchanged in time, and the mechanical pressure of the accumulating dense rust layer on concrete is taken into consideration by additional pressure p on the steel boundary, further discussed in Section 2.3.Thus, in reality, the flux of ferrous ions prescribed on the steel boundary is only the portion of this flux that did not precipitate into the dense rust layer.
Equation ( 7) is solved on an evolving domain 0, t r + t c , which is separated into the dense rust layer subdomain 0, t r ) with diffusivity D r,II and the concrete subdomain t r , t r + t c ) with diffusivity D c,II , where t r changes in time (see Fig. 5).For the sake of simplicity, the thickness of the dense rust layer is considered to be equal to the thickness of the corroded steel layer predicted by Faraday's law.
It is possible to solve equations (7a) and (7b) independently and analytically, considering Faraday's law dictated flux J II,F ar at x = 0, c II = 0 at x = t r + t c and a given unknown flux value at the current interface x = t r .From the assumption of flux continuity at x = t r , a solution for k f is obtained.If D c,II is further replaced by S l D c,II to consider the effect of pore-clogging with precipitates on Fe 2+ diffusivity, the resulting flux can be expressed as where the flux reduction coefficient is calculated using rust-related and concrete-related constants Based on the results of our previous study [29], it is set t c = 2 mm.

Phase-field description of precipitation-induced cracks
In this model, the phase-field cohesive zone model (PF-CZM) of Wu [34,35] is employed, which allows for the calibration of the softening curve such that the quasi-brittle fracture nature of concrete is captured.Here a brief summary of the governing equations of the PF-CZM model and the necessary background is provided.For the derivation and more detailed explanation, see section 2.4 in Ref. [29].The fracture problem is solved exclusively on the concrete domain Ω c while the steel in Ω s is assumed to remain linear elastic.The primary unknown variables in Ω c are the displacement vector and the phase-field variable where Γ u ⊂ (Γ c ∪ Γ s ) is the portion of concrete boundary with prescribed displacement u.In (11) and ( 12), d = 2, 3 is the geometric dimension of the problem and function space W 1,2 (Ω c ) d is the Cartesian product of d Sobolev spaces W 1,2 (Ω c ) consisting of functions with square-integrable weak derivatives.The phase-field variable φ characterises the current damage of concrete such that φ = 0 for undamaged material and φ = 1 for a fully cracked material.Note that φ is not the traditional damage variable that directly reflects the relative change of stiffness, but it is linked to the material integrity through the degradation function, to be introduced later.
The small-strain tensor ε = ∇ s u = (∇u + (∇u) T )/2 is used, and the total strain ε = ε e + ε ⋆ is additively decomposed into the elastic part ε e and the precipitation eigenstrain ε ⋆ .Because damage limits the stress that can be carried by the material, the Cauchy stress tensor is calculated as where C c e is the fourth-order isotropic elastic stiffness tensor of concrete, characterized by Young's modulus E and Poisson's ratio ν c , and g(φ) is the degradation function that describes the remaining integrity of the material.Typically, it is a function that satisfies conditions g(0) = 1 and g(1) = 0 and is non-increasing and continuously differentiable in [0, 1].
Assuming tractions t prescribed on Γ t = (Γ c ∪ Γ s ) \ Γ u and a volume force b acting on Ω c , the governing equations and the related boundary conditions for the coupled damage-displacement problem read where G f is the fracture energy and ℓ is the characteristic phase-field length scale governing the size of the process zone [36].To make the results independent of ℓ, it is recommended to use ℓ ≤ min (8ℓ irw /3π, L/100 ∼ L/50) where L is a characteristic length of the structure (e.g., the beam depth) and ℓ irw = EG f /f 2 t is the Irwin internal length, evaluated from the fracture energy, tensile strength f t and elongation modulus In the crack process zone, the characteristic size of the finite elements must be sufficiently small to provide good resolution of high strain gradients, typically 5-7 times smaller than ℓ [36].In (14c) the crack driving force history function is employed to secure the irreversibility of damage [37] ).The softening curve resulting from the phase-field model is affected by the choice of the degradation function g(φ).For the PF-CZM model, Wu [34] proposed to use This function is calibrated by an appropriate choice of parameters p ≥ 2, a 1 > 0, a 2 and a 3 , based on the required shape of the softening curve, characterized by the ratios where w c is the crack opening at zero stress (full softening), k 0 is the initial slope of the softening curve (i.e., its negative slope at the onset of cracking), and are the reference values that would correspond to linear softening.
For softening laws with zero stress attained at a finite value of crack opening, the parameter p is set to 2, and for laws that approach zero stress only asymptotically, a value greater than 2 is appropriate.The other parameters of the degradation function (16) are then evaluated as For concrete, the experimentally measured softening can usually be reasonably represented by the Hordijk-Cornelissen cohesive law [38], for which and thus β w = 5.136/2 = 2.568 and β k = 1.3546/0.5 = 2.7092.Wu [34] demonstrated that his version of the phase-field model indeed provides a close approximation of the Hordijk-Cornelissen curve if the parameters of the degradation function ( 16) are set to These are the values adopted in the present study.The parameter a 1 , linked to the choice of the characteristic length ℓ, will be evaluated from the first formula in Eq. ( 19).Rust has a lower density than steel and the dense rust layer thus tends to occupy a larger volume (characterised by its thickness t r ) than was vacated by the corrosion process (characterised by t cor ).

Steel
However, the rust volumetric expansion is constrained by the surrounding concrete matrix.Thus, the rust layer exerts pressure p on the concrete boundary, which leads to displacement u c such that the dense rust layer occupies volume characterised by thickness t r = t cor + u c .

Pressure of the accumulated rust layer
As was discussed in previous sections, Fe 2+ and Fe 3+ ions eventually precipitate into rust that either forms a dense rust layer in the vicinity of steel rebar or gradually fills the concrete pore space.Rust has a significantly lower density than original iron, typically by 3 -6 times [39].
For this reason, if the volume characterised by thickness t cor is vacated by the corrosion process, the dense rust layer would occupy 3 -6 times larger volume under unconstrained conditions (see Fig. 6).However, the rust volumetric expansion is constrained by the concrete matrix and the expansion of the rust layer thus exerts pressure p on the concrete boundary.Constrained rust volumetric expansion in the pore space also causes pressure on pore walls, which is described by the precipitation eigenstrain, as explained in the following section.
To estimate the pressure p of a dense rust layer on concrete, let us virtually isolate the cylinder containing a rebar of radius a and the adjacent layer of concrete with radius b = αa from the studied concrete specimen (see Fig. A.1). Assuming the length of the cylinder L ≫ max(a, b), isotropic linear elasticity of concrete and symmetric boundary conditions, the thick-walled cylinder can be described as an axisymmetric linear elastic plane-strain problem.If the damage is simplified to be uniformly distributed across the cylinder, the displacement on the inner concrete boundary u c for a given pressure p reads (the derivation is provided in Appendix A) where ν c is Poisson's ratio of concrete.E c,d = g(φ)E c is the damaged secant modulus of concrete, calculated as the product of Young's modulus of concrete E c and degradation function g(φ) (see Section 2.2 for more details) evaluated on the inner concrete boundary.The concrete compliance constant C c thus encapsulates the impact of concrete sample geometry and material properties.
Both u c and p in Eq. ( 23) are unknown, and thus it is necessary to link them by an additional equation.For this reason, let us now analyse the relation between pressure p compressing the rust layer and the volume V r,d , which is (for small thicknesses) proportional to t r = u c + t cor .Let us assume that the bulk modulus of rust K r is constant and equal to E r /(3(1 − 2ν r )), with E r and ν r being respectively the Young modulus and Poisson ratio of rust.The bulk modulus is the reciprocal value of the compressibility of rust, which is under isothermal conditions expressed as Integration of (24) leads to where V 0 is an integration constant, which has the meaning of the volume that would be occupied by the rust at zero pressure.By inversion, it is possible to express the pressure as The rust layer of thickness t r = t cor + u c occupies volume V r,d = π((a + u c ) 2 − (a − t cor ) 2 )L.The corrosion penetration is usually much smaller than the bar diameter.The assumption that a ≫ t r leads to V r,d ≈ 2π(t cor + u c )aL = 2πt r aL.Furthermore, the rust volume at zero pressure can be expressed as is a volumetric expansion coefficient resulting from the disparity of molar volumes of steel and rust.This coefficient is evaluated from the mass fractions r h and 1 − r h of iron hydroxy-oxides and iron oxides in rust produced at the given corrosion current density (see Fig. 3) and the corresponding molar volume ratios κ o and κ h of oxide rust/iron and hydroxy-oxide rust/iron, respectively.Substituting the expressions for V 0 and V r,d into Eq.( 26) results in p = K r ln κt cor t cor + u c (28) Using the previously derived relation (23), it is easy to eliminate the unknown pressure p from Eq. ( 28) and construct an equation for a single unknown, u c : This can be rewritten in the simple form in which F = κt cor e tcor/(CcKr) /(C c K r ) is a known right-hand side and the left-hand side is a nonlinear function of unknown G = (u c + t cor )/(C c K r ).For F ≥ 0, Eq. ( 30) has a unique real solution G = W (F ), where W is the so-called Lambert function, also known as the product logarithm.Since F is always positive for our problem, the solution of Eq. ( 29) can be written as The corresponding pressure p is easily expressed as p = u c /C c .
Let us note that even though the formula (31) has been derived assuming axial symmetry, the real distribution of damage is not uniform on the boundary between steel and concrete due to the gradual localisation of cracks.For this reason, the assumption of axisymmetry will be violated and the formula will be evaluated point-wise on the steel-concrete boundary.

Precipitation eigenstrain
Fe 2+ and Fe 3+ ions transported into the pore space gradually precipitate there into iron oxides and hydroxy-oxides.As was the case for the dense rust layer adjacent to the rebar, rust accumulates in the confined conditions of concrete pore space and exerts pressure on concrete pore walls [3,39] (see Fig. 2 for a graphical illustration).This statement is supported both by the similarity with the well-described damage mechanism of precipitating salts in porous materials [40][41][42][43][44][45][46] and the recent simulations of Korec et al. [29], which revealed that this mechanism could considerably contribute to corrosion-induced fracture in its early stages until the pore space surrounding rebar is locally filled.Macroscopic stress resulting from the accumulation of rust in the pore space under confined conditions is simulated with the precipitation eigenstrain proposed by Korec et al. [29], based on the previous works of Coussy [43], Krajcinovic et al. [47] and Mura [48].The precipitation eigenstrain is calculated as where is the bulk modulus of rust-filled concrete calculated from Young's modulus of rust-filled concrete E and concrete Poisson ratio ν c .Because the mechanical properties of rust and concrete are very different, E = (1 − θ r )E c + θ r E r is interpolated between Young's modulus of rust-free concrete E c and rust Young modulus E r following the rule of mixtures.

Damage-dependent diffusivity tensor
Dissolved Fe 2+ and Fe 3+ ions are transported much more easily through emerging cracks than through the surrounding concrete matrix.For this reason, the damage-dependent diffusivity tensor [29,49] is considered, and defined as where D m,α is the diffusivity of the considered species in concrete and D c,α ≫ D m,α controls the diffusivity of the cracked material.

Overview of the governing equations
The solved coupled problem comprises six differential equations (34a)-(34f) for six unknown variables -displacement vector u, phase-field variable φ, Fe 2+ ions concentration c II , Fe 3+ ions concentration c III , hydroxy-oxide rust volume fraction θ h and oxide rust volume fraction θ o : Differential equations (34a)-(34d) are associated (in the presented order) with boundary conditions Boundary conditions for equations (34e) and (34f) are not required because these equations do not contain any spatial derivatives.
The resulting coupled system of differential equations is numerically solved with the finite element method using COMSOL Multiphysics software.The domain Ω c ∪Ω s is discretised with linear triangular elements and a staggered solution scheme is employed.To achieve mesh-independent results, ℓ = 3 mm is chosen such that ℓ would be five times larger than the maximum element size [36].

Results
After the presentation of parameter values employed for numerical simulations (Section 3.1), the observed characteristic features of the model are discussed in Section 3.2, together with the implications for the understanding of corrosion-induced cracking.Then, in Section 3.3, impressed current tests with varying corrosion current density from the comprehensive experimental study of Pedrosa and Andrade [11] are simulated.The predicted relations between the surface crack and the thickness of the corroded steel layer are fitted with linear functions and their slopes are compared with the experimental measurements, revealing the ability of the model to capture the decrease of the slope with increasing corrosion current density.Based on the obtained results, a correction factor to be applied on the crack width slope obtained from accelerated impressed current tests is proposed, allowing for extrapolation of results to naturally corroding structures.

Choice of model parameters
The mechanical parameters listed in Table 1 correspond to the concrete specimens employed in the comprehensive experimental study by Pedrosa and Andrade [11], in which the corrosioninduced cracking under different values of applied corrosion current density was analysed.The 15 x 15 x 50 cm samples reinforced with a 16 mm rebar were separated into two groups with curing times of 28 and 147 days, which led to different mechanical properties.Pedrosa and Andrade [11] measured the tensile and compressive strength.The fracture energy is estimated using the formula of Bažant and Becq-Giraudon [50].For the steel rebar, Young's modulus of 205 GPa and Poisson's ratio of 0.28 are assumed, as these are common values for steel.Fracture energy Tab. 1: Model parameters: mechanical properties of concrete, based on the measurements by Pedrosa and Andrade [11] and literature data.
Only the capillary porosity of cement paste is considered to be relevant in the vicinity of rebar, and it is estimated using the model of Powers and Brownyard [52], assuming the degree of hydration of 0.9.The values of the initial diffusivity of Fe 2+ ions and Fe 3+ ions in concrete are adopted from the study of Stefanoni et al. [31], where they were estimated from the known diffusivity in solution assuming that the ratio between the diffusivity in water and concrete is the same as for chloride ions.The diffusivity of iron ions in cracks was considered identical to the diffusivity in the pore solution.
Molar volume ratios of iron oxide rust/iron and hydroxy-oxide rust/iron characterising rust volumetric expansion are chosen in the range of values reported by Vu et al. [12], Zhao and Jin [27] and Zhang et al. [17] for commonly present iron oxides and iron hydroxy-oxides.The measurements of Young's modulus and Poisson's ratio of rust are infamously scattered in the current literature, so two intermediate values within the range reported by Zhao and Jin [27] are considered; they were identified to lead to an accurate prediction of crack width in the previous study of the authors [29].The diffusivity for iron ions in rust was adopted from the study of Ansari et al. [53].
Let us also stress that in numerical simulations conducted within this study, uniform corro- sion and thus the constant value of corrosion current density i a are considered.This is because impressed current tests typically ensure uniform corrosion by artificially introducing high chloride concentration into concrete, as was done in tests by Pedrosa and Andrade [11] considered in this study.For the results on the simulation of a non-uniform chloride-induced corrosion under natural conditions with a similar model as in this study, see Ref. [54].

General aspects of the model (the impact of the dense rust layer)
In the proposed model, rust precipitates either in the dense rust layer adjacent to the corroding rebar surface or deeper in the pore space of concrete.The flux of Fe 2+ ions released from the corroding steel surface according to Faraday's law is divided between these two processes based on the flux reduction coefficient k f ∈ [0, 1].As can be seen in Fig. 8a, k f = 1 if there is no dense rust layer and all released ferrous ions are transported into the concrete pore space.With the advance of the corrosion process, the thickness of the dense rust layer increases, decreasing the portion of Fe 2+ ion flux that is released to the concrete pore space.For a thick rust layer, J II eventually vanishes, leaving all released Fe 2+ ions to contribute to the build-up of the dense rust layer.Also, as saturation of the pore space with precipitates increases, the Fe 2+ ion diffusivity in concrete is of the corroded steel layer but decreases when concrete becomes more compliant due to damage.In many currently available models, rust is often simplified to be incompressible.However, the stiffness of rust has a profound impact on pressure predictions (b).Results in (c) reveal that the assumption of rust incompressibility leads to nearly identical results as if rust had the same elastic properties as steel.
reduced, which decreases the ability of Fe 2+ ions to escape the dense rust layer.
The accumulation of precipitated rust in (A) a dense layer of corrosion products in the space vacated by steel corrosion, and (B) the concrete pore space under constrained conditions, induces pressure on the concrete matrix.Pressure p in the dense rust layer increases with its thickness t cor (see Figs. 7a and 7c).However, it also non-linearly decreases with damage, since concrete becomes more compliant and the expansion of rust is then less constrained (Figs.7a and 7b).The mechanical properties of rust have a strong influence on p (see the dependence on the bulk modulus of rust in Fig. 7b).If the bulk modulus of rust was similar to steel, pressure p would decrease almost in proportion to the degradation factor g(φ).However, for lower bulk moduli of rust, in tens or hundreds of MPa, the decrease deviates from linearity and the induced pressure p is significantly smaller.Many currently available models neglect the impact of rust compressibility and simplify the steel-concrete boundary displacement to be u c = t cor (κ − 1) as if rust was incompressible.As demonstrated in Fig. 7c, this assumption leads to a nearly equivalent displacement prediction as if rust elastic properties were equivalent to those of steel.Since the experimental measurement of Young's modulus and Poisson's ratio of rust are quite scattered [27] in the range from 100 MPa to 360 GPa, our results stress the need for accurate characterization studies.
From the conducted simulations, it was observed that the distribution of rust in pore space differed for natural-like low corrosion current densities (about 1 µA/cm 2 and lower) compared to high current densities (tens or hundreds of µA/cm 2 ) typical for corrosion induced by accelerated impressed current (see Fig. 8b).In the case of high current densities, the rust was distributed closer to the steel surface than for low current densities, which led to a significantly higher rust saturation of pores in the vicinity of steel rebar.Thus, pores are blocked with rust earlier than in natural conditions, not allowing for additional transport of dissolved iron species into pore space.This suggests that during the natural corrosion process, significantly more dissolved iron species can escape to the concrete pore space and precipitate there.The rapid increase of the rust saturation for high current densities probably results from significantly higher concentrations of released Fe 2+ ions, which accelerates the precipitation reactions (1c) and (1d).Similar results were obtained in the computational study of Stefanoni et al. [31].For high corrosion current densities in accelerated tests, rust accumulates closer to the steel surface than during corrosion in natural conditions, and thus the maximum rust saturation of pores is higher.

The impact of the magnitude of applied corrosion current density on crack width-analysis and validation
Fig. 9c shows the geometry and the typical fracture pattern observed for the simulated impressed current tests.The main crack propagates across the shortest distance between the concrete surface and the steel rebar, perpendicularly to the concrete surface, and opens wide with the ongoing corrosion process.In addition, two other lateral cracks develop.The surface crack width is calculated as the integral of the x-component (i.e., the normal component in the direction perpendicular to the crack) of the inelastic strain tensor ε d = ε − ε e − ε ⋆ over the upper concrete surface [55].The inelastic strain is obtained by subtracting the eigenstrain ε ⋆ and the elastic strain ε e = (C c e ) −1 : σ = g(φ)(ε − ε ⋆ ) from the total strain.For the same thickness of corroded steel layer t cor , the crack width decreases with increasing applied corrosion current density.This trend agrees with the experimental observations of [10-12, 14-16, 56] and is the result of the hypothesis that the chemical composition of rust, in particular the mass ratio of iron oxides and iron hydroxy-oxides, changes with the magnitude of the applied current density (see Fig. 3) such that the rust density increases with increasing current density, reducing the pressure caused by constrained expansion.
The predicted crack width is linearly fitted by w = βt cor + ζ, so that the crack width slope β can be compared with its experimental counterpart measured by Pedrosa and Andrade [11], who employed the same fitting procedure on their experimentally measured crack widths.To prevent the distortion caused by nonlinearity for zero and very small crack widths, only w ≥ 1 µm was considered for fitting.Though linear fitting approximates the results very accurately for the thinner concrete of 20 mm (Fig. 9a), the crack width increases superlinearly with t cor for the thicker cover of 40 mm, which hinders the accuracy of the linear fitting procedure.For thicker concrete covers, the predicted crack width slope β is thus affected by the maximum evaluated t cor and linear fitting is not optimal.In the discussed simulations, the maximal t cor was 8 µm and 12 µm for the cases of 20 and 40 mm thick concrete cover, respectively.
The comparison of the predicted crack width slope β with the results recovered from the fitting of experimental measurements [11] reveals that the predicted data lie in the experimental range and that the trend characterised by the decay with a negative power of current density [11] is reproduced very well.For a better understanding of how well the various sets of predicted data agree with this trend, the predicted results are fitted with β = γ 1 (i a /i a,ref ) γ 2 , where i a,ref = 1 µA/cm 2 is the reference corrosion current density.In Fig. ??, an excellent agreement with the data for the smaller concrete cover of 20 mm can be observed.On the other hand, the predictions for the thicker concrete cover of 40 mm are significantly more scattered, even though the overall trend still agrees.Fig. 9a and Fig. 9b suggest that the accuracy of the predictions is clearly related to how well can the data be fitted with a linear function.While data for the 20 mm cover are fitted very well, the non-linear dependency for the thicker 40 mm cover diminishes the representative value of β.Also, it was observed that with the increasing thickness of the concrete cover, the model becomes numerically more sensitive such that there are several similar crack patterns that lead to a slightly different evolution of the surface crack width.The realization of a particular crack pattern appears to be triggered by small numerical differences.
Critically, the proposed model captures the initially rapid decrease of β with increasing current density i a suggested by the experimental data of Pedrosa and Andrade [11].For current densities larger than approximately 50 µA/cm 2 , the decrease of β is more moderate.As can be seen in Fig. ??, this agrees qualitatively very well with the experiments.However, the magnitude of slopes observed by Pedrosa and Andrade [11] for current densities smaller than 10 µA/cm 2 were even several times larger than those recovered from numerical simulations.Though the paper of Pedrosa and Andrade [11] is currently the most comprehensive experimental study on the impact of the magnitude of impressed current on the crack width slope, the total number of tested specimen was arguably small (11 in total and every test was conducted only once which does not allow to evaluate experimental error).Thus, it is not clear whether the observed discrepancy between the quantitative values of numerical predictions and experiments can be attributed to the experimental error or some other phenomena unaccounted by the model.Clearly, more comprehensive experimental campaigns are needed.For current densities larger than 10 µA/cm , although the role of the pressure of rust precipitating in the pore space of concrete is also not negligible.However, one has to bear in mind that the proposed model overestimates the corrosion-induced pressure of the dense rust layer, as the entire volume of steel vacated by corrosion is assumed to be filled with rust regardless of the portion of ferrous ions escaping into pore space.In reality, the pressure of the dense rust layer can be reasonably expected to be smaller, raising the relative impact of the pressure of precipitates accumulating in the pore space.
In Fig. 11b, the relative error of slope β, caused by the acceleration of the applied current compared to corrosion in natural conditions, is evaluated in terms of newly proposed crack width slope correction factor k β inspired by the loading correction factor of Vu et al. [12] for time-tocracking.The correction factor k β was calculated as k β = (i a /i a,ref ) γ 2 .The factor k β can thus be interpreted as the relative error caused by the acceleration of the impressed current tests compared to results that could be expected in natural conditions when i a ≈ 1 µA/cm 2 .One can see that the error induced by current acceleration strongly increases with the thickness of the concrete cover.
While for 147 days cured samples and i a = 500 µA/cm 2 , k β = 0.58 for the 20 mm thick cover, it drops to k β = 0.32 for a 40 mm thick cover.The same trend was experimentally documented by Alonso et al. [10] who reported even six times larger β if current density was accelerated from 10 to 100 µA/cm 2 for the concrete cover of 70 mm.Also, it can be observed that the impact of curing time on k β is smaller than the impact of the thickness of the concrete cover.The knowledge of the accurate slope of crack width with respect to corrosion penetration in natural conditions is important for the accurate predictions of the length of the corrosion propagation period (time since corrosion initiation to critical corrosion-induced delamination/spalling of concrete cover).Currently, durability estimates in national design codes are based on time until corrosion initiation, which is dictated by the time needed for the diffusion of aggressive species such as chlorides from exposed concrete surfaces through concrete cover to steel rebars.This approach is overly conservative, as even though corrosion initiation is usually the longest part of the corrosion process, the length of the propagation period can be substantial.For example, the results for a 40 mm thick concrete cover and concrete cured for 28 days indicate that under natural conditions (typically i a ≤ 1 µA/cm 2 ) the propagation period would last for at least 7 years until the surface crack width of 0.3 mm (a commonly used Eurocode criterion of limit serviceability state) was reached.That said, even the crack opening of 0.3 mm can underestimate the time when and the material behaviour is characterized by Hooke's law for plane strain, The assumption of axial symmetry makes all considered quantities independent of ϕ, and also leads to u ϕ = 0. Therefore, from (A.3) follows that ε rϕ = 0 and (A.The derivation has been done for concrete considered as an elastic material characterized by Young's modulus E c and Poisson's ratio ν c .If the concrete is damaged, E c should be replaced by the secant modulus g(φ)E c where g is the degradation function and φ is the phase-field variable that describes the state of damage.In general, φ varies in the radial as well as circumferential direction, and the problem is no longer axially symmetric and does not have an analytical solution.
For simplicity, the effect of damage is considered only approximately by using a uniformly reduced modulus E c,d = g(φ(a))E c , where the degradation function g(φ) is evaluated on the inner concrete boundary, i.e., at r = a.This means that the damage of concrete is simplified to be uniformly distributed over the thick-walled cylinder and the value of the damage is the same as on the inner concrete boundary.The final formula for the compliance factor to be used in (A.16) is then .17)

= 2 •
(2)-(3), D II and D III are the second-order diffusivity tensors of Fe 2+ and Fe 3+ ions.R II ,R o and R h are the reaction rates of reactions (1b), (1c) and (1d), respectively, with k II→III r reaction constants.Since the ratio of the mass fractions of iron oxides and hydroxy-oxides in rust changes with the corrosion current density, reaction constants k II→o r and k III→h r have to depend on the corrosion current density too.Previous studies provided estimates of k II→III r = 0.1 mol −1 m 3 s −1 [31] and k III→h r 10 −4 s −1 [32], but k II→o r [s −1 ] was not known.For this reason, k III→h r is considered as fixed, and k II→o r is varied with

Fig. 5 :
Fig.5: Schematic illustration of rust precipitation in the vicinity of the steel rebar.A portion of the released Fe 2+ ions and emerging Fe 3+ ions precipitates to form a dense rust layer adjacent to the yet uncorroded steel surface while the remainder is transported further in concrete pore space where it eventually precipitates.
. The function H(t) stores the maximum reached value of the crack driving force H(t) in time.Crack nucleation occurs when H(t) exceeds the damage nucleation threshold H calculated from the basic properties of concrete.The function σ1 = (σ 1 + |σ 1 |)/2 is the positive part of the maximum principal value of the effective stress tensor σ = C c e : (ε − ε ⋆

Fig. 6 :
Fig. 6: Schematic illustration of the pressure induced by rust precipitation in the vicinity of the rebar.

Fig. 7 :
Fig. 7: Pressure p exerted by a dense rust layer on concrete -(a) Pressure p increases with the thickness

Fig. 8 :
Fig. 8: The impact of rust layer thickness and concrete rust saturation on the flux reduction coefficient k f and the impact of the magnitude of corrosion current on the distribution of rust in pore space -(a) Flux reduction coefficient k f ∈ [0, 1] reduces the Faraday law dictated flux J II,F ar , such that the flux released to concrete pore space is J II = k f J II,F ar .This reduction is the result of the deposition of a part of the precipitates in a dense rust layer adjacent to the corroding steel surface.The coefficient k f decreases with the thickness of the dense rust layer and with the saturation of the pore space by precipitates, such that J II eventually vanishes and all released Fe 2+ ions contribute to the formation of the dense rust layer.(b)

Fig. 9 :
Fig. 9: Growth of the upper surface crack width w depending on the thickness of the corroded steel layer t cor for three values of applied corrosion current density and for concrete cover (a) 20 mm and (b) 40 mm.For identical t cor , w decreases with the increasing applied current density as the result of the changing chemical composition of rust, causing its increasing density and thus decreased corrosion-induced pressure.The predicted values in the range where w ≥ 0.5 µm are linearly fitted by w = βt cor + ζ.While w(t cor ) is nearly linear for the smaller concrete cover of 20 mm, non-linear dependence is observed for the larger concrete cover of 40 mm.In (c), the typical fracture patterns characterised by the contours of phase-field variable φ are depicted for the simulated impressed current test conducted on 28 days old concrete with a cover of 20 mm and the corrosion current of i a = 1 µA/cm 2 .

Fig. 10 :
Fig.10: The slope β of linearly fitted predicted crack width (w = βt cor + c) is decreasing with the increasing applied corrosion current density.The reason is that the magnitude of applied current affects the composition of rust such that its density increases with the increasing applied current density, decreasing corrosion-induced pressure in the process.The comparison with experimental measurements (a) in log-log scale reveals that even though the decreasing trend of crack width slope with current density is captured well, experiments suggest even higher slopes for lower current densities such as i a = 5 µA/cm 2 .Figure (b) allows to more easily observe that crack width slope rapidly decreases with increasing current density up to approximately i a = 50 µA/cm 2 .From this, a more moderate decrease was observed.

Fig. 11 :
Fig. 11: (a) The pressure of the dense rust layer accumulating in the volume vacated by steel corrosion is responsible for a considerable part of corrosion-induced pressure but the pressure resulting from the constrained accumulation of rust in pore space is not negligible.(b) The correction factor k β represents the relative error caused by the current acceleration compared to the same test corroding under natural conditions.k β decreases with the applied corrosion current density and the thickness of the concrete cover.

2 )
in which the strains are linked to the displacements by
[29]umerical predictions agree quantitatively better with experimental measurements, though they mostly underestimate the experimental values.The comparison of data for the four combinations of two curing times (leading to different mechanical properties) and two concrete covers in Fig.?? indicate that β increases with increasing curing time (i.e.enhanced mechanical properties).Also, a larger cover can lead to larger β for small natural-like current densities (within the evaluated case studies, this held for the combination of 40 mm cover and 147 days curing time), as was previously observed in Ref.[29].The modelling results in Fig.?? indicate that the influence of concrete cover is significantly more important than of curing time (i.e.concrete strength).Also, it appears that the larger the current density, the smaller the impact of curing time (i.e.concrete strength) on the crack width slope.The pressure generated by the constrained accumulation of the dense rust layer has a considerable impact on corrosion-induced fracture (see Fig11a)