A multi-model approach toward understanding iron fouling at rock-fill drainage sites along roadways in New Hampshire, USA

Factors affecting iron fouling in wet areas adjacent to roadways were investigated by collecting field rock cut and aqueous physicochemical data; developing exploratory predictive models; and developing geochemical models. Basic data included the identification of iron fouling from aerial imagery and field visits at 374 New Hampshire rock cut locations, and their associated rock-fill sites. Based on field water quality measurements from wet areas at 36 of the rock-fill sites, the occurrence of iron fouling was associated with higher values of specific conductance, lower concentrations of dissolved oxygen and lower pH compared to areas without iron fouling. A statistical model, using boosted regression trees, was developed to predict the occurrence of iron fouling in wet areas adjacent to roadways where rock-fill from nearby rock cuts was used in roadway construction. The model was used to develop a continuous iron fouling probability map for the state of New Hampshire that can be used to better understand the occurrence of iron fouling. Geochemical models illustrate how iron fouling of waters increases along roadways built with fill from sulfidic rock cuts as a result of acid generation from pyrite dissolution and ferrous iron (Fe2+) oxidation and increases in areas with greater specific conductance from deicing runoff caused by cation exchange. More iron is precipitated as goethite in simulations that include pyrite, and in simulations with deicing salts added, indicating that rock-fill sites with rocks that contain pyrite and water with greater salt content could have enhanced iron fouling.


Introduction
Metasedimentary rocks with iron-rich biotite mica and (or) iron sulfide minerals, such as pyrite, are potential sources of iron in New Hampshire (NH) groundwater [1]. When rock cuts are made during road construction the excavated rock is commonly used as rock-fill in the nearby roadway base and along roadways to channel and drain storm water [2,3]. These freshly excavated rocks may contain iron-rich minerals, which when exposed to oxygen, water, and bacteria, become weathered. Subsequently, deposits of red-orange iron oxide minerals and biofilms may form (collectively referred to as iron fouling). Naturally occurring bacteria commonly catalyze iron reactions and form a biofilm. The iron-associated water quality changes and microbial deposits can result in adverse impacts to aquatic organisms, water bodies, streambeds, and roadway structures. The NH Department of Transportation (NHDOT) has identified iron fouling as a continuous maintenance challenge in stream crossings, culvert outlets, and other areas [4]. In this study, we develop a statistical model to predict iron fouling occurrence from road fill throughout NH and use mechanistic geochemical models to understand the geochemical conditions contributing to iron fouling.
Statistical models are commonly used to predict environmental conditions including groundwater and surface-water quality. Machine-learning methods such as boosted regression trees (BRT) are well suited to prediction, fit complex, non-linear relations, handle different types of predictor variables (numeric, binary, categorical, etc.) and do not require data transformation or elimination of outliers to satisfy assumptions of traditional statistical methods [5]. BRTs combine regression trees, which explain the variation of a single response variable by repeatedly splitting the data into more homogeneous groups; and boosting, a machine-learning technique that combines many simple models to get improved predictive performance [5,6]. These methods have been used to develop predictive models for a range of environmental applications including groundwater quality [7][8][9][10], roadside air quality [11], landslide susceptibility [12], and ecological studies [5]. Studies comparing models developed with parametric techniques to those developed with machine-learning techniques have found that the machine learning tends to produce models with better predictive capabilities [8,13]. This study is a new application of BRT models.
Acid rock drainage (ARD) is common at mine drainage sites [14,15] and is often caused by the oxidation of iron sulfide minerals such as pyrite, and the subsequent release of acidity, sulfate, and dissolved metals including iron [16][17][18][19]. Iron fouling in runoff from rock cuts [20][21][22] and adjacent to blasted rock-fill used for roadway construction [23][24][25] is caused by a similar geochemical mechanism to the one causing ARD at mine drainage sites [26]. Variations in temperature, dissolved oxygen (DO), pH, and specific conductance (SC) values in surface water can also affect the concentrations of iron [27] and subsequently, iron fouling. Higher salinity water can increase iron mobility [28,29] and aggravate ARD at these sites.
The oxidation of 1 mol of pyrite by O 2 produces 1 mol of Fe 2+ , 2 mol of SO 4 2and 2 mol of H + [30]: The Fe(II) initially released may be oxidized to Fe(III), which is a more effective oxidant than oxygen [31][32][33], but the reaction is limited by the low aqueous Fe 3+ (1) concentrations at neutral pH. Further, the oxidation rates change by orders of magnitude depending on the pH and other factors such as microbial catalysis [17,33,34]. Precipitates of Fe(III) oxyhydroxides may form, such as schwertmannite, ferrihydrite, or goethite, and can be colloidal [35,36]. The oxidation of Fe 2+ initially causes a rise in pH because the oxidation of Fe(II) consumes protons, but then the Fe 3+ hydrolyzes and decreases the pH, and finally, the hydrolyzing ion precipitates as a mineral phase and further decreases the pH [17]: The overall reaction of pyrite oxidation coupled to the precipitation of iron oxyhydroxides can be described as: and results in the release of 4 mol of H + for each mole of pyrite oxidized. The resulting low-pH conditions can have adverse effects on aquatic life [37,38]. The pH and chemistry of ARD waters are commonly buffered by carbonates, and iron solubility can be limited by precipitates of Al, Fe, Mn, Ba, and Sr, and by SO 4 complexation and sorption reactions [16,39]. Furthermore, the effect of cations such as Na + (aq) exchanging for Fe 2+ in exchangeable cations represented by FeX 2(ex) (Eq. 4) and for H + in HX (ex) (Eq. 5) can result in an increase in dissolved iron concentrations and a decrease in pH [28,40,41] and exacerbate iron fouling.
Cations such as Na + can compete with Fe 2+ for sorption or cation exchange sites on goethite and other iron minerals [28]; changes in pH, redox, temperature, biological activity, and photoreduction can affect the capacity of iron minerals for sorption or exchange [27,[42][43][44].
Geology and geochemistry in NH at the state scale are associated with high iron concentrations in groundwater [53]. However, iron has many mineralogical sources and mobilizing conditions [54]. For example, iron can be mobilized under oxic, low pH conditions typical of ARD sites, or under anoxic (reducing) conditions with a neutral or higher pH [55][56][57]. Sulfides including pyrite commonly occur as secondary minerals in NH bedrock along fractures, joints, faults, or in vuggy zones [58][59][60]. Reducing conditions are common in older groundwater [53], which can dissolve and transport iron from aquifers, rock-fill, and stream bed material before becoming oxygenated in surface water and precipitating iron. Also, independent redox microzones can form within the hyporheic zone of the streambed [61] providing an alternative source of reducing water. Iron fouling affects water bodies and highway structures and is a concern for many transportation related jurisdictions [18,20,26,62]. Although iron geochemistry is well understood, investigations of iron fouling related to rock cuts and rock-fill usage at roadway sites in NH are rare and are mostly qualitative. A study of iron fouling along roadways constructed with rock-fill in NH included a site in Meredith, where rock-fill consisting of biotite schist was a likely iron source to drainage water [4]. Iron fouling occurs at other sites across the state, but no systematic research has been undertaken. NHDOT has indicated that more quantitative research is needed to understand the occurrence and to help predict and mitigate iron fouling. A probability map is developed here to predict the likelihood of iron fouling at untested locations throughout NH. In addition, geochemical models are presented to better understand the relations between high SC and iron fouling observed in the field.

Iron fouling data collection
A database of large roadcut locations maintained by the NHDOT [63] identified 374 rock cuts that are well distributed along the roadways throughout NH (Fig. 1). U.S. Geological Survey (USGS) National Hydrography Dataset (NHD) flowlines and waterbodies [64] were used to identify adjacent or nearby areas where rock-fill associated with each rock cut could potentially contact water and have the potential for iron fouling. If there were multiple rockfill sites near a rock cut, the closest rock-fill site was paired with the rock cut. Data used in this study were collected at both the rock-fill site and the associated rock cut site. Data from the rock cut sites included the presence or absence and color of rock weathering or staining (if present) and the identification of the bedrock for the statistical model. Data from the road fill locations included the presence or absence of iron fouling, all other variables used in the statistical model, and field measurements.
Imagery analyses using images from Google Earth (both aerial view and street views) and Esri ArcGIS World Imagery Map Services were conducted at all 374 sites to identify evidence of iron fouling at the rock-fill sites and weathering and staining at rock cuts. Leaf-off imagery provides a better view of the land surface beneath forest canopies at and near rock-fill sites and was available at 60% of the sites. A rock-fill site was identified as having iron fouling if the images indicated orange and/or reddish colors at streambeds, retention ponds, culverts, or wetlands near roadways. The imagery was also used to identify rocks that had minerals weathering in place to a red and brown color at the associated rock cuts, and the presence and color of staining of the rocks Fig. 1 Presence or absence of iron fouling at 374 rock-fill sites inNH from intermittent seepage of groundwater discharge at or below any visible fractures in the rock cut. "Rustyweathering" serves as an indicator of sulfidic-bearing bedrock units [65].
Field visits to 206 of the 374 rock-fill sites were conducted in 2016 and 2017 to look for additional evidence of iron fouling or to confirm evidence of iron fouling documented from the imagery analysis (Table 1). Sites visited were chosen to provide a wide spatial distribution and did not depend on the presence or absence of iron fouling during the imagery analysis. During field visits, retention ponds, culverts or wetlands near roadways with fill locations were inspected for evidence of iron fouling. A site that exhibited iron fouling in an aerial photo may not have exhibited fouling at the time of the field visit due to seasonal differences or changes over time at the site. Conversely, iron fouling was sometimes seen during the field visit at sites that did not show iron fouling in the aerial photos. The temporal variability in the presence of iron fouling or the utility of aerial photos for the identification of iron fouling due to contrasts with ground cover and foliage, weather, and forest canopy leaf-on or leaf-off conditions was not assessed in this study. Thus, evidence of iron fouling during either the desktop reconnaissance or the field reconnaissance resulted in a positive identification of iron fouling.
Physicochemical data including pH, DO, SC and temperature were measured in situ in surface water in contact with rock-fill at 36 of the 374 rock-fill sites in October and November 2016 [66]. DO measurements were used primarily to distinguish oxic from anoxic conditions that affect iron mobility. Sites with iron fouling and where water was discharging from springs or seeps downgradient of rock-fill were targeted preferentially for water quality measurements. The sites generally had small drainage areas, and many were surrounded by forested land. Physicochemical data were quality assured and approved in accordance with USGS protocols [67].

Statistical modeling
A boosted regression tree (BRT) statistical model was developed to predict the probability of iron fouling occurring at any location within NH. The model was fitted to the rock-fill dataset using the caret [68] and generalized boosted modules (gbm) [69] packages in the R computing environment [70].
The dependent variable in the model was the presence or absence of iron fouling at the 374 rock-fill locations as determined from desktop and field reconnaissance. Many potential predictor variables were assembled and evaluated for use in the final model. The variables were primarily derived from Geographic Information System (GIS) datasets and varied in scale from 1:24,000 to 1:500,000 (Supplementary Materials Table S1). In the model development, 111 individual predictor variables were tested and describe the geology, stream sediment and water chemistry, hydrology, and topography of the rock-fill site. Variables evaluated include 24 variables characterizing the lithology and lithogeochemical character of near-surface bedrock in NH [71]; stream sediment concentrations of 7 major oxides and 10 trace elements [72]; 3 surface-water chemistry variables [72], and the bedrock geologic unit at each site [73]. Predictor variables that are categorical, such as bedrock geologic unit, were converted to one variable for each geologic unit. For example, there are 50 bedrock geologic units represented in the model training dataset and each bedrock unit is a variable. The presence or absence of that bedrock unit at a location is represented as a binary value. Hydrologic and topographic variables including mean annual precipitation [74], predicted bedrock well yield [75], and elevation and slope characteristics obtained from 1:24,000 and 1:250,000 digital elevation models (DEMs) [76] were also tested. Each data source contains inherent uncertainties and sources of error; however these were not assessed in regard to the model outcome. Other potential predictor variables were identified but excluded in the final model, and these include soil chemistry information [77] and the soil survey geographic database for NH (SSURGO) [78]. The soil chemistry data were excluded due to the course resolution of the data (approximately 10 km 2 grids) and the SSURGO data were excluded due to large areas of missing data for NH. The 374 rock-fill sites were divided into a model training dataset containing 262 sites (70%) and a testing dataset containing 112 sites (30%). The model training dataset was used to tune the model using tenfold cross-validation methods while the testing dataset was used to test the model prediction capabilities of the most accurate model from model tuning and select a simpler model within one standard error of the most accurate model (1SE model). The testing dataset was also used to evaluate the prediction capabilities of the 1SE model, however this cannot be considered an independent validation because the same dataset was used to select the model. Model prediction results for the training and testing datasets are a probability value ranging from 0 to 1 of the likelihood of iron fouling. A desktop reconnaissance or field reconnaissance positive identification of iron fouling resulted in an iron fouling presence designation in the model. A probability cut point of 0.5 was used to assign a presence or absence of iron fouling, i.e., probabilities greater than or equal to 0.5 are assigned a value of 1 indicating the presence of iron fouling, probabilities less than 0.5 are assigned a value of 0 indicating the absence of iron fouling.
The model was tuned using the caret package [68] in the R computing environment [70]. Metaparameters of the model adjusted during tuning were number of trees (number of iterations), interaction depth (maximum nodes per tree, where 2 means no variable interaction), shrinkage (learning rate), and minobsinnode (minimum number of observations in each terminal node of a tree). The model tuning grid consisted of 640 combinations of the metaparameters. The most accurate model from tenfold cross-validation tuning in caret was determined. Machine-learning models tend to overpredict to model training data and less complex models typically have better prediction metrics when applied to new data [8]. Less complex models that were within one standard error of the most accurate tuning model (1SE models) were tested on the testing data to determine if better model prediction was possible. Following the selection of a 1SE model, the model was run through a variable reduction loop that sequentially removes the variable of least influence and tests the predictive capability of the model on testing data with the remaining variables. The number of variables in the final model was chosen with the goal of creating a simpler model that does not result in a decrease in predictive performance of the model.
Metrics used to evaluate model prediction were total accuracy, sensitivity, specificity, Kappa, and area under the receiver operating characteristic curve (ROC) ( Table 2). Total accuracy is the ratio of correct model predictions to known values. Sensitivity is the ratio of correct predictions of iron fouling to known locations of iron fouling, or true positives. Locations with predicted probabilities of iron fouling greater than 0.5 are assigned to the presence of iron fouling. Specificity is the ratio of correct predictions of no iron fouling to known locations without iron fouling (probabilities < 0.5), or true negatives. The Kappa statistic is also a measure of agreement between model predictions and observations but includes expected accuracy under chance agreement. Kappa values range from 0 to 1, with a value of 1 indicating complete agreement [79]. The area under the ROC is an indicator of model prediction and considers using all possible cut-points (not only 0.5) to compare predictions with observations. Values of ROC range from 0 to 1 with 1 indicating total agreement [80].
After the selection of a final model a probability map that covers the state of NH was produced. The map indicates the probability of iron fouling occurring at any location in the state if bedrock from that location were quarried, used as road fill, and then came into contact with water. To create this map from the model, each predictor variable was rendered into a GIS raster of 30 m resolution, exported as an ASCII file, and read into R using the raster package [81]. The predict function was used to compute the model predictions and the result was converted to a raster file.

Geochemical modeling
Geochemical models using the PHREEQC (pH Redox Equilibrium in C language) [82,83] a computer program for simulating chemical reactions and transport processes, were used to investigate the hypothesis that increased iron fouling occurs in areas with greater SC of drainage water (road-salt deicer application is the assumed source of elevated SC). Hypothetical models, based on data collected during this and other related studies, were set up to simulate the effects of pyrite oxidation and cation exchange processes caused by runoff with high salinity flowing along and/or infiltrating rock-fill generated from sulfidic rock cuts in NH. The models include a hypothetical starting groundwater solution (GW) and a hypothetical road-salt deicer effluent solution (D1) as described in Sect. 3.4 and Table 4. GW was mixed with D1 in a 70:30 ratio to simulate road deicer applications. Phases used in the model include the minerals pyrite, calcite, and biotite, the gases O 2 (gas) and CO 2 (gas), and an exchange composition of 250 mM.

Physicochemical results
Iron fouling was identified at 71 out of the 374 (19.0%) rock-fill sites [66] ( Table 1). The 71 sites included 47 sites that had evidence of iron fouling during the imagery analysis (31 of which also had evidence of iron fouling during the field visit), and an additional 25 that had evidence of iron fouling during the field visit but did not have evidence of iron fouling on the imagery, thus additional field visits may have increased the number of iron fouling observations (Table 1). A total of 303 sites did not have evidence of iron fouling in either the imagery analysis or the field visit.
Of the 374 rock cut sites, 153 had some red and brown weathering, and all but 2 also had iron staining from groundwater discharge below fractures in the rock cut. An additional 157 sites had staining below fractures but did not consist of rock that had red or brown weathering of minerals in place. Thus, 308 sites had staining of one or more colors at or below fractures including 289 rock cuts with red and brown staining, 202 with black staining, 23 with yellow staining and 27 with white staining. Many more of the rock cut sites (310) had evidence of iron (rock weathering, staining from discharge or both) than rock-fill sites that had evidence of iron fouling (71), suggesting that more sites have iron than end up being iron fouled.
There was quite a bit of variability in water quality properties measured in springs, seeps, and surface water adjacent and downgradient of 36 rock-fill sites. Values of pH ranged from 3.3 to 7.4; DO ranged from 0.3 to 11.9 mg/L; temperature ranged from 3.5 to 13 °C; and SC ranged from 17 to 3560 µS/cm. Water quality properties were compared to the presence or absence of iron fouling at these rock-fill sites (Fig. 2). The two sample Wilcoxon rank sum test [84] indicates that SC is higher at sites where iron fouling was observed (median value of 837 µS/cm) than at sites where it was not observed (median = 168 µS/cm, p value = 0.012) (Fig. 2a). Iron fouling was only observed at sites that had SC greater than 320 µS/cm. DO and pH were both lower at sites with observed iron fouling (p-values of 0.036 and 0.059 respectively, as calculated with the Wilcoxon rank sum test) (Fig. 2b, c). Further, there was a tendency toward anoxic water at the sites with iron fouling as indicated by the low values of DO in the lower 50th percentile of the data (Fig. 2c).
Analyses of rock weathering and staining observed from aerial and street view photos indicate that there is an increased chance of iron fouling at rock-fill sites associated with rock cuts that have red and brown weathering (p-value = 0.01) or red and brown staining below fractures line (p-value = 0.02). There is an increased probability of iron fouling at rock-fill sites associated with rock cuts that have yellow staining (p-value = 0.07) although the statistical significance is less pronounced likely due to the relatively small number of sites with yellow staining (n = 23). Sites with black and white staining did not have a statistically significant increased presence of iron fouling.

Statistical modeling
A model that was within 1SE of the most accurate model from the cross-validation tuning was selected as the final model (Table 2). In this study, the selected model had better prediction metrics for the testing data compared to the most accurate model from cross-validation tuning.
The number of predictor variables used in the final model was reduced from 111 to 30 by sequentially removing the variable of least influence and testing the predictive capability of the model on testing data with the remaining variables. For this dataset the 1SE model with 30 variables predicting to testing data had slightly higher accuracy, specificity, Kappa, and ROC values than the model with all 111 variables ( Table 2).
Seventeen of the 30 variables in the final model include stream sediment concentrations of metal oxides and trace metals, and steam water pH, alkalinity, and SC from Robinson et al. [72]. These stream sediment and stream water values were interpolated between sampling points for a wall to wall coverage of NH [72]. Ten of the variables in the final model are elevation and topographic characteristics derived from DEMs at the 1:24,000 and 1:250,000 scale [76]. The three remaining variables in the model are the mean annual precipitation, occurrence of metamorphic bedrock, and the modeled probability of obtaining 40 gallons of water per minute from a 400-foot-deep bedrock well. The model variables and their relative influence are presented in Table 3.
The field measurements of pH and SC do not correspond to interpolated values from stream water used to make the statistical model. The SC values from Robinson et al. [72] are generally much lower than the field measurements. The highest value for SC from Robinson et al. [72] that corresponds to field measurement locations is 208 µS/cm. Only four of the 36 field locations had measured SC values below 208 µS/cm. The pH values also vary between the values used to make the model and those measured at field locations with no discernible pattern between the values. These differences are not surprising as the values used in the model are generally representative of background conditions, were not collected to identify hot spots, and were interpolated over large areas. Thus, the values used to make the model will not necessarily agree with measurements at specific locations and, in fact, were generally inversely related.
Partial dependence plots for the variables in the final model provide insights into the relationships between each variable and the model probability outcomes [85]. The aluminum oxide concentration in stream sediments has the largest relative influence on the final BRT model (8.97%), followed by the silicon dioxide concentration in stream sediments (8.36%) ( Table 3). The partial dependence plots for these two variables show a similar pattern with respect to model predictions where higher concentrations are associated with a higher probability of iron fouling (Fig. 3).
The remaining 28 variables have between 7.44 and 1.18% relative influence on the model ( Table 2). The partial dependence plots for these variables are available in the Supplementary Materials Figures S1A-S1AB. General patterns shown by these plots indicate that the probability of iron fouling increases with increases in stream sediment concentrations of lead, barium, arsenic, strontium, and mercury (Figures S1A, S1G, S1M, S1V, and S1X). The probability of iron fouling also increases with increases in surface elevation, mean annual precipitation, and the presence of metamorphic bedrock (Figures S1D, S1AA, S1F, and S1T). The probability of iron fouling increases with decreases in stream water pH and SC and stream sediment concentrations of zinc ( Figures S1E, S1I, S1W). The planform curvature variables at the local (1:24,000) and regional (1:250,000) scales have higher probabilities of iron   [72] fouling as they become more convex (Figures S1H, S1O). The local angles of slope variables have different relationships with the probability of iron fouling at the regional and local scales. At the regional scale, the probability of iron fouling is greatest at the lowest slope values ( Figure  S1R), whereas at the local scale, the probability of iron fouling is greatest at the highest slope values ( Figure S1Z). The remaining variables have more complex, non-linear relationships with the probability of iron fouling.

Probability mapping
A continuous iron-fouling probability map covering NH was developed using the BRT model (Fig. 4). The map shows the probability of iron fouling occurring at any location across the state where rock-fill used in road construction projects comes into contact with drainage. Locations that have a greater than 50% chance of iron fouling are shown in red (Fig. 4). A digital raster file of Fig. 4, which can be used with a GIS, is available as a part of an associated data release [66].

Geochemical modeling
Water in contact with sulfidic rock-fill along highways can cause ARD given oxic conditions. Geochemical simulations with PHREEQC were used to model these conditions under several hypothetical scenarios informed by data collected mg/kg milligrams per kilogram for this study. The oxidative dissolution of pyrite causes increased acidity, and under oxic conditions, the oxidation of ferrous iron, and precipitation of Fe (hydr)oxides such as goethite. Field data from this study suggest that areas with locally high SC and low pH which could result from pyrite oxidation and deicer salts, have higher incidences of iron fouling. SC and pH are also variables in the statistical model. We focus on the pH and SC variables in the geochemical model because they are easily measured in the field. With increased deicer concentrations, elevated dissolved Na + and other cations replace exchangeable cations (including FeX 2 and HX), resulting in decreased pH, and increased dissolved Fe(II), which oxidizes to Fe(III). The dissolved Fe(III) precipitates out as the water is buffered and pH rises, as simulated using calcite as an equilibrium phase. Carbonate and sulfide minerals are highly reactive and have a disproportionately large effect on water chemistry compared to other bedrock minerals in this region [65]. Highway runoff from road deicers commonly contains high concentrations of Na, Ca, and other ions that can cause cation exchange [47]. The models include a starting groundwater solution (GW) based on median major constituent concentrations for 28 samples collected from supply wells drilled in calcareous metasedimentary rocks in NH and Maine as a part of the USGS National Water-Quality Assessment (NAWQA) New England Coastal Basins (NECB) Major Aquifer Study (lithogeochemical group M c in [7]; Table 4). These wells were used because it was noted that rock-fill locations in this type of rock tended to have iron fouling. This lithogeochemical group is primarily located in the southeastern part of NH and the southern part of Maine. The groundwater is assumed to be oxic, with a DO concentration of 7 mg/L (or 0.44 mM (mM)) and low dissolved iron concentrations (0.2 mg/L or 0.002 mM). The models used an initial exchange composition of 250 mM, but then the exchange sites were equilibrated with the GW solution or the mixture of GW with highway deicing runoff (D1).
The primary model demonstrates the following:  (1) Initial increase in aqueous Fe and SO 4 by pyrite oxidation; (2) Equilibrium exchange that moves some aqueous Fe to solid exchange sites; (3) Equilibration with the atmosphere, goethite, and calcite, which oxidizes aqueous Fe and precipitates goethite; (4) Mixing GW with highway deicing effluent results in increased Na cation exchange; and (5) Equilibration of the GW-D1 mixture with exchange sites and the same phases in (3) results in release of exchangeable Fe and additional goethite precipitation.
Model scenarios used in this study are described in Table 5, and more detailed results, including the (1) primary model with incremental pyrite reactions, (2) the primary model without calcite, and (3) the primary model  with biotite instead of pyrite, are provided in Supplementary Materials Table S5, S5A, and S5B. Stormwater runoff sample, D1, represents a water sample taken from I-95 in Connecticut after halite and CaCl 2 deicers were applied ( Table 5). The primary model includes calcite as an equilibrium phase, because calcite is common in the M c lithogeochemical group and acts to buffer the pH; 208 mM of calcite were allowed to dissolve. A default pe of 4.0 was used to define redox for the initial solution; although the DO measurements were available from field data, other redox couples required to define a representative pe were not Table 5 Results of forward geochemical modeling using PHREEQC All simulation steps react groundwater (GW) with pyrite in five steps (in increments of 0, 1, 5, 10, and 50 mM); equilibrium phases (CO 2(g), O 2(g) , goethite, and calcite); and 0.25 M of exchange sites. GW is mixed with deicer (D1) in 70:30 concentration ratio. See Table 4  measured. This is reasonable, however, since the defined pe is only used in the initial speciation of a solution and all subsequent "reaction" calculations for a given simulation react to redox equilibrium. Therefore, the pe of 4 (Table 5) was not significant compared to that in the primary reactions of the simulation. For each model, the exchange sites, and then O 2(gas) , CO 2(gas) , goethite, and calcite (if included) were set to reach equilibrium with the water, and most of the dissolved iron was precipitated as goethite. GW was then incrementally reacted with pyrite, which resulted in the oxidative dissolution of pyrite (or biotite as shown in SI Table S5B), the release of Fe 2+ , H + (acidity), and sulfate (Eq. 1), and, ultimately, the oxidation of Fe 2+ to Fe 3+ and precipitation of goethite. Siderite (FeCO 3 ) is supersaturated in the GW solution, but is undersaturated under the atmospheric conditions simulated for the rock-fill sites. Note that these models assume that all reactions progress to equilibrium.
For the primary model, the GW solution reacts with pyrite (steps 1-5), pH decreases from 7.75 to 5.9, and Fe 2+ , SO 4 2− , and H + are released (Tables 5 and S5; Eq. 1). Next, the solution reacts with the exchange sites (step 6), and equilibrates with atmospheric conditions, calcite, and goethite (step 7), which results in buffered pH and the precipitation of dissolved iron as goethite. Aqueous Fe 2+ adsorbs as FeX 2 (69 mM) in step 6 and oxidizes and precipitates as goethite (11 mM) in step 7. FeX 2 is the most abundant exchangeable cation, followed by other divalent cations, CaX 2 and MgX 2 . (Table S5). Some calcite is allowed to dissolve (190 mM), which helps buffer the pH to 7.69. The typical concentrations of sodium and chloride in M c bedrock groundwater are only at 0.48 to 0.28 mM (Table 4) and cation exchange is a negligible contribution to the exchange of Fe 2+ or H + . The next model step (8) includes the mixing of GW with deicer effluent at a ratio of 70:30 and equilibrium with exchange sites ("Exchange_2"), which causes additional release of exchangeable Fe. With deicer added, the dissolved Na + exchanges off other cations, including Fe 2+ and H + , and the concentrations of FeX 2 and dissolved Fe 3+ decreases.
Finally, the GW-D1 mixture is set to equilibrium with O 2 , CO 2 , calcite, and goethite, aqueous Fe 2+ is oxidized and precipitated as 12 mM Fe, which represents the exchanged FeX 2 . The goethite precipitated in step 10 is even greater than that precipitated in step 6; although the Fe in both cases was ultimately derived from the pyrite. The exchangeable iron on solid surfaces, therefore, represents a sizeable iron source that can potentially be mobilized by interaction with saline water and subsequent cation exchange.
Two additional models were run to determine the effects of conditions: (1) without calcite (SM Table S5A), and (2) with reaction of biotite instead of pyrite (SM Table S5B). The model without calcite indicates that, without buffering of acidity generated by pyrite dissolution and the oxidation of Fe 2+ , pH drops to 1.8 in the final step and slightly less goethite precipitates (10 mM) compared to the model with calcite (12 mM), as some of the iron remains as dissolved Fe 3+ due to the acidic conditions fewer calcium ions (from calcite dissolution) are available to mobilize exchangeable iron (SM Table S5A). The model run with biotite as the reacting iron source indicates high pH (up to 11.4), since 6 mol of H + are consumed for every mole of biotite; as a result, very low concentrations of dissolved and exchangeable iron are formed (SM Table S5B). This model uses an impure biotite with a log K from [86] (SM Table S5B).

Discussion
The presence or absence of iron fouling at rock-fill locations associated with NHDOT rock cuts determined through imagery analysis and field visits was the basis for this study. There was good agreement (86%) among iron fouling occurrence documented in the imagery analysis and from field visits. The 14% of sites that had differences between the imagery and the field visit (primarily iron fouling observed in the field that was not observed in the imagery) might be attributed to temporal iron occurrence factors or image quality (leaf off, ground cover) not examined as a part of this study. The possible change in iron fouling over time at a site could be investigated in future work.
A spatially continuous iron fouling probability map was developed for NH using a boosted regression tree model that identified geochemical, topographic, hydrologic, and geologic variables that best predicted the occurrence of iron fouling at rock-fill sites throughout the state. The presence of aluminum oxide (Al 2 O 3 ) and silicon dioxide (SiO 2 or quartz) in stream sediment as predictor variables in the model are notable. Quartz is typically abundant in siliceous, non-carbonate rocks and tends to monopolize the elemental composition in non-calcareous soils or sediments, a process known as "quartz dilution" [88,89]. The presence of quartz (SiO 2 ) as well as aluminum oxide (Al 2 O 3 ) as important predictor variables in the model suggests that a paucity of carbonate rocks could be associated with iron fouling.
The probability map developed from the model provides a preliminary desktop tool that can be used at any location in NH to evaluate the probability that iron fouling will be a concern at rock-fill sites where rock is sourced from the bedrock at that location. The model was developed using sites that are along roadways and may be less predictive in areas away from roadways and in areas with limited rock cut and rock-fill locations such as the northern portion of the state. At sites that have a high probability of iron fouling based on the map, additional parameters can be evaluated through field visits. Iron fouling is more likely to occur where associated rock cuts have red and brown weathering or red and brown staining below fracture lines. These are easy parameters to observe and require no additional equipment. While only 3 sites had iron fouling without observable staining or weathering at the associated rock cut, 246 sites had weathering or staining without having evidence of iron fouling, suggesting that more rock cut sites have iron present than end up containing iron fouling at fill sites. This provides additional evidence that the presence of iron in the rock-fill may not be the limiting factor in the occurrence of iron fouling at many sites. Local hydrology and (or) geochemical conditions at a site may be limiting the occurrence of iron fouling in these cases.
Iron fouling was more common at sites with low pH that result from iron oxidation and (or) pyrite dissolution, as supported by geochemical modeling results. Iron fouling also was more common at sites with high SC, often resulting from highway deicing runoff. Geochemical models using PHREEQC show that the increased iron fouling of waters in areas with low pH and high SC along sulfidic rock cuts are likely caused by pyrite dissolution and exacerbated by added salts that result in greater cation exchange. Geochemical models with biotite as the iron source (not pyrite) result in alkaline conditions with low dissolved iron or goethite precipitated. The typical concentrations of sodium and chloride in M c bedrock groundwater are low and cation exchange is a negligible contribution to the exchange of Fe 2+ or H + . Increased salt concentrations from deicing runoff and subsequent cation exchange results in greater amounts of displaced Fe 2+ that become oxidized and precipitated as goethite. This understanding could be helpful in planning road management alternatives with regard to deicers in locations where there is a high likelihood of iron fouling.
Several methods are identified for understanding preliminary factors related to iron fouling at rock-fill sites along roadways. The methods employed here may be used at other locations that experience iron fouling along roadways. The iron fouling probability map is a broad-scale tool for understanding the potential for iron fouling at sites in NH. Additionally, identifying the processes involved in iron fouling and understanding the role that local hydrology, geochemical conditions, and deicers play in the process is the first step in the minimization of iron fouling. This first look at fouling statewide suggests several factors and processes that could be further explored with additional field-scale investigation. Field or bench testing of geochemical controls of deicers and pH are two important potential drivers of iron mobility and fouling. Additionally, testing alternate deicers in iron-fouling prone areas may also be warranted. An area of future investigation is to explore the timing of application and subsequent fouling. Management decisions regarding road cuts, rock-fill placement, site design and the use of deicers will benefit from this understanding.

Conclusion
Iron fouling in runoff from rock cuts and adjacent to blasted rock-fill used for roadway construction is caused by the oxidation of iron sulfide minerals such as pyrite, and the subsequent release of acidity, sulfate, and dissolved metals including iron. Factors affecting iron fouling in wet areas adjacent to roadways were investigated by collecting field rock cut and aqueous physicochemical data; developing geochemical models; and developing exploratory predictive models. Road salt was examined as a primary concern related to iron fouling at these road construction sites.
Field water quality measurements indicated that the occurrence of iron fouling was associated with higher values of SC, lower concentrations of DO and lower pH compared to areas without iron fouling. Geochemical models illustrated the mechanisms for how iron fouling of waters increase along roadways built with fill from sulfidic rock cuts due to acid generation from pyrite dissolution and ferrous iron (Fe 2+ ) oxidation; and increase in areas with greater SC due to deicing runoff caused by cation exchange. Models demonstrate how higher salinity water can potentially increase iron mobility and aggravate ARD at these sites.
The study showed that the identification of iron fouling from aerial imagery is a useful screening technique for identifying potential iron fouling sites that can be used in similar studies in the future. The iron fouling probability map developed here will give NHDOT a tool to help predict and mitigate the occurrence of iron fouling at untested locations throughout NH. In addition, geochemical models provide a mechanism to better understand the relations between high SC and iron fouling observed in the field.
At sites that have a high probability of iron fouling based on the map, additional parameters can be evaluated through field visits. Iron fouling is more likely to occur where associated rock cuts have red and brown weathering or red and brown staining below fracture lines. There is evidence that the presence of iron in the rock-fill may not be the limiting factor in the occurrence of iron fouling at many sites. Local hydrology and (or) geochemical conditions at a site may be limiting the occurrence of iron fouling in these cases.