Arsenic in Petroleum-Contaminated Groundwater near Bemidji, Minnesota Is Predicted to Persist for Centuries

: We used a reactive transport model to investigate the cycling of geogenic arsenic (As) in a petroleum-contaminated aquifer. We simulated As mobilization and sequestration using surface complexation reactions with Fe(OH) 3 during petroleum biodegradation coupled with Fe-reduction. Model results predict that dissolved As in the plume will exceed the U.S. and EU 10 µ g/L drinking water standard for ~400 years. Non-volatile dissolved organic carbon (NVDOC) in the model promotes As mobilization by exerting oxygen demand, which maintains anoxic conditions in the aquifer. After NVDOC degrades, As re-associates with Fe(OH) 3 as oxygenated conditions are reestablished. Over the 400-year simulation, As transport resembles a “roll front” in which: (1) arsenic sorbed to Fe(OH) 3 is released during Fe-reduction coupled to petroleum biodegradation; (2) dissolved As resorbs to Fe(OH) 3 at the plume’s leading edge; and (3) over time, the plume expands, and resorbed As is re-released into groundwater. This “roll front” behavior underscores the transience of sorption as an As attenuation mechanism. Over the plume’s lifespan, simulations suggest that As will contaminate more groundwater than benzene from the oil spill. At its maximum, the model simulates that ~5.7 × more groundwater will be contaminated by As than benzene, suggesting that As could pose a greater long-term water quality threat than benzene in this petroleum-contaminated aquifer.


Introduction
Arsenic (As) is a common geogenic groundwater contaminant that originates in the sediments and rocks that comprise aquifer matrices. It is a known toxin and carcinogen associated with development of skin, bladder, and lung cancer [1]. Arsenic is regulated in drinking water at 10 µg/L in the United States [2] and in the European Union.
Although geogenic As can be released from sediments and rocks to groundwater through several mechanisms, the dominant process attributed to widespread groundwater contamination by As is a reductive dissolution of As-bearing ferric (Fe(III))-oxides coupled with microbial oxidation (i.e., biodegradation) of organic carbon (C) [3][4][5]. Recent work has highlighted four factors that are critical for accurate long-term assessments of groundwater contaminated by As via this mechanism [6]: (1) the source and biodegradability of organic C that drive Fe-reduction and subsequent As mobilization [7][8][9][10][11]; (2) the mineralogy, reactivity, and As content of Fe-oxide minerals [12][13][14][15]; (3) the retention of mobilized As by incorporation in secondary mineral phases and/or resorption of As to sediment mineral surfaces [7,16,17]; and (4) the groundwater flow rate [4,[17][18][19]. the study site [46]. After initial clean-up efforts, approximately 25% of the spilled oil infiltrated into the subsurface and settled on the water table. This oil body remains a source of dissolved organic C in groundwater. Biodegradation of the dissolved organic C is coupled with anaerobic TEAPs including Fe-reduction and methanogenesis [47,48]. Since 1984, the site has been studied extensively to monitor the long-term evolution of the petroleum hydrocarbon plume, the natural attenuation capacity of the aquifer, and secondary water quality impacts such as depleted dissolved oxygen (DO), elevated Fe 2+ , and pH changes triggered by biodegradation of organic C (for a review of research at the Bemidji oil spill site, see [36]). is the approximate center of the oil body and is the zero-reference point for this study. Values for background aqueous phases used for initial conditions and constant concentration fluxes are determined from the upgradient well labeled "Background". The 10 µg/L As contour is from Cozzarelli et al. [29].
Recent work from our group documented the mobilization of geogenic As (primarily As(III)) caused by biodegradation coupled with Fe-reduction at the oil spill site [29]. Results from groundwater sampling showed a secondary plume of dissolved As that overlapped with plumes of benzene, dissolved Fe 2+ , and depleted DO. Sediment and groundwater surveys showed Fe and As depletion in sediment near the oil, corresponding to elevated dissolved Fe 2+ and As in groundwater. Dissolved Fe 2+ and As decreased along the groundwater flowpath to the leading edge of the petroleum hydrocarbon plume (i.e., transition zone), where As accumulated in sediment by association with Fe(III) in a Ferich "curtain" [30,31]. Comparisons of data from 1993 [38], 2008 [49], and 2013-2015 [30] have suggested that Fe(III) in the curtain is reduced as the plume expands, releasing Fe 2+ and As from the curtain into groundwater [30].

Methods
The model presented in this study adds to previous modeling studies of the Bemidji oil spill site that have investigated secondary processes affecting water quality [28,41,[50][51][52][53]. Here, we use a surface complexation model with As sorption and desorption from Fe(OH)3 during biodegradation to simulate As cycling in a two-dimensional cross-section along the plume transect ( Figure 1). Model simulations were performed using PHT3D is the approximate center of the oil body and is the zero-reference point for this study. Values for background aqueous phases used for initial conditions and constant concentration fluxes are determined from the upgradient well labeled "Background". The 10 µg/L As contour is from Cozzarelli et al. [29].
In 1979, a pipeline rupture released 10,700 barrels of crude oil to the land surface of the study site [46]. After initial clean-up efforts, approximately 25% of the spilled oil infiltrated into the subsurface and settled on the water table. This oil body remains a source of dissolved organic C in groundwater. Biodegradation of the dissolved organic C is coupled with anaerobic TEAPs including Fe-reduction and methanogenesis [47,48]. Since 1984, the site has been studied extensively to monitor the long-term evolution of the petroleum hydrocarbon plume, the natural attenuation capacity of the aquifer, and secondary water quality impacts such as depleted dissolved oxygen (DO), elevated Fe 2+ , and pH changes triggered by biodegradation of organic C (for a review of research at the Bemidji oil spill site, see [36]).
Recent work from our group documented the mobilization of geogenic As (primarily As(III)) caused by biodegradation coupled with Fe-reduction at the oil spill site [29]. Results from groundwater sampling showed a secondary plume of dissolved As that overlapped with plumes of benzene, dissolved Fe 2+ , and depleted DO. Sediment and groundwater surveys showed Fe and As depletion in sediment near the oil, corresponding to elevated dissolved Fe 2+ and As in groundwater. Dissolved Fe 2+ and As decreased along the groundwater flowpath to the leading edge of the petroleum hydrocarbon plume (i.e., transition zone), where As accumulated in sediment by association with Fe(III) in a Fe-rich "curtain" [30,31]. Comparisons of data from 1993 [38], 2008 [49], and 2013-2015 [30] have suggested that Fe(III) in the curtain is reduced as the plume expands, releasing Fe 2+ and As from the curtain into groundwater [30].

Methods
The model presented in this study adds to previous modeling studies of the Bemidji oil spill site that have investigated secondary processes affecting water quality [28,41,[50][51][52][53]. Here, we use a surface complexation model with As sorption and desorption from Fe(OH) 3 during biodegradation to simulate As cycling in a two-dimensional cross-section along the plume transect ( Figure 1). Model simulations were performed using PHT3D (v2.10) [54],

Model Domain and Hydrogeologic Parameters
The model domain and hydrogeologic parameters are identical to those in the original model [28]. The model simulation covers a 400-year period starting from the occurrence of the spill in August 1979. The model analysis domain is a cross-section of the groundwater plume that extends horizontally from 45 m upgradient from the center of the oil body to 215 m downgradient. The vertical dimension of the domain covers an elevation range of 417 to 424 m above mean sea level, with the upper boundary being the water table. In order to avoid boundary effects during model simulations, the computational grid was further extended to 475 m downgradient and lowered to an elevation of 410 m. Individual cells in the computational grid have dimensions of 4.3 m (horizontal) by 0.47 m (vertical).
Hydrogeologic parameters were informed using grain-size analyses from 269 sediment samples from the aquifer. The arithmetic mean for porosity from these samples was 0.38 (± 0.04) [42]. Given the narrow range in variability, a homogeneous porosity of 0.38 was applied uniformly to the model grid. Hydraulic conductivity (K) in the aquifer is spatially heterogeneous, ranging from 10 −3 to 10 m/d, with a median of 6.3 m/d based on 58 slug tests [43]. K was applied heterogeneously to the model grid based on geostatistical interpolation of the grain size analyses previously performed on Bemidji sediments [42], with geometric and arithmetic means for K of 5.1 m/d and 6.3 m/d, respectively. Longitudinal and transverse dispersivities were applied uniformly at values of 1 m and 0.04 m, respectively, based on a field tracer test [41]. A steady-state flow field was imposed with a constant head boundary of 424 m above sea level at the left boundary to produce the observed average hydraulic gradient of 0.0035. A constant recharge rate of 4.88 × 10 −4 m/d was applied uniformly to the upper water table boundary.

Geochemical Formulation
The geochemical model uses a partial equilibrium approach, which assumes that organic C oxidation is the rate-limiting step, and the associated reduction reaction proceeds as an equilibrium reaction [58]. Using this approach, an organic molecule is irreversibly broken down into its constitutive elements and oxidation states (in parentheses) at a specified rate, r. An example for benzene is shown in Equation (1): In this example, one mole of benzene is transformed into six moles of C(−1) and H(+1), which are then reintroduced into the equilibrium solution. C(−1) is converted into HCO 3 − by equilibrium reactions using the most thermodynamically favorable TEAP. All subsequent geochemical reactions then proceed by thermodynamic equilibrium.
The original reactive transport model [28] described in detail the model formulation for biodegradable organic molecules. The chemical composition of the non-aqueous phase liquid (NAPL) is reported in Table S1. Benzene, ethylbenzene, and xylenes were modeled as one component (BEX); toluene, non-volatile dissolved organic C (NVDOC), short-chain n-alkanes, and long-chain n-alkanes were the other modeled organic constituents. Each constituent had its unique dissolution rate from the NAPL into groundwater and firstorder biodegradation rate under different terminal electron accepting processes. These rate parameters were calibrated to match field data from 1987, 1993, and 2008. For a detailed description of the organic biodegradation formulation in the model, we refer the reader to [28].
Surface complexation was modeled using the generalized two-layer model and database described in [59] for As complexed to hydrous ferric oxide (HFO). The twolayer model defines "Type 1" and "Type 2" sites, alternatively termed strong and weak sites, respectively. Type 1 sites, defined in the PHREEQC database as Hfo_s ("_s" denotes "strong"), correspond to a smaller subset of surface sites with a high affinity for cation sorption. Type 2 sites, defined as Hfo_w in the PHREEQC database ("_w" denotes "weak"), are the total reactive sites available for sorption of protons, cations, and anions as determined by observed sorption maxima [59]. To simplify the model, Type 1 sites were excluded because they do not sorb As oxyanions. In our model, HFO was user-defined as Fe(OH) 3 . The specific surface area (600 m 2 /g) and surface site density for Hfo_w (0.2 mol/mol Fe) were from the same database [59].
Complexation of As(V) and As(III) to Fe(OH) 3 is described thermodynamically by empirically derived intrinsic surface complexation constants and acid dissociation constants ( Table 1). All complexes are represented as monodentate complexes. Using Reaction 7 in Table 1 as an example, Hfo_wOH represents a single Type 2 sorption site on HFO, and the formation of Hfo_wH 2 AsO 4 represents the formation of a monodentate complex of a single arsenate molecule to one Type 2 site on HFO. The best-estimate logK values for surface complexation reactions, as determined by previous studies (see references in Table 1), were used in the model. Upper and lower 95% confidence intervals from those previous studies were used in sensitivity tests (see Section 2.5). Carbonate complexation is also considered in the surface complexation model (Table 1) because carbonate has been shown to outcompete As for sorption sites and to displace sorbed As [62]. Sulfate and phosphate, which are other potential sorption competitors, are not included because sulfate is only present sporadically at low concentrations (generally 0.5-5 mg/L) and phosphate is below detection (<0.02 mg/L) in the Bemidji aquifer.
Spatiotemporal model output was constrained by multiple dissolved and sediment chemical constituents. Arsenic concentrations were constrained by measured groundwater and sediment concentrations collected along the plume transect from 2011-2015 [29,30,63]. Other parameters (pH, DO, Fe, NVDOC, BEX) were constrained by groundwater and sediment chemistry data collected in 1987, 1993, and 2008 (data reported in [64]). Additional concentration data for pH, DO, and Fe from 2011-2015 were used to constrain the model (data reported in [29,30,63]).

Initial Conditions
Initial parameters for organic C constituents, including the NAPL density, initial total C mass, and mass fractions of BEX, toluene, NVDOC, short-chain n-alkanes, and long-chain n-alkanes, are presented in Table S1. For additional details and references for these parameters, we refer the reader to [28].
Initial concentrations for dissolved constituents (total carbonates, CH 4 , Ca 2+ , Cl − , Fe 2+ , Mg 2+ , Na + , DO, Mn 2+ , N 2 , pH) were determined using measurements from an uncontaminated background well upgradient from the oil spill site (see Figure 1 for the background well location). Concentrations from the background well from six sampling periods between 1986-1995 were averaged and then equilibrated with calcite and charge-balanced using PHREEQC-2 before application in the model. The average background concentrations for these constituents have remained relatively constant; background concentrations measured in the 2010s were similar to those measured in the 1980s and 1990s. Initial model input concentrations and concentrations after equilibration are found in Table 2. Constant concentration flux boundaries with the initial solution chemistry (excluding As) were imposed at the left boundary and water table boundary of the model domain. Table 2. Initial conditions and background concentrations of inorganic aqueous species. The second column shows average concentrations from an uncontaminated background well, 310e (see Figure 1), located upgradient; the third column shows PHREEQC-2 batch simulation results after equilibrating with calcite, Fe(OH) 3 (a), and MnO 2 and establishing charge balance. The equilibrated concentrations were used for constant concentration fluxes at the left and top model boundary. Arsenic was removed from the constant concentration flux boundaries because it was non-detectable (<1 µg/L; <1.35 × 10 −8 M).
Extracted background Fe(III) ranged from 0.0228 to 0.0657 mol/L v (775-2224 mg/kg), with a mean concentration of 0.0412 mol/L v (1394 mg/kg), where L v is the bulk aquifer volume (sediment plus pore space). To determine how much Fe(III) was reducible, we compared HCl-extractable Fe(III) from uncontaminated background sediment with Fe(III) in highly contaminated sediment near the oil body [30]. Low, but measurable, amounts of HCl-extractable Fe(III) remain near the oil body, where the predominant TEAP is methanogenesis, suggesting that some residual HCl-extractable Fe(III) may not be easily reducible by microbes. Thus, for the model, the difference between the average background concentration and the residual Fe(III) near the oil (0.0370 mol/L v , or 1254 mg/kg) was applied homogeneously to the model domain as the initial reducible Fe(OH) 3 . This is an adjustment from the original model, which used an initial concentration of 0.0288 mol/L v based on higher measured residual Fe(III) near the oil body observed in 1993 [38] compared to residual Fe(III) measured near the oil in 2015, the basis for the current study.
The surface complexation model used to describe As cycling was equilibrated with the dissolved As in the initial solution (Table 2) to approximate observed concentrations of As sorbed to Fe(OH) 3 . The initial complexed As concentration was calibrated to the range of sorbed As concentrations determined from a phosphate extraction of background sediment [30], which quantifies operationally defined sorbed As. Based on preliminary model simulations, the initial sorbed As concentration was set to 1.62 × 10 −2 mmol/L v in the model, which was 60% of the mean phosphate-extractable concentration from background sediment. We determined this was a reasonable initial concentration because the phosphate extraction is not specific to As sorbed to Fe(OH) 3 ; it will also extract As sorbed to more recalcitrant minerals that would likely not be reduced during biodegradation. Thus, the total measured average phosphate-extractable As would likely be an overestimate of As mobilizable due to biodegradation reactions. After equilibration, two As(V) forms were the predominant complexes: Hfo_wAsO 4 2− and Hfo_wOHAsO 4 3− , due to the circumneutral pH, with concentrations of 6.43 × 10 −3 mmol/L v and 9.76 × 10 −3 mmol/L v , respectively. These two forms accounted for 99.8% of the initial sorbed As. The forms Hfo_wH 2 AsO 4 , Hfo_wHAsO 4 − and Hfo_wH 2 AsO 3 were negligible.

Calibration Approach
Initial concentrations of sorbed As and Fe(OH) 3 were calibrated so that their concentrations (1) were reasonably similar to observed measurements in background sediments (see description in Section 2.3), and (2) allowed the model to match three key characteristics of the As plume to field measurements: the maximum concentration of dissolved As, the total dissolved As mass, and the location of the 10 µg/L dissolved As front. The first two calibration characteristics were chosen to reflect the overall magnitude of As mobilization. The 10 µg/L front was chosen because it is the regulated maximum contaminant level (MCL) for As in the United States and the European Union. Because of its relevance for groundwater regulation and safety of drinking water, our highest priority for calibration was to accurately simulate the location of the leading edge of the 10 µg/L front.
While we also attempt to accurately simulate spatially distributed As concentrations measured in the As plume, we chose to prioritize the three coarser-scale metrics listed above because As concentrations can fluctuate in individual wells for a variety of reasons that are not related to aquifer processes, including sampling practices [66] and well construction [67][68][69]. Thus, comparing point-to-point As concentrations is less effective in providing insight into the biogeochemical reactions that influence As than examining these larger-scale plume characteristics.

Model Sensitivity
We examined the sensitivity of the model output to the three parameters that were modified from, or added to, the original model. They were the initial concentration of Fe(OH) 3 in aquifer sediments, the initial concentration of sorbed As, and intrinsic surface complexation constants for As sorption onto Fe(OH) 3 . For the sensitivity tests, the initial concentration of Fe(OH) 3 was adjusted to the minimum and maximum measured values in background sediment. The concentration of sorbed As used in the model was halved and doubled, which was within the range of measured background sorbed As [30]. The intrinsic surface complexation constants for As sorption were adjusted to the upper and lower 95% confidence intervals reported by Dzombak and Morel [59]. Parameter values described in previous sections were used as the "base case" to assess model sensitivity.
We evaluated the model's sensitivity using the sensitivity coefficient approach [70], which calculates the sensitivity of a dependent variable to a model parameter as the partial derivative of the dependent variable with respect to the model parameter. For easier comparison of different model parameters, the sensitivity coefficient can be normalized to a dimensionless form (Equation (2)): where X i,k is the sensitivity coefficient of the dependent variable y with respect to the kth independent variable (or parameter) at the ith time point, which in this case is 37 years. y i,1 is the base case dependent variable and y i,2 is the altered dependent variable; α k,1 is the base case independent variable; α k,2 is the manipulated independent variable.

Aqueous Species
Simulated aqueous species that are important for As cycling at the Bemidji site (pH, DO, Fe 2+ ) are shown for 2016, 37 years after the oil spill ( Figure 2). pH is important because the predominant dissolved and sorbed As species are pH-dependent, as is the charge on Fe(OH) 3 surface sites. Over time, pH decreases in the model, which is attributed to the production of H + from methanogenic biodegradation and cation exchange of dissolved cations with surficial H + [28] (Figure 2A). For an in-depth discussion of reactions affecting pH in the Bemidji aquifer, see [44]. DO is depleted from groundwater near the oil due to the onset of aerobic biodegradation shortly after the spill occurred ( Figure 2B). Fe 2+ is a byproduct of biodegradation coupled with reductive dissolution of Fe(OH) 3 . The Fe 2+ plume in the 2016 simulation ( Figure 2C) has extended 30 m farther downgradient than the previously simulated Fe 2+ plume for 2007 [28] and is consistent with the observed Fe 2+ plume shape determined from groundwater sampling in 2013 [29].
to targeted field sampling aimed at capturing higher concentrations of As in groundwater (i.e., the sampling was not random), causing a shift in the cumulative distribution of field concentrations to higher values. We also note that the simulated maximum concentration of As is nearly twice as high as the observed maximum (428 vs. 230 µg/L). However, it is likely that field sampling did not capture the true dissolved As maximum in groundwater. The model results show that this high As concentration is very localized; only a single grid cell showed a concentration greater than 175 µg/L ( Figure 3A).
Our second approach for comparing modeled vs. measured As concentrations was point-specific. Measured field concentrations along the plume transect were compared to simulated As concentration profiles at two elevations in the plume: 422 m and 419 m above mean sea level. These elevations approximate the vertical center and the bottom of the plume, respectively. The modeled concentration profiles show a reasonably close fit to measured concentrations ( Figure 3B) [29,31], with the exception of the simulated As maximum at 422 m elevation, which was not observed during field sampling. Concentrations from other field samples were within ~5 µg/L of the simulated concentration.   2+ , and (D) total dissolved As (As T ) for 37 years after the oil spill (2016). The solid white line for dissolved As T represents the simulated 10 µg/L front. The dashed white line for As T represents the approximate extent of the 10 µg/L front based on field sampling (data available in [29]). The gray box is the reference for the oil body extent. Top boundary is the water table.
Simulated concentrations of dissolved As are shown as total dissolved As (As(III) + As(V)) ( Figure 2D). In the model, nearly all dissolved As is As(III); As(V) does not exceed 1 µg/L anywhere in the model domain in 2016, in general agreement with field measurements from 2013 that showed that most As in groundwater (>80%) is As(III) [29]. The simulated 10 µg/L front, which was prioritized during calibration, matches the location of the front observed from field sampling within two meters. The total dissolved As mass from the model (0.14 moles) is within the range estimated from groundwater sampling (0.11 to 0.15 moles) [31]. The center of As mass occurs slightly further downgradient (65 m) in the model simulation ( Figure 2) compared to the estimated center of mass determined from field sampling (~40 m) [29].
We compared the calibrated model output for As to field samples in two ways. The first was a broad approach, in which we generated cumulative distribution functions comparing the model output concentrations in each cell in the As plume and actual As concentrations measured in the field [63]. This comparison shows that the model reasonably simulates measured As concentrations in groundwater, with nearly all concentrations less than 175 µg/L ( Figure 3A). One difference is that the model output has a larger proportion of lower (but non-zero) concentrations compared to field measurements. We attribute this to targeted field sampling aimed at capturing higher concentrations of As in groundwater (i.e., the sampling was not random), causing a shift in the cumulative distribution of field concentrations to higher values. We also note that the simulated maximum concentration of As is nearly twice as high as the observed maximum (428 vs. 230 µg/L). However, it is likely that field sampling did not capture the true dissolved As maximum in groundwater. The model results show that this high As concentration is very localized; only a single grid cell showed a concentration greater than 175 µg/L ( Figure 3A).  Open circles designate samples that were reported as <10 µg/L, and thus are plotted at half the detection limit. The dashed horizontal line shows the 10 µg/L MCL for As.

Solid Phase Species
Model results show that the simulated depletion of Fe(OH)3 in the aquifer occurs in two distinct zones, which we term a "dual zone effect" (Figure 4). Close to the oil, Fe(OH)3 is almost entirely depleted from aquifer sediments. Surrounding this area is another zone where Fe(OH)3 has been only partially depleted. This dual zone effect corresponds almost exactly to the pattern of elevated Fe 2+ in groundwater ( Figure 2C). These patterns are influenced by groundwater pH, which is, in turn, controlled by the predominant TEAP. Near the oil, methanogenic conditions produce H + that lower pH (Figure 2A), which thermodynamically promotes Fe-reduction. Immediately downgradient from this zone of Fe(OH)3 depletion, only partial depletion of Fe(OH)3 occurs because the process of Fereduction raises pH, which thermodynamically disfavors Fe-reduction. Our second approach for comparing modeled vs. measured As concentrations was point-specific. Measured field concentrations along the plume transect were compared to simulated As concentration profiles at two elevations in the plume: 422 m and 419 m above mean sea level. These elevations approximate the vertical center and the bottom of the plume, respectively. The modeled concentration profiles show a reasonably close fit to measured concentrations ( Figure 3B) [29,31], with the exception of the simulated As maximum at 422 m elevation, which was not observed during field sampling. Concentrations from other field samples were within~5 µg/L of the simulated concentration.

Solid Phase Species
Model results show that the simulated depletion of Fe(OH) 3 in the aquifer occurs in two distinct zones, which we term a "dual zone effect" (Figure 4). Close to the oil, Fe(OH) 3 is almost entirely depleted from aquifer sediments. Surrounding this area is another zone where Fe(OH) 3 has been only partially depleted. This dual zone effect corresponds almost exactly to the pattern of elevated Fe 2+ in groundwater ( Figure 2C). These patterns are influenced by groundwater pH, which is, in turn, controlled by the predominant TEAP. Near the oil, methanogenic conditions produce H + that lower pH (Figure 2A), which thermodynamically promotes Fe-reduction. Immediately downgradient from this zone of Fe(OH) 3 depletion, only partial depletion of Fe(OH) 3 occurs because the process of Fe-reduction raises pH, which thermodynamically disfavors Fe-reduction. two distinct zones, which we term a "dual zone effect" (Figure 4). Close to the oil, Fe(OH)3 is almost entirely depleted from aquifer sediments. Surrounding this area is another zone where Fe(OH)3 has been only partially depleted. This dual zone effect corresponds almost exactly to the pattern of elevated Fe 2+ in groundwater ( Figure 2C). These patterns are influenced by groundwater pH, which is, in turn, controlled by the predominant TEAP. Near the oil, methanogenic conditions produce H + that lower pH (Figure 2A), which thermodynamically promotes Fe-reduction. Immediately downgradient from this zone of Fe(OH)3 depletion, only partial depletion of Fe(OH)3 occurs because the process of Fereduction raises pH, which thermodynamically disfavors Fe-reduction. The effect of the "dual-zone" of Fe-reduction has confounding effects for As cycling. In the zone closest to the oil, As(V) initially complexed to Fe(OH)3 (as Hfo_wAsO4 2− , and Hfo_wOHAsO4 3− ) is entirely stripped from sediment due to extensive Fe-reduction (Figure 5), creating a plume of dissolved As ( Figure 2D). This modeled depletion of As has been observed by chemical extractions of sediments collected from this zone [30]. However, in the zone of partial Fe-depletion, sediments actually appear to be enriched in As via sorption of As(III) as Hfo_wH2AsO3 ( Figure 5C), despite the fact that Fe(III) is actively The effect of the "dual-zone" of Fe-reduction has confounding effects for As cycling. In the zone closest to the oil, As(V) initially complexed to Fe(OH) 3 (as Hfo_wAsO 4 2− , and Hfo_wOHAsO 4 3− ) is entirely stripped from sediment due to extensive Fe-reduction ( Figure 5), creating a plume of dissolved As ( Figure 2D). This modeled depletion of As has been observed by chemical extractions of sediments collected from this zone [30]. However, in the zone of partial Fe-depletion, sediments actually appear to be enriched in As via sorption of As(III) as Hfo_wH 2 AsO 3 ( Figure 5C), despite the fact that Fe(III) is actively being reductively dissolved. We attribute this As(III) enrichment to the fact that enough residual Fe(OH) 3 remains in this region of the "dual-zone" to sorb As(III) that was stripped from the upgradient part of the "dual-zone". Sediment samples from cores collected in this part of the plume (n = 55) similarly showed spatially heterogeneous sediment-bound As concentrations, ranging from 8.9 × 10 −3 to 0.13 mmol/L v , with a mean of 0.045 mmol/L v (range = 0.4 to 5.8 mg/kg; mean = 2.05 mg/kg [30]). A subset of these samples (n = 25) showed that~75% of As was sorbed (~0.033 mmol/L v ; 1.5 mg/kg), slightly higher than the modeled values (up to 0.025 mmol/L v ), though in general agreement ( Figure 5).
In the model, some dissolved As breaks through the zone of elevated Hfo_wH 2 AsO 3 and is transported further downgradient, where it is oxidized to form As(V) complexes as Hfo_wAsO 4 2− ( Figure 5). Interestingly, the modeled DO in this region is~0 ( Figure 2B), indicating that the oxidation of As(III) perhaps contributes to DO depletion, or some other unidentified oxidant is responsible for As(III) oxidation. In this zone, the concentration of Hfo_wAsO 4 2− is two times greater than background concentration due to resorption ( Figure 5). Other As(V) complexes do not contribute to As retention.

Long-Term Cycling of As in the Petroleum Hydrocarbon Plume
After calibrating the model to field data collected 37 years after the spill, a simulation was run for 400 years to provide insight into the long-term cycling of As. The 400-year simulation results are presented as cross-sectional concentration maps for pH, DO, benzene, ethylbenzene, and xylenes (BEX), non-volatile dissolved organic carbon (NVDOC), Fe 2+ , Fe(OH) 3  and is transported further downgradient, where it is oxidized to form As(V) complexes as Hfo_wAsO4 2− ( Figure 5). Interestingly, the modeled DO in this region is ~0 (Figure 2B), indicating that the oxidation of As(III) perhaps contributes to DO depletion, or some other unidentified oxidant is responsible for As(III) oxidation. In this zone, the concentration of Hfo_wAsO4 2− is two times greater than background concentration due to resorption (Figure 5). Other As(V) complexes do not contribute to As retention.

Long-Term Cycling of As in the Petroleum Hydrocarbon Plume
After calibrating the model to field data collected 37 years after the spill, a simulation was run for 400 years to provide insight into the long-term cycling of As. The 400-year simulation results are presented as cross-sectional concentration maps for pH, DO, benzene, ethylbenzene, and xylenes (BEX), non-volatile dissolved organic carbon (NVDOC),  Long-term simulation results show that H + produced from methanogenic biodegradation lowers the groundwater pH for 100 years, at which point a pH minimum of~6.5 occurs near the oil body ( Figure 6). Between 100 and 150 years, pH begins to revert to its background condition after NVDOC biodegradation is complete. Without the added acidity from methanogenic NVDOC biodegradation, recharging groundwater from the constant flux boundary at the left model boundary and the water table recharge boundary allows pH to return to near-background conditions.
The plume of depleted DO also expands downgradient as the petroleum hydrocarbon plume evolves over time ( Figure 6). DO does not return to background conditions as quickly as pH. In fact, DO remains below background concentrations through the 400-year simulation (data not shown), and the depleted DO plume extends well past the 400 m model domain beyond 200 years. The slower return to background DO is likely due to DO consumption via multiple processes, including aerobic respiration of residual BEX and NVDOC near the oil body. A more important consideration for DO consumption is that the reducing conditions in the plume have created substantial chemical oxygen demand in the aquifer by building up a repository of electron donors in groundwater and sediments. For example, recharging DO from the top and left boundaries oxidizes aqueous Fe 2+ , sorbed Fe(II), and FeCO 3 to Fe(OH) 3 . This can be observed by a "halo" of a relative increase in Fe(OH) 3 (note the halo of light blue in the zone of dark blue at 150 years in Figure 6W). This halo region corresponds to a decrease in concentrations of aqueous Fe 2+ (compare concentrations from 0 to~75 m in Figure 6R,S), sorbed Fe(II) (compare concentrations from 0 to~75 m in Figure 6Z,AA), and FeCO 3 (compare concentrations from 0 to~75 m in Figure 6DD,EE). By 200 years, the halo of increased Fe(OH) 3 extends to greater depths and further downgradient in the model domain. By 200 years, sorbed Fe(II) and FeCO 3 substantially decrease in concentration due to oxidation. Based on these simulations, we suspect that background DO concentrations will not be reestablished until both biologic and chemical oxygen demand in the aquifer have been satisfied by oxidizing remaining organic and inorganic electron donors in the plume. Model simulations suggest this will take more than 400 years.  In the long-term simulations, As transport displays a "roll front" behavior, in which As is reductively dissolved, released into groundwater, and transported to the plume's leading edge, where it then sorbs to Fe(OH) 3 (Figure 7). Over time, the plume expands and As is remobilized as As-bearing Fe(OH) 3 reductively dissolves. Then dissolved As is re-attenuated at the plume's new leading edge further downgradient. This can be seen in the model output by observing changes in the location of the center of As mass, the 10 µg/L front, and the zone of accumulated sorbed As over time. For example, between 50 and 100 years, the center of As mass migrates from~80 m downgradient from the oil body to 200 m, the 10 µg/L front migrates from 128 m to 225 m downgradient, and the zone of accumulated sorbed As migrated from 100-200 m to 200-300 m downgradient (Figure 7). These simulations suggest that the center of As mass migrates faster (average rate of 2.4 m/y) than the 10 µg/L front (average rate of 1.9 m/y) over this timeframe. By 100 years, the dissolved As mass reaches its simulated maximum of 0.4 moles (Figure 8). From 100 to 150 years, the advancement of the center of mass and the 10 µg/L front slows, corresponding to a substantial loss of dissolved As mass; by 150 years, the total dissolved As mass decreases to 0.08 moles, 20% of the maximum (Figure 8), though concentrations in the center of mass still remain >100 µg/L in some cells. After 150 years, the dissolved As plume begins to stabilize, and the 10 µg/L front reaches 300 m downgradient, where it remains fixed until 387 years, at which time dissolved As in all cells finally reaches concentrations of <10 µg/L. During this period, the total mass of dissolved As also stabilizes at 0.02 moles, a 95% decrease from the simulated maximum.
Similar to the 37-year simulation, the long-term model results suggest that the primary mechanism responsible for As retention is sorption of As(III) to Fe(OH)3 as Hfo_wH2AsO3, which forms a ring around the plume fringes ( Figure 7M-P). Beyond the ring of sorbed As(III), some oxidized As(V) is also sorbed to Fe(OH)3 ( Figure 7E-H); the predominant sorbed As(V) species is Hfo_wAsO4 2− . Close examination of the time series evolution of As surface species provides hints that the As(III) complex is not thermodynamically stable with recharging background groundwater and will not persist indefinitely. Rather, As(III) complexes will eventually be oxidized to As(V) complexes, serving as another example of chemical oxygen demand. From 50 to 200 years, there is a loss of Hfo_wH2AsO3 and a corresponding increase in Hfo_wAsO4 2− and Hfo_wOHAsO4 3− complexes at the water table boundary, suggesting that DO in recharging groundwater transforms the reduced As(III) complexes to oxidized As(V) complexes. Eventually, we suspect By 100 years, the dissolved As mass reaches its simulated maximum of 0.4 moles (Figure 8). From 100 to 150 years, the advancement of the center of mass and the 10 µg/L front slows, corresponding to a substantial loss of dissolved As mass; by 150 years, the total dissolved As mass decreases to 0.08 moles, 20% of the maximum (Figure 8), though concentrations in the center of mass still remain >100 µg/L in some cells. After 150 years, the dissolved As plume begins to stabilize, and the 10 µg/L front reaches 300 m downgradient, where it remains fixed until 387 years, at which time dissolved As in all cells finally reaches concentrations of <10 µg/L. During this period, the total mass of dissolved As also stabilizes at 0.02 moles, a 95% decrease from the simulated maximum.
Similar to the 37-year simulation, the long-term model results suggest that the primary mechanism responsible for As retention is sorption of As(III) to Fe(OH) 3 as Hfo_wH 2 AsO 3 , which forms a ring around the plume fringes ( Figure 7M-P). Beyond the ring of sorbed As(III), some oxidized As(V) is also sorbed to Fe(OH) 3 (Figure 7E-H); the predominant sorbed As(V) species is Hfo_wAsO 4 2− . Close examination of the time series evolution of As surface species provides hints that the As(III) complex is not thermodynamically stable with recharging background groundwater and will not persist indefinitely. Rather, As(III) complexes will eventually be oxidized to As(V) complexes, serving as another example of chemical oxygen demand. From 50 to 200 years, there is a loss of Hfo_wH 2 AsO 3 and a corresponding increase in Hfo_wAsO 4 2− and Hfo_wOHAsO 4 3− complexes at the water table boundary, suggesting that DO in recharging groundwater transforms the reduced As(III) complexes to oxidized As(V) complexes. Eventually, we suspect that all As(III) will be transformed back to As(V) complexes, with the predominant complex determined by the groundwater pH. From these simulations, we conclude that the final step in As cycling in the aquifer will be the return of As to As(V) sorbed to Fe(OH) 3 , though this process is not yet complete by the end of the 400-year simulation.

21, 13, x FOR PEER REVIEW
15 of 24 final step in As cycling in the aquifer will be the return of As to As(V) sorbed to Fe(OH)3, though this process is not yet complete by the end of the 400-year simulation. Simulations suggest that the presence of NVDOC is the dominant cause of the persistence of As in groundwater for the first 100 years (Figure 8). While NVDOC is present in groundwater (0-100 years), As remains in groundwater at high concentrations, and the plume migrates downgradient at a steady rate. However, once NVDOC has been biodegraded completely (between 100 and 150 years; Figure 6N,O), the amount of dissolved As decreases ( Figures 7B,C and 8). From these observations, we propose that the biodegradation of NVDOC sustains the reducing conditions needed to promote Fe-reduction and As mobilization, resulting in continued release of Fe 2+ and As into groundwater. However, after NVDOC becomes (nearly) completely biodegraded after 100 years, the dissolved As mass rapidly declines, most likely because Fe(OH)3 is no longer reductively dissolved due to NVDOC biodegradation. Therefore, more As is being removed from groundwater by sorption at the plume fringes than is being introduced into groundwater via Fe-reduction, causing a net removal of As from groundwater.

Model Sensitivity
Of the three dependent plume characteristics considered in the sensitivity analysis, the location of the 10 µg/L front was relatively insensitive to any of the independent variables assessed (Table 3). Doubling and halving the initial sorbed As concentration moderately influenced the total dissolved mass and the maximum dissolved As, but they had little effect on the location of the 10 µg/L front. By far, the model was most sensitive to the initial Fe(OH)3 concentration. In particular, the model is sensitive to underestimations of initial Fe(OH)3. Because sorption to Fe(OH)3 is the main retention mechanism for As, a poorly quantified initial concentration for Fe(OH)3 can cause substantial errors in simulating As transport by either over-or underestimating sorption sites for As, which then either incorrectly enhances or lessens As retention, respectively. Collecting a robust background dataset will minimize errors that could result from natural heterogeneity and a Simulations suggest that the presence of NVDOC is the dominant cause of the persistence of As in groundwater for the first 100 years (Figure 8). While NVDOC is present in groundwater (0-100 years), As remains in groundwater at high concentrations, and the plume migrates downgradient at a steady rate. However, once NVDOC has been biodegraded completely (between 100 and 150 years; Figure 6N,O), the amount of dissolved As decreases ( Figure 7B,C and Figure 8). From these observations, we propose that the biodegradation of NVDOC sustains the reducing conditions needed to promote Fe-reduction and As mobilization, resulting in continued release of Fe 2+ and As into groundwater. However, after NVDOC becomes (nearly) completely biodegraded after 100 years, the dissolved As mass rapidly declines, most likely because Fe(OH) 3 is no longer reductively dissolved due to NVDOC biodegradation. Therefore, more As is being removed from groundwater by sorption at the plume fringes than is being introduced into groundwater via Fe-reduction, causing a net removal of As from groundwater.

Model Sensitivity
Of the three dependent plume characteristics considered in the sensitivity analysis, the location of the 10 µg/L front was relatively insensitive to any of the independent variables assessed ( Table 3). Doubling and halving the initial sorbed As concentration moderately influenced the total dissolved mass and the maximum dissolved As, but they had little effect on the location of the 10 µg/L front. By far, the model was most sensitive to the initial Fe(OH) 3 concentration. In particular, the model is sensitive to underestimations of initial Fe(OH) 3 . Because sorption to Fe(OH) 3 is the main retention mechanism for As, a poorly quantified initial concentration for Fe(OH) 3 can cause substantial errors in simulating As transport by either over-or underestimating sorption sites for As, which then either incorrectly enhances or lessens As retention, respectively. Collecting a robust background dataset will minimize errors that could result from natural heterogeneity and a small sample size. Table 3. Sensitivity coefficients for model independent variables (logK for surface complexation reactions, initial Fe(OH) 3 concentration, and initial sorbed As concentration) on the model dependent variables (location of the 10 µg/L As front, the total dissolved As mass, and the maximum dissolved As concentration) 37 years after the spill. Sensitivity coefficients were calculated relative to the base case using Equation (2). Concentration profiles of dissolved As for elevations of 422 m and 419 m above mean sea level (Figure 9) show the sensitivity of our model to alterations in intrinsic surface complexation logK values, initial Fe(OH) 3 , initial sorbed As, and two binary scenarios (i.e., "on" vs. "off"): (1) formation of HCO 3 − surface complexes, and (2) the incorporation of a commonly disregarded As(V) surface complexation reaction (Reaction 9 in Table 1).
The As concentration profiles show that doubling and halving the initial sorbed As and increasing and decreasing complexation logK values had a minimal effect on the distribution of As at 422 m and 419 m elevations (Figure 9), whereas changing the initial concentration of Fe(OH) 3 substantially impacts As concentrations. When the initial concentration is overestimated, the As plume becomes over-attenuated, and the dissolved As maximum at 422 m elevation occurs upgradient from the center of the oil compared to the base model, where the maximum occurs 65 m downgradient. Additionally, the maximum As concentration is an order of magnitude less than the base model and barely exceeds the 10 µg/L threshold. At 419 m elevation, dissolved As never exceeds the 10 µg/L threshold. This under-representation of dissolved As occurs because there is an excess of Fe(OH) 3 , and thus, an excess of surface sites to sorb As, which amplifies As attenuation. Underestimating initial Fe(OH) 3 has the opposite effect. At 422 m, the dissolved As maximum occurs at 145 m, 23 m further than the base model, though the As maximum concentration is roughly similar. At a depth of 419 m, the magnitude of As in groundwater is grossly overestimated (note log-scale). These trends occur because there is less Fe(OH) 3 , and thus fewer surface sites to sorb As, which allows for more persistent As in groundwater and further As transport. These two scenarios illustrate the critical importance of accurate background sediment chemistry, as an overestimation or underestimation of initial Fe(OH) 3 can cause substantial errors in predicting As mobility.
plex. When Reaction 9 is omitted, As is pervasive in groundwater beyond 200 m at both 422 m and 419 m elevations. Because As(III) sorption is the dominant attenuation mechanism in these geochemical conditions, omission of this reaction caused only slight deviations from the base model. However, we suspect model inaccuracies would be amplified if this reaction was omitted from model formulations for scenarios in which As(V) sorption was the main attenuation mechanism. , doubled initial sorbed As (pink, solid), halved initial sorbed As (pink, dashed), upper 95% confidence interval (C.I.) for intrinsic As surface complexation constants (green, solid), lower 95% C.I. for intrinsic As surface complexation constants (green, dashed), maximum measured (orange, solid) and minimum measured (orange, dashed) background Fe(OH)3 concentration, no complexation of HCO3 − (gold), and omission of an As(V) surface complexation reaction (Reaction 9 in Table 1) (blue). Field measured As concentrations shown in red circles are from [29]. Open circles designate samples <10 µg/L and are plotted at half the detection limit. Note that the y-axes are on the logarithmic scale. Sensitivity scenarios shown include the model base case (black), doubled initial sorbed As (pink, solid), halved initial sorbed As (pink, dashed), upper 95% confidence interval (C.I.) for intrinsic As surface complexation constants (green, solid), lower 95% C.I. for intrinsic As surface complexation constants (green, dashed), maximum measured (orange, solid) and minimum measured (orange, dashed) background Fe(OH) 3 concentration, no complexation of HCO 3 − (gold), and omission of an As(V) surface complexation reaction (Reaction 9 in Table 1) (blue). Field measured As concentrations shown in red circles are from [29]. Open circles designate samples <10 µg/L and are plotted at half the detection limit. Note that the y-axes are on the logarithmic scale.
The omission of HCO 3 − as a sorption competitor for As can lead to incorrect surface complexation modeling results for As sorption on HFO [62]. In our model, omission of HCO 3 − complexation increases dissolved As at 422 m and 419 m elevation ( Figure 9); the effect is particularly pronounced at 422 m, where maximum As concentrations are approximately 80 times greater than the base case. At both elevations, the location of the 10 µg/L front is approximately 45 m further downgradient compared to the base model. These results corroborate the findings of other studies that have identified sorption of HCO 3 − as a critical factor to accurately simulate As cycling in groundwater [20,61,62,71]. Some reactive transport models in the published literature that simulate As surface complexation on Fe(OH) 3 do not include the formation of the Hfo_wAsO 4 2− complex (Reaction 9 in Table 1), likely because a logK for this reaction does not exist in the Dzombak and Morel surface complexation database [59], which is the standard database for surface complexation modeling on HFO. However, based on thermodynamic data (Table 1), the dominant dissolved and sorbed As(V) species at the pH range of groundwater are HAsO 4 2− and Hfo_wAsO 4 2− , respectively (see Figure S1 for relative abundances of As aqueous phases and surface complexes vs. pH). This suggests that the addition of the logK value for this reaction may improve the simulation of As surface complexation in aquifers. For this reason, we included Reaction 9 and its corresponding logK reported in [61].
When Reaction 9 is omitted, there is no deviation in the modeled concentration profiles up to 140-150 m downgradient (Figure 9), likely because As(III) is the dominant sorbed form of As up to this distance. However, downgradient from 140-150 m, As is fully removed from groundwater in the base model due to formation of the Hfo_wAsO 4 2− complex. When Reaction 9 is omitted, As is pervasive in groundwater beyond 200 m at both 422 m and 419 m elevations. Because As(III) sorption is the dominant attenuation mechanism in these geochemical conditions, omission of this reaction caused only slight deviations from the base model. However, we suspect model inaccuracies would be amplified if this reaction was omitted from model formulations for scenarios in which As(V) sorption was the main attenuation mechanism.

Model Limitations
While the model used in this study accounts for biodegradation and secondary reactions that affect water quality, there are multiple limitations that could have confounding effects on As cycling.
The fate of dissolved Fe 2+ in the transition zone (~125-135 m downgradient from the oil at 37 years), remains unclear in the Bemidji aquifer. Sediment samples collected in the transition zone show spatially heterogeneous Fe(III) concentrations that generally correlate with particle size; sand-and-gravel-sized sediments have similar Fe(III) concentrations as uncontaminated background bulk sediments, whereas silt-sized sediments from heterogeneous lenses were substantially elevated in Fe(III) relative to uncontaminated bulk sediments [30]. Two different interpretations have been proposed to explain this phenomenon. The first is that the fine-grained sediments have low permeability and may be diffusion-controlled. If this is true, groundwater from the plume likely flows around these lenses, and they retain their initial geochemistry, which, compared to coarse sediments, tends to be more abundant in secondary mineral phases like iron oxides due to the larger surface area-to-volume ratio. This interpretation is supported by the observation of persistent Fe(III) in highly contaminated, fine-grained sediment near the oil [72], which suggests that these fine-grained materials may be recalcitrant to plume reactions. A second interpretation is that Fe 2+ oxidation and precipitation of Fe(OH) 3 occur heterogeneously in the transition zone, with the bulk of reactivity occurring in or adjacent to the fine-grained lenses, catalyzed by the more reactive surface area of fine-grained minerals. This interpretation is supported by the observation of low but detectable DO in the transition zone, generally from 100-1000 µg/L [29], which spans the range that can oxidize dissolved Fe 2+ and precipitate Fe(OH) 3 .
In its current formulation, the model does not produce the heterogeneous zones of elevated Fe(OH) 3 observed in sediments collected from the transition zone [30]. The postulated higher Fe(III) background signature of fine-grained sediments from the first interpretation was not incorporated into the model because there are insufficient data to constrain background Fe(III) concentrations in fine-grained sediments. Additionally, the model was not able to simulate the low DO concentrations observed in the transition zone, and therefore it does not precipitate Fe(OH) 3 in the transition zone as would be suggested from the second interpretation. Regardless of the mechanism that causes elevated Fe(III) in the transition zone, the mass of Fe(OH) 3 in the model's transition zone is likely an underestimation of the true mass in the aquifer, which would result in an underrepresentation of the aquifer's overall attenuation capacity for As.
The implementation of a homogeneous initial concentration of Fe(OH) 3 in the model is a simplification, as there is natural heterogeneity not just with respect to sediment grain size, but also with respect to the Fe(OH) 3 concentrations within the coarse fraction of sediments themselves. Based on sampling of background bulk sediment, Fe(OH) 3 concentrations can span a wide range (e.g., 775-2224 mg/kg in one background core [30]), though average background Fe(III) concentrations from two different locations were in close agreement, 1330 mg/kg [38] and 1390 mg/kg [30]. By invoking a homogeneous distribution of Fe(OH) 3 as an initial condition, the model likely neglects some of the finer-scaled nuances in Fe cycling caused by natural heterogeneity, which could affect As cycling.
Given the modeling result that the prevalence of As in groundwater is controlled by the presence of NVDOC, a discussion of the constraints of NVDOC biodegradation is warranted. NVDOC is a complex mixture of organic C of varying chemical structure and bioavailability. Ng et al. [28] and this study generically represent NVDOC as C 19 H 24 O 6 and characterize its biodegradation by a first-order rate constant calibrated to measurements from two observation periods. However, given the actual complexity of NVDOC, its biodegradation may follow a different kinetic rate model than that used in the model, which may influence how long As remains in groundwater.

Human Health Consequences
To assess the relative threats of benzene and As in groundwater over time, we compared changes in their respective plume volumes over the model simulation period (Figure 10). Because benzene, ethylbenzene, and xylene are modeled under a common parameter (BEX), we assumed that 90% of the BEX mixture was benzene based on previous groundwater saturation experiments [73]. In 2013, field observations showed that the plume of benzene exceeding its MCL (5 µg/L) extended~15 m further downgradient than the plume of As exceeding its MCL (10 µg/L). Otherwise, the plumes were roughly similar in shape and size [29]. Long-term model simulations show that the benzene plume remains relatively stable and eventually retracts toward the oil source. Meanwhile, the As plume continues to migrate downgradient. Generally, the volume of aquifer exceeding the As MCL is greater than the aquifer volume exceeding the benzene MCL (Figure 10). At its maximum (100 years),~5.7× more groundwater volume exceeds the As MCL than the benzene MCL. Based on these model results, it is possible that As could contaminate more groundwater over a longer period of time than benzene in the Bemidji aquifer.
Water 2021, 13, x FOR PEER REVIEW 19 of 24

Human Health Consequences
To assess the relative threats of benzene and As in groundwater over time, we compared changes in their respective plume volumes over the model simulation period (Figure 10). Because benzene, ethylbenzene, and xylene are modeled under a common parameter (BEX), we assumed that 90% of the BEX mixture was benzene based on previous groundwater saturation experiments [73]. In 2013, field observations showed that the plume of benzene exceeding its MCL (5 µg/L) extended ~15 m further downgradient than the plume of As exceeding its MCL (10 µg/L). Otherwise, the plumes were roughly similar in shape and size [29]. Long-term model simulations show that the benzene plume remains relatively stable and eventually retracts toward the oil source. Meanwhile, the As plume continues to migrate downgradient. Generally, the volume of aquifer exceeding the As MCL is greater than the aquifer volume exceeding the benzene MCL (Figure 10). At its maximum (100 years), ~5.7× more groundwater volume exceeds the As MCL than the benzene MCL. Based on these model results, it is possible that As could contaminate more groundwater over a longer period of time than benzene in the Bemidji aquifer. This modeling study suggests that over the long-term, the legacy of As in groundwater may become a greater threat than petroleum hydrocarbons from a human health perspective for three reasons: (1) simulated As is transported further downgradient and contaminates groundwater further from the oil source than petroleum hydrocarbons, which biodegrade and retract to the oil source; (2) the aquifer volume that exceeds the As MCL This modeling study suggests that over the long-term, the legacy of As in groundwater may become a greater threat than petroleum hydrocarbons from a human health perspective for three reasons: (1) simulated As is transported further downgradient and contaminates groundwater further from the oil source than petroleum hydrocarbons, which biodegrade and retract to the oil source; (2) the aquifer volume that exceeds the As MCL is greater than the aquifer volume that exceeds the benzene MCL over the 400-year simulation period; and (3) after biodegradation, petroleum hydrocarbons are fully removed from the aquifer, whereas sorbed As can be remobilized if a new source of biodegradable organic carbon is introduced. Regulatory agencies are charged with monitoring petroleum hydrocarbons, including benzene, in groundwater near oil spill sites, but As is not generally considered as a contaminant of concern at these sites, and therefore, agencies may not monitor for geogenic As mobilization as a consequence of petroleum hydrocarbon biodegradation. In effect, this might allow As to routinely go undetected at oil spill sites.
Arsenic is widespread in aquifer sediments throughout the U.S. and worldwide [5,74]. Recent calculations show that essentially any aquifer has enough natural As to cause unsafe concentrations in drinking water if geochemical conditions promote As release from aquifer sediments to groundwater [31]. Although few other studies have quantitatively examined As fate in petroleum hydrocarbon plumes, several studies have found that As is released from aquifers under enhanced bioremediation conditions where organic carbon, such as whey, is added to stimulate reduction of chlorinated hydrocarbons [75][76][77][78]. These results, combined with our modeling results, underscore the importance of examining the long-term water quality and health consequences of As at sites with biodegradable organic carbon, including, but not limited to, sites with oil spills, injections of carbon for enhanced bioremediation, landfills, sewage plumes.

Conclusions and Opportunities for Future Research
Using a reactive transport model, we simulated the release of geogenic As that occurs during Fe-reduction coupled with the biodegradation of petroleum hydrocarbons resulting from a crude-oil spill near Bemidji, Minnesota (USA). The reactive transport model was informed by decades of field observations for organic and inorganic constituents from the petroleum-contaminated aquifer. Our modeling results suggest that As concentrations in groundwater within the plume will exceed 10 µg/L for~400 years. Anoxic conditions are maintained due to oxygen demand exerted by biodegradable organic carbon in the plume. However, after the organic carbon degrades and oxygenated conditions are re-established, As will re-associate with Fe(OH) 3 . Our simulations show that As will contaminate more groundwater than benzene from the oil spill, suggesting that As should be closely monitored in petroleum-contaminated aquifers. Investigations into the toxicological effects of combined As and benzene exposure are warranted, as this study suggests the two can co-occur in petroleum-contaminated aquifers for centuries. Such a study would provide crucial knowledge of the human health risks at oil spill sites.