Modeling porosity loss in Fe0-based permeable reactive barriers with Faraday’s law

Solid iron corrosion products (FeCPs), continuously generated from iron corrosion in Fe0-based permeable reactive barriers (PRB) at pH > 4.5, can lead to significant porosity loss and possibility of system’s failure. To avoid such failure and to estimate the long-term performance of PRBs, reliable models are required. In this study, a mathematical model is presented to describe the porosity change of a hypothetical Fe0-based PRB through-flowed by deionized water. The porosity loss is solely caused by iron corrosion process. The new model is based on Faraday’s Law and considers the iron surface passivation. Experimental results from literature were used to calibrate the parameters of the model. The derived iron corrosion rates (2.60 mmol/(kg day), 2.07 mmol/(kg day) and 1.77 mmol/(kg day)) are significantly larger than the corrosion rate used in previous modeling studies (0.4 mmol/(kg day)). This suggests that the previous models have underestimated the impact of in-situ generated FeCPs on the porosity loss. The model results show that the assumptions for the iron corrosion rates on basis of a first-order dependency on iron surface area are only valid when no iron surface passivation is considered. The simulations demonstrate that volume-expansion by Fe0 corrosion products alone can cause a great extent of porosity loss and suggests careful evaluation of the iron corrosion process in individual Fe0-based PRB.

www.nature.com/scientificreports/ The model described the decline of iron reactivity as a function of space and time by observing the degradation of chlorinated ethanes, but did not consider iron corrosion processes. Numerous detailed geochemical models were proposed to simulate the chemical reactions and flow transport inside Fe 0 -PRB and the effect of mineral precipitation on hydraulic properties of PRBs 20,30,[41][42][43] . The iron corrosion rate in these models is expressed with a first-order dependence on iron surface area, and the rate coefficient was derived from the report of Reardon 44 . According to the modeling results, it is perceived that the permeability loss is mainly caused by foreign precipitates (e.g. CaCO 3 ) or mixed precipitates (e.g. with FeCO 3 ) 20 . Moreover, iron surface passivation is not considered in previous studies. However, above pH ~ 4.8 and in an oxygen containing aqueous environment the generated porous oxide layers contains three-valent iron causing strong inhibition effect on the iron corrosion 45 . In this study, Iron surface passivation is described as linear or parabolic growth of corrosion products 46,47 . The aging behavior due to corrosion of iron particles was investigated most of the time using column tests or by measuring the hydrogen pressure build-up in long-term batch studies 44,[48][49][50] . According to previous studies, the estimated life-time of iron granular particles ranged from several years to several decades 51,52 . Some studies assumed that iron particles will be completely consumed in the reaction of groundwater and estimated the lifetime of Fe 0 -based PRB by the iron mass and iron corrosion rate [53][54][55][56] . This approach to estimate the PRB service life is only valid, if the initial PRB pore volume and the used Fe 0 /aggregate ratio enable complete Fe 0 depletion 35,57 . In this study, the residual amount of Fe 0 and the residual pore space of the PRB are calculated on a time-line to evaluate this life-time estimation method. Therefore, in this study, a mathematical model is formulated to study the porosity loss of Fe 0 -based PRB solely caused by the volumetric expansive corrosion of iron based on Faraday's Law including iron surface passivation. For simplification, iron corrosion in the deionized (DI) water is considered. Based on the results of Luo et al. 58 , which show the porosity change of iron exposed to deionized water, our model is calibrated in order to simulate the porosity loss for long-term operation.

Fundamental of Fe 0 /H 2 O system
Since Fe 0 is not stable under environmental conditions, and the redox couple H + /H 2 (E 0 = 0.00 V) is higher than that of Fe II /Fe 0 (E 0 = − 0.44 V) at a H+ = 1 59,60 , a transfer of electrons from the Fe 0 body (solid state) to the Fe/H 2 O interface occurs whenever a Fe 0 specimen is immersed in an aqueous solution 59,61,62 . Equations (1a) and (1b) show that the oxidative dissolution of Fe 0 by protons (H + ) from water (H 2 O ⇔ H + + OH − ) forms Fe 2+ and Fe(OH) 2 by increasing the pH. In the presence of dissolved oxygen, Fe 2+ and Fe(OH) 2 can be oxidized to less soluble Fe(OH) 3 (Eqs. 2a, 2b). Fe(OH) 2 and Fe(OH) 3 are polymerized and further transformed to various oxyhydroxides (Eq. 3) 59,63-65 . Equation (4) summarizes the process of aqueous iron corrosion.
Comprehensive research on Fe 0 for water treatment revealed that the generation of iron oxyhydroxides (iron corrosion products or FeCPs) is the basis of contaminant removal in Fe 0 /H 2 O systems 7,8,60,66 . Figure 1 depicts the principle of contaminant removal process in Fe 0 /H 2 O system. The electrochemical corrosion of immersed Fe 0 induces the generation of reducing agents, i.e. Fe 2+ , Fe(OH) 2 and H 2 (Eqs. 1a, 1b). The generated iron oxyhydroxides on Fe 0 (red layer) is an adsorbent for contaminants 64 , as well as a contaminant scavenger (Eqs. 2a, 2b).
However, the generation of FeCPs has the side effect of being expansive. Thus, modeling the volume-expansion process is essential to design an efficient and sustainable Fe 0 /H 2 O filter system. where m is the mass of the substance liberated or deposited at the electrode (g); Q is the total electric charge (Coulombs or Amperes seconds)); M is the molar mass of the substance (g/mol), F is the Faraday constant; and z is the valency of an ion formed from the reacting substance. For Fe 0 oxidative dissolution electrochemical reaction (Eqs. 1a, 1b), Eq. (5) can be transformed by using Eqs. (6) and (7) into Eq. (8): where I is the current (Ampere), t the reaction time (s), i the current density (Ampere/m 2 ) and A the surface area of iron (m 2 ). ∂V iron is the volume depletion of iron, ρ the density of iron = 7.85 × 10 3 kg/m 3 , M = 55.85 g/ mol, F = 96,500 C/mol; z is taken equal to 2 (Eq. 1a).
Assuming the iron particle is a sphere, we obtain for one iron particle where ∂r iron is the radius depletion of the iron particle. Combining Eq. (8) and Eq. (9), we get: where ∂r iron ∂t is the corrosion rate (in mm/year) and M zFρ is a constant. So the corrosion rate (in mm/year) and the current density are linearly related.
Calculation of the coefficient of volumetric expansion. As discussed above, the generation of FeCPs is a volumetric-expansion process. A coefficient of volumetric expansion ( η ) is introduced to describe this behavior. Equation (11) states the definition of η: where V oxide is the volume of the generated FeCPs. The change of volume and radius can be described as follows: where V expansion is the expansion volume, r oxide is the increased radius with the generation of FeCPs and r expansion is the expansion radius of the iron particle.
If we combine Eq. (8) and Eq. (12), we have the expression of the total volume change ( V ) over time as: Due to the complexity of the iron corrosion product, η varies with different corrosion environment and Fe 0 intrinsic reactivity 65 . Table 1 depicts the volumetric expansion coefficients of different possible corrosion products based on the study of Caré et al. 33 .
Estimate of growth of corrosion products and passivation. The growth of corrosion products may follow a linear or parabolic law depending on the metal properties and geochemical conditions 46,47 . Due the complexity of iron corrosion, both linear and parabolic law are considered in this study.
For a metal that the generated oxide film is not protective, which means that no passivation occurs, the rate of growth of oxide film remains constant: where k is a constant, y is the film thickness (m) and t is the corrosion time.
For a metal that forms a protective oxide film, Q = It, www.nature.com/scientificreports/ However, the relationship between corrosion rate and time is not so simple. The following equation is usually used 68 : where n is the coefficient of passivation. The value of n is usually larger than 1. In case of iron corrosion in air or corrosion in soil, n ranges from 1 to 3 depending on the suppression of diffusion of oxygen through the formed oxide film 68 .
The coefficients of passivation n are taken equal to 1, 1.5 and 2 respectively in this study. For n = 1, the rate of growth of oxide film on the iron surface remains constant and the passivation of Fe 0 is not considered. For n = 2, it means that iron passivation occurs during the corrosion process and the generated corrosion products form a protective oxide film. For the case n = 1.5, iron passivation occurs with time but the oxide film on the surface of iron is not completely protective.
If we consider iron passivation occurs in the system, Eqs. (17) and (18)  where β is a constant.
Estimate of surface area and porosity change. The iron surface area is constantly changing during the operation. The area changes can be calculated by the following equations: where N is the total iron particle number, r 0 is the initial radium of iron particles. The total iron particle number can be calculated by Eq. (25): www.nature.com/scientificreports/ where V 0 solid is the initial volume occupied by the solid particles (i.e. iron and sand particles), and τ iron is the initial iron volume ratio.
The porosity of the system ( ) can be described as follow: where V reactive zone is the volume of the reactive zone, and i is the total volume change (Eq. 14).

Model assumptions.
A simplified illustration of the model is shown in Fig. 2. As iron constantly transforms in the water, the radius of the iron particle decreases with time. In the meanwhile, the generated FeCPs, which have larger volumes than the Fe 0 , fill the pore space and cause the porosity loss in the system. The following assumptions are made in this model.
• All iron particles are spheres and have identical radium, which is taken equal to 1 mm.
• Uniform Fe 0 corrosion: the radius reduction of spherical Fe 0 particles is the same for all particles.
• The volume of the reactive zone remains constant.
• Fe 0 corrosion products progressively fill the available pore space.
Model calibration. The results of Luo et al. 58 are used to calibrate the current density in this model. The corrosion rate (in mm/year) can be obtained with the calibrated current density value (Eq. 10). In Luo et al. 's study, iron/sand mixtures were packed between two layers of sand in the column experiments. Three iron mixing ratios (100%, 50%, 10%, w/w) of the barrier material were tested for different water type (a synthetic groundwater, acidic drainage and deionized (DI) water). The porosity data for only 100% Fe 0 columns which tested with DI water were reported. Since this study considers the condition that iron corrodes in DI water, the porosity change data for 100% Fe 0 column reacted with DI water were used to calibrate the parameters of the model. These parameters were then used to simulate the porosity loss in the systems with different iron mixing ratios.
The authors reported a porosity loss from 57 to 32% when 100% Fe 0 column are exposed to deionized water within 16 days under 2 mL/min flow velocity. The flow condition is 56 times greater than the average flow rate at a typical PRB installation 69 . In order to estimate the real operation of PRBs, the aging of the fast flow experiment have to be scaled to real time conditions. The surface-loading rate was used in the calculation. At this fast-flow rate, 16 days represents an equivalent reaction period of 2.2 years (800 days) under typical field conditions, as the surface-loading rate is proportional to flow velocity 48 . This prediction does not include the effect of kinetics on precipitation due to increased velocity and reduced residence time, and the operating life of PRBs will be overestimated using results from fast flow rate experiment 48,49 .
The grain size was reported as between 2.38 and 0.30 mm. Since a uniform particle radius is assumed in this study, the grain size taken was equal to 2 mm.
Model results. Corrosion rates for different coefficients of passivation. The current density (i) is the calibrated parameter in this study. The corrosion rate value can be then calculated (Eq. 10). The corrosion rate (in mm/year) is either a constant value or a function of the corrosion time when different coefficients of passivation (n) are taken. The derived corrosion rate values are shown in Fig. 3. It is assumed that goethite (FeOOH) is the only corrosion product, which is reported by Luo et al. 58 on the basis of SEM images and EDX spectra results.
The corrosion rate (in mm/year) is a constant (= 0.06 mm/year) if the passivation of iron corrosion is neglected. When the passivation of iron is considered ( n > 1 ), dramatic variations of the calibrated corrosion rate values are detected in the beginning phase of corrosion. The initial corrosion rate values are significantly Relative porosity loss for different coefficients of passivation. Figure 4 depicts the simulations of long-term porosity changes for different corrosion patterns. It is assumed that goethite is the only corrosion product. The Y-axis in the figure represents the relative porosity of the system, which can be described as where is the porosity at time t and 0 is the initial porosity of the system. Significant porosity losses can be detected in all three simulations (Fig. 4). The simulation with constant corrosion rate (in mm/year) shows the most remarkable porosity reduction, which decreases to 0 on day 2464. Zero relative porosity means there is no pore space left in the barrier, i.e. no underground water can flow through the PRB. Thus, the PRB has no water remediation effect after day 2464.
The relative porosity values of three simulations show dissimilar features along corrosion pathway and also divergent results after long-term simulation. After 10 years of simulation, the relative porosity values decrease to 5.94% for coefficient of passivation of n = 1.5 and 19.50% for n = 2. The simulation with higher coefficient of passivation (n) shows a more rapid porosity loss in the beginning phase but only slight porosity change after long-term corrosion. This indicates that the rate of diffusion process decreases with the increase of the thickness of the generated corrosion products. The differences among the three simulation results indicate that the iron passivation is an important factor determining porosity for Fe 0 -based PRBs' long-term performance estimation.  www.nature.com/scientificreports/ Relative porosity loss for different iron mixing ratios. The calculated corrosion rates were utilized to simulate the porosity loss in the systems with different iron mixing ratios. Figure 5 depicts the relative porosity change along time with iron mixing ratios of 10%, 50% and 100% (W/W). It can be seen that a lower percentage of Fe 0 within the barrier shows less porosity reduction during long-term operation.
Relative porosity loss for different corrosion products. The previous simulations in this study all assume that goethite is the only corrosion product of iron corrosion. However, in the real corrosion process, other FeCPs may form. Figure 6 illustrates the porosity loss simulations for different possible corrosion products with no iron passivation considered. All possible iron corrosion products have a larger volume than the corroded iron. In general, the simulations with higher coefficients of volumetric expansion ( η ) show stronger porosity reduction. The results of the simulations imply that iron corrosion products have an important effect on the porosity reduction of the PRB system. Figure 7 shows the calculated porosity and Fe 0 volume decrease with time by the formation of goethite. According to simulation results, when the porosity value of the simulation with constant corrosion rate (in mm/year) reaches 0, there is still 0.09 m 3 /m 3 Fe 0 volume fraction left in the system. It means the PRB system loses its capability to remove contaminants before the iron is completely consumed. Therefore, the previous method to estimate the lifetime of Fe 0 on basis of the corrosion rate 56 , which assumes the iron will be totally oxidized, cannot be used to estimate the service lifetime of iron-based PRB systems. Moreover, this study simulates only the contact of Fe 0 and deionized water and considers merely the effect of expansive volume of iron corrosion. If the geochemical condition changes, e.g. the solution has high calcium concentration, which will cause additional mineral precipitation in the iron zone and trigger even larger porosity loss, an earlier failure of PRB technique can be expected.

Comparison of corrosion rates in different studies
It is difficult to compare the corrosion rate results of this study with those of previous studies, since different units of corrosion rate were used (e.g. mmol/(kg day), mol/(m 2 day), etc.) in former studies. The most frequently used corrosion rate value in Fe 0 -based PRB modeling studies 20,30,41,50 is derived from the report of Reardon et al. 44 . Reardon measured the iron corrosion rates by monitoring the hydrogen pressure increase in sealed cells containing iron granules and water. The corrosion rate was given in mmol/(kg day). The units of corrosion rates in this study are converted to mmol/(kg day) and the results are shown in Fig. 8a,b. Differences of more than an order of magnitude are shown between the calibrated corrosion rates in this study and the data reported by Reardon 44 at the beginning of corrosion. The initial corrosion rates with different coefficients of passivation (n) in this study are 5.0 mmol/(kg day) (n = 1), 31.4 mmol/(kg day) (n = 1.5) and 73.2 mmol/ (kg day) (n = 2) respectively. The corrosion rates strongly decrease over time and the average corrosion rates  www.nature.com/scientificreports/ of a 10 years simulation are 2.60 mmol/(kg day) (n = 1), 2.07 mmol/(kg day) (n = 1.5), and 1.77 mmol/(kg day) (n = 2) respectively. These rates are higher than the rate of 0.4 mmol/(kg day) for deionized water published by Reardon 44 but fall within the range of reported corrosion rates of 0.2-50 mmol/(kg day) 49 .
The relatively low corrosion rates reported by Reardon 1995 contradicts to the porosity reduction (57-32%) observed by Luo et al. 58 . In both cases, deionized water was in contact with granular iron. The possible reasons for this remarkable difference are (1) Reardon measured the corrosion rates with a batch experiment with no water flow (rate = 0). But in typical PRB systems, the underground water flows through the barrier under the natural hydraulic gradient, which increases the iron corrosion rate 70 . (2) Reardon monitored the hydrogen pressure increase in sealed cells. The increasing hydrogen partial pressure may inhibit the iron corrosion reaction. For a real site PRB, the generated hydrogen can escape immediately from the system. (3) The intrinsic reactivity of iron materials varies significantly which may cause an order magnitude difference in corrosion rate 71 .
Model approaches on the basis of data from Reardon 44 underestimate the iron corrosion process as well as the influence of the iron corrosion products on the porosity of the PRB system during the long-term operation. With higher corrosion rates of iron, more iron will react in the water, and larger amount of iron corrosion products are generated. The increasing generated corrosion products can fill the pore in the barrier quickly and cause the early failure of the PRB technique.
Moreover, the reaction process is very complicated in the real PRB systems. The iron corrosion rate is easily influenced by many factors 20 and can vary dramatically in different parts of the PRB. For example, dissolved oxygen (DO) is consumed once it enter the iron zone, and it can accelerate the iron corrosion process, which induces more corrosion products, and thus higher porosity reduction in the entrance zone of PRBs 27 .

Considerations on reactive surface area change versus time
The reactive surface area of Fe 0 in this study can be calculated by the depletion of the radius of the iron particles (Eq. 10) and the assumption of uniform corrosion. Plots of calculated reactive surface area of different coefficients of passivation (n) over time are shown in Fig. 9.
For most of previous modeling studies on simulating the operation of Fe 0 -based PRBs, the reaction rate for iron corrosion by water was assumed to have a first-order dependency on the iron surface area 20,30,41,42,50,[72][73][74][75][76] . A plot of calibrated corrosion rates versus the calculated reactive surface area values under different coefficients of passivation is shown in Fig. 10.
When no iron passivation considered, which means that the current density (i) or the corrosion rate (in mm/year) are constant, the corrosion rate in mmol/(kg day) has a complete first-order dependency on the iron  www.nature.com/scientificreports/ surface area (Fig. 10). However, when the growth of corrosion products follows a parabolic law, the first-order dependence assumption is no longer applicable. The first order model underestimates the iron corrosion rate at large reactive surface area in the system. As shown in Fig. 10, the corrosion rate results for coefficient of passivation (n) is 1.5 or 2 are significantly larger than that of first-order dependency (n = 1) when the system has a large reactive surface area. This is valid especially during the beginning phase with strong increase of iron corrosion products, which causes the porosity reduction in the barrier. Thus, a more accurate expression of corrosion rate should be applied in modeling simulations in order to have a better estimation of PRB endurance. Table 2 summarized the porosity loss simulations of previous studies with the simulations in this study with assumptions that the corrosion rate (in mm/year) is a constant and goethite is the only corrosion product. The porosity loss after 1 year simulation in this study is over one order of magnitude larger than the simulation results from former studies. The simulated porosity loss from these studies contradict to the measured porosity loss reported by Luo et al. 58 . Possible reasons for the divergence are the different iron corrosion rates utilized in the model as discussed in "Comparison of corrosion rates in different studies". In addition, this study simulated the condition that the granular iron is in contact with deionized water and only the iron corrosion process is considered. If the deionized water is replaced by underground water with multiple dissolved ions, the porosity loss, i.e. by precipitation of carbonates, can be even more significant. Figure 5 shows the simulated long-term porosity loss for systems with different Fe 0 mixing ratios. Table 3  The results from column experiments in Luo et al. 58 confirm that the system with lower Fe 0 mixing ratio in the barrier can remain higher porosity after exposure to water. The higher porosity can be explained that a lower percentage of Fe 0 generates less corrosion products, which reduce the likelihood of pore clogging in the system. Therefore, mixing Fe 0 and less reactive materials (e.g. sand) is a solution for long-term porosity loss in Fe 0 -based PRBs 77 . However, a low Fe 0 mixing ratio might reduce the ability of a PRB system to remove contaminants 58 . Thus an appropriate ratio between Fe 0 and less reactive materials is important.

Conclusion
A mathematical model is presented to simulate the long-term porosity loss of Fe 0 -based PRBs as induced by deionized water. It is assumed that only the volumetric expansive corrosion of iron contributes to the porosity loss of the system. Faraday's law was applied to describe the correlation of the amount of corroded iron and the iron corrosion rates. Different coefficients of passivation were taken into account to describe different growth features of corrosion products. Measured porosity results from Luo et al. 58 were used to calibrate the parameters in the model. Based on experimental findings from literature and the simulations here, the following major conclusions can be dawn.
(a) There are iron residues in the system (0.09 m 3 Fe 0 /m 3 ) when the porosity reduces to 0, which means the groundwater can no longer flow through the Fe 0 -based PRB before the Fe 0 is completely consumed. Thus, it is not correct to assume that the iron in Fe 0 -based PRB is totally consumed and that the endurance of PRB can be estimated from the amount of iron and iron corrosion rate. (b) The derived iron corrosion rates in presented model (2.60 mmol/(kg day), 2.07 mmol/(kg day) and 1.77 mmol/(kg day)) are significantly larger than the corrosion rate used in previous studies (0.4 mmol/  www.nature.com/scientificreports/ (kg day)). Higher iron corrosion rate means more iron can dissolve in the water, which leads to more significant porosity loss caused by larger amount of generated iron corrosion products. Thus, the previous simulations with low iron corrosion rate may underestimate the porosity loss in PRB. Moreover, we propose, a uniform unit of iron corrosion rate (e.g. mm/year) for Fe 0 -based PRB systems in order to improve the comparability of the different studies. (c) The assumption in previous modeling studies, which describes the iron corrosion rate (in mmol/(kg day)) as a first-order dependency on iron surface area, is accurate only when iron passivation is neglected. When iron passivation is considered, such an assumption underestimates the corrosion rates especially at the beginning phase of operation. (d) The modelled porosity loss in this study (0.12/year with assumptions that the corrosion rate is a constant and goethite is the only corrosion product) is larger than the simulation results from previous studies (average 0.02/year). Our study demonstrates that iron corrosion products can cause large porosity loss in the filter. Iron passivation features and possible corrosion products are responsible for large differences between the simulation results. Therefore, iron corrosion processes need to be properly considered in order to accurately estimate the long-term operation of Fe 0 -based PRB systems.