Towards Understanding Factors Affecting Arsenic, Chromium, and Vanadium Mobility in the Subsurface

Arsenic (As), chromium (Cr), and vanadium (V) are naturally occurring, redox-active elements that can become human health hazards when they are released from aquifer substrates into groundwater that may be used as domestic or irrigation source. As such, there is a need to develop incisive conceptual and quantitative models of the geochemistry and transport of potentially hazardous elements to assess risk and facilitate interventions. However, understanding the complexity and heterogeneous subsurface environment requires knowledge of solid-phase minerals, hydrologic movement, aerobic and anaerobic environments, microbial interactions, and complicated chemical kinetics. Here, we examine the relevant geochemical and hydrological information about the release and transport of potentially hazardous geogenic contaminants, specifically As, Cr, and V, as well as the potential challenges in developing a robust understanding of their behavior in the subsurface. We explore the development of geochemical models, illustrate how they can be utilized, and describe the gaps in knowledge that exist in translating subsurface conditions into numerical models, as well as provide an outlook on future research needs and developments.


Introduction
Potentially toxic metal(loid)s in groundwater pose a significant threat to human health. A number of potentially toxic elements can be released from geomaterials into groundwater, including arsenic (As), chromium (Cr), vanadium (V), selenium (Se), fluorine (F), molybdenum (Mo), manganese (Mn), and boron (B) [1][2][3][4][5][6][7][8][9][10]. The most well-known geogenic contamination of groundwater occurs in Southern Asia, where up to 75 million people are estimated to be at risk of increased disease burden due to drinking-wells drawing from a shallow aquifer containing high As concentrations, a situation that has been described as "the largest mass poisoning in history" [11]. Worldwide, more than 150 million people drink water from wells with arsenic (As) concentrations greater than the USEPA and WHO maximum level of 10 μg/L [12,13]. Arsenic is a well-documented carcinogen, affecting the skin, liver, kidneys, lungs, and bladder. Chronic exposure is frequently associated with skin lesions, cardiovascular and respiratory disorders, reproductive failures, insulin resistance, and neurotoxicity [12][13][14]. Although As is often used as an example of an element that exhibits chronic toxicity and is globally widespread, related oxyanion-forming trace elements (e.g., Cr and V [15][16][17][18], which, along with As, we focus on in this work), are also naturally elevated in specific regions across the globe. For example, Cr (as hexavalent Cr) is a dangerous carcinogen, teratogen, and mutagen that has been found to occur naturally in groundwater systems, including in the continental US [7,19]; although the toxicology of V is less understood [20], this element can be very mobile under naturally occurring conditions [21] and is emerging as an environmental concern [22]. Despite the catastrophic size of the problem, global sustainable development goals do not consider As [23], or other geogenic contaminants, as a factor in supplying drinking water to communities worldwide.
Drinking water is a primary route of exposure for As, often from groundwater, which more than 2 billion people worldwide use as their primary domestic water source [24,25]. This dependence, along with the fact that >99% of accessible freshwater resides in groundwater [26], makes understanding, predicting, and managing of groundwater of foremost importance. However, subsurface processes are notoriously poorly understood, and controlled by coupled hydrogeochemical and biogeochemical processes, often involving kinetic or transport limitations, microbial processes, and interfacial reactions. In general, the solubilization of the potentially toxic contaminants in aquifers, be they geogenic or anthropogenic in origin, are thought to be controlled by phase transformations that are mediated thought redox, precipitation, dissolution, and (de)sorption reactions, all influenced by pH, ionic strength, and the presence of other constituents [3,5,9,16,18]. Beyond chemical reactions, physical considerations, such as lithologic heterogeneity or hydrologic complexity, limit our understanding of subsurface processes impacting contaminant behavior.
These complex physicochemical factors limit our ability to produce comprehensive and robust conceptual models, and hence reliable quantitative models, that represent our understanding of geogenic contamination in the subsurface and guide management to protect human health. The goal of this work is to review the current state of knowledge, identify key challenges, and provide an outlook for resolving issues associated with assessing and modeling inorganic geogenic contaminants, including As, Cr, and V, in subsurface environments that are sources of drinking water. This article aims to: (1) consider the groundwater hydrology and its effect on transport and distribution; (2) discuss geogenic sources of these contaminants; (3) investigate geochemical mechanisms affecting contrasting contaminant behavior of metal(loid)s; (4) describe methods for determining the concentration and speciation of metal(loid)s in water and geomaterials; (5) discuss approaches for deriving mechanistic models of contaminant transport and fate; and (6) conclude with a reflection on our needs and an outlook for future activities. Understanding and creating conceptual models of the release, retention, and transport of geogenic contaminants in groundwater systems, will ultimately enable the creation of quantitative models that can be used to predict future states and guide efforts to remediate contaminated groundwater to protect human health.

Groundwater Flow and Its Relationship to Chemical Evolution
To understand the chemical release and mobility of contaminants in aquifers, it is first necessary to develop conceptual models of groundwater flow and how physical factors such as flow velocity, mixing, and topography influence its geochemical evolution ( Figure 1). Different scenarios are needed to understand how the master variables Eh and pH change over space and time, as these parameters often exert dominant control over the geochemical processes that result in contaminant solubilization or sequestration in the subsurface.
Groundwater in confined and unconfined aquifers occupies distinct environments for oxicanoxic relationships and trace element behavior. In unconfined aquifers, the water table is an interface between saturated and unsaturated zones where inputs of fresh oxic water and soil gases can occur along the flow path. Therefore, as an interface between saturated and unsaturated conditions, the water table is also an interface that modulates the redox state and pH of shallow groundwater via open system inputs of oxygen (O 2 ) and carbon dioxide (CO 2 ). However, anoxic groundwater is not unusual in unconfined aquifers where electron donors such as organic carbon can scavenge O 2 [28][29][30]. In contrast, the upper interface of a confined aquifer is defined by a relatively impermeable formation, rather than the water table. Therefore, age and redox conditions in confined aquifers are associated with distance along the flow path, with groundwater becoming more reducing away from the recharge zone as it encounters electron donors.
Groundwater chemistry evolves in terms of pH, redox state, and major ion concentrations. From the water table, recharging waters are predominantly oxic, with acidic pH and low alkalinity [31,32] but become more alkaline and reducing after encountering subsurface materials (including organic matter or reduced mineral phases), which in turn impact the speciation and movement of As, Cr and V ( Figure 2). However, specific flow paths may be complex, incorporating mixtures of groundwater from different sources.
The topography of a landscape impacts groundwater flow patterns and may also have some association with groundwater redox gradients. Where confining units are absent, groundwater recharge augments flow along a flowpath. Accordingly, waters at mid-to high elevations on a landscape are dominated by recharge water whereas waters at mid-to low elevations on a landscape represent a mix between recharge and existing flowpath. Increasing distance blow the water table can promote development of more suboxic to anoxic and more alkaline conditions, as reductants and other reactive minerals are encountered ( Figure 1). In discharge zones characterized by low topography in unconfined aquifers, groundwater flow can be significantly upward-directed, bringing more evolved, anoxic waters closer to the land surface [29].
Oxic and anoxic groundwater can be separated by centimeters to tens of meters vertically -capable of mixing across the scale of a well screen-and meters to hundreds of meters along a flow path (although key microscale redox variation may occur at the pore scale [37,38]). Groundwater evolves toward redox interfaces both vertically where groundwater flows predominantly downward (e.g., hilltop areas) and laterally where groundwater predominantly flows toward local stream valleys. These trends suggest the presence of the oxic-anoxic interface as a third critical interface along with the regolith-bedrock interface and the water table, which together contribute to landscape-scale heterogeneity of groundwater chemical conditions [15,39]. This pattern also suggests that depth and topographic position (hilltop vs. valley) could be master variables for approximating groundwater chemical evolution in an unconfined aquifer [40]. In addition to flow direction, there is nominal evidence that groundwater flow velocity can impact the concentration of redox constituents in groundwater; for example, a study in a shallow aquifer in Greece found that slower groundwater flow was correlated with higher concentrations of Cr(IV) due to a longer period of water-rock interaction [41].
Although the geology and hydrology of an area sets up the transition between oxic surficial inputs and electron donors hosted in the aquifer and forms the framework for groundwater redox evolution, anthropogenic inputs can alter groundwater redox evolution. Mixing can introduce oxic waters into a reduced environment, or reduced waters into an oxic environment, altering the overall chemistry. In pumped aquifers, mixing of disparate redox states can occur across long well screens [28]. Additional sources of electron acceptors (e.g., O 2 or nitrate (NO 3 − )) or donors can result from surficial inputs by human land use activities.
These can include organic materials, fertilizers, inorganic contaminants, petroleum, etc., and each will impart their own effect on the redox environment by their reactions and interactions with aqueous substrates, microorganisms, and existing groundwater conditions [42][43][44][45].

Limitations and Challenges
Recent studies on occurrence and health effects suggest that a large number of wells and aquifers are of concern for As or another naturally occurring element [12]. This situation is suited to probabalistic modeling approaches (e.g., probability of exceeding a water quality standard [15]), but mechanistic interpretation and modeling that could inform prediction at the scale of a water well are significantly more challenging. Several gaps in knowledge and practice are apparent for translating element occurrence to improved understanding of expected trace element behavior in well water. First, solid-phase geochemistry and mineralogy may only be known from a few sites in a region that contains many wells.
In practice, the solid-phase level of an element of interest might be represented only approximately, for example by a map-scale geologic unit in which a well is located. Second, vertical heterogeneity can be considerable at the scale of a well's open interval, which would create difficulty for translating element mobility concepts derived at the microscale/pore scale to the redox evolution that occurs along the scale of a realistic flowpath. Third, groundwater in a shallow unconfined aquifer represents a mix of multiple water sources of different residence times and flow lengths. This problem of along-flow and vertical heterogeneity is apparent in studies showing spatial variability of As, Cr, or other elements across the landscape, as discussed above. Heterogeneity's effects on water chemistry might be especially significant when analyzing a pumped water sample across a large open interval, which can induce mixing of different water sources and chemistries. Groundwater flowpaths and mixing relationships can be inferred broadly from regional and topographic generalizations [40], but the flow path or mixing relationships of a pumped water sample may not be known for wells outside of well-characterized research sites. Collectively, these factors imply that the multiple geochemical parameters within a water sample (e.g., redox-sensitive solutes, trace element concentrations, speciation) are encoded with richer geochemical process information than the other more general attributes in field-based studies, such as well depth, landscape position, map-scale geological unit, or solid-phase composition from limited solids analysis. Therefore, there is a gap between understanding the cumulative net geochemical processes that resulted in a water sample's observed chemistry, and predicting where in an aquifer, and at what rate these effects are imparted on groundwater.

Mineralogical Sources of Geogenic Contamination
Anthropogenic activities, such as electroplating, metallurgy, chemical production, energy production, pest control, and waste disposal, are common activities that result in contamination by metal(loid)s [46]. However, in the last two decades, there has been a growing recognition of the widespread importance of geogenic sources of these contaminants in groundwater [8]. This evolution in thinking is critical; a relatively low concentration of an element in geomaterial has the potential to contaminate groundwater for thousands of years under normal conditions [9]. However, the relative mobility of these elements is strongly controlled by what phases harbor contaminants, as well as the water chemistry of aquifer. Mineralogical information has not historically been considered by regulators, who are primarily concerned with the composition of the groundwater, and the presence of aqueous-phase hazards [47]. Generally, aquifer solids consist of multiple minerals, each containing different levels of As, Cr, or V, and having a wide range of weathering rates. We thus discuss the mineralogical sources of geogenic As, Cr, and V as a prelude to reviewing biogeochemical processes associated with their solubilization and transport. It should be noted that the factors that solubilize the elements may also redistribute them to more dilute, labile solid phase associations, such as adsorbed to clay mineral or (oxhydr)oxide phases, resulting in potentially new, potentially diffuse sources of contamination [48].

Arsenic
As is a ubiquitous element in the Earth's crust and the 20th most abundant; concentrations in igneous and sedimentary rocks are typically 2-3 mg/kg, but concentrations can be higher in clayey sediments, phosphorites, and reduced marine sediments [49]. Geogenic As is often hosted mainly by mineral phases in metavolcanics, mudstones and their alteration products. There are over 500 As bearing minerals, including sulfides, oxides, arsenates and arsenites [50]. Because As can be coprecipitated with iron (Fe) (oxyhydr)oxide and sulfides in sedimentary rocks, and sedimentary Fe ores and Mn nodules are considered As-rich materials [49]. The most common and most studied are arsenopyrite (FeAsS) and As-bearing pyrite (FeS 2 ), where As can substitute for S in the pyrite crystal structure. Globally, As rich ore deposits are found in China and Morocco, as well as the US, Russia, Japan, Iran and Belgium [51,52].
Arsenic contamination to groundwater derives not just from mineralized deposits, but also from hydrothermal activity (as could be found in Yellowstone, USA) and from aquifer sediments associated with silicic volcanic aquifers, and tuffaceous sediments and in sedimentary basins derived from these volcanic rock types [4,[52][53][54]. The conditions under which As in released from minerals into a groundwater system are varied based on the source material. Thus, contamination that results from reducing groundwater environments has occurred in Northern China, Taiwan, Vietnam, Hungary, Romania, Slovakia, Bangladesh, southwestern US, and India. Conversely, in Mexico, Chile, Argentina, and the southwestern US, arid and oxidizing environments trigger the release of As to groundwater, highlighting the need to understand the biogeochemical processes associated with contamination in an aquifer [51,52].

Chromium
Cr, the 21st most abundant element in the Earth's crust [55], occurs at an average concentration of around 100 mg/kg [41], typically as oxides and silicates [7,19]. In rock forming minerals, Cr pyroxenes, amphiboles, diopside, serpentine and chlorite are the chief hosts of Cr, where up to 1.5% Cr 2 O 3 have been reported [56,57]. In Cr-bearing spinels, commonly chromite (FeCr 2 O 4 ), magnesiochromite (MgCr 2 O 4 ), and magnetite (Fe 3 O 4 ), the percentage Cr could be up to 35%, often presenting as Cr 2 O 3 [41,56]. Alteration of primary silicate minerals during metamorphic events may produce minerals such as fuchsite, uvarovite, Cr-mica and tawmawite, where Cr substitutes for Al in the octahedral site through isomorphism [19]. Globally, chromite is the most common deposit of Cr, found circumpacific and the Mediterranian [7] as well in South Africa, Khazahstan, Zimbabwe, and India [55]. It is generally associated with ultramafic rocks and their hydrothermally altered equivalents [19,41]. Though considered to be of low solubility and relatively chemically inert, dissolution of a small fraction of chromite can release elevated levels of Cr.
There have been reports of elevated levels of Cr in groundwater systems worldwide, most at levels exceeding the 50 μg/L maximum contaminant level (MCL) set by the World Health Organization (WHO). Concentrations up to 700 μg/L have been reported in soil solutions in phosphate-fertilized soils in Greece, but geogenic Cr levels around the world are usually at least an order of magnitude below this value [56,58]. In the US, high concentrations of total Cr, mostly occurring as Cr(VI) have been detected at levels up to 60 μg/L in the Mojave Desert [59]. A combination of pH as high as 9 and dissolved O 2 > 1 mg/L and fluctuating redox conditions dictated Cr speciation in groundwater sourced from alluvial aquifer materials weathered from mafic rocks in such areas [59].

Vanadium
Vanadium occurs in the Earth's crust at an average content of around 2 mg/kg and is generally widely dispersed; concentrated mineral deposits are rare [60]. V tends to be present in large quantities in mafic igneous rocks, in lesser quantities in ultramafic rocks, and in very small quantities in felsic rocks [61]. There are more than 80 V-bearing minerals present in bedrock materials occurring as sulfides, sulfates, silicates, oxides, phosphates and vanadates [21]. Due to its affinity for organic rich sediments, V is also found in high concentrations in crude petroleum and shales [62]; an average concentration of 790 mg/kg have been reported in black shales around the world [63]. Major V deposits tend to occur as sandstone-hosted V, shale-hosted V, vanadate, or vanadiferous titanomagnetite, occurring in China, Russia and South Africa, as well as the US, Canada, and Australia [61]. Dispersion of V in the environment occurs due to its sensitivity to O 2 and acidity, as well as its ability to readily substitute for both Fe and Al in both primary and secondary materials, such as clays, due to the size and charge of its trivalent ion in the environment [62]. Anthropogenic activities such as steel production, coal, petroleum production, rock phosphate fertilizers, superphosphate fertilizers and mining are also major anthropogenic sources of V into the environment [64].
In groundwater, the occurrence of V is highly dependent not only on potential anthropogenic sources, but also on the geology and aqueous conditions, and as such will be highly variable. For example, the average V concentration in California groundwater was found to be 5 μg/L. However, up to 70 μg/L V have been reported in groundwater in the northeastern San Joaquin Valley [5,6]. In some drinking water sources in Italy, values as high as 180 μg/L have been detected [5,65]. In general, higher concentrations in groundwater occur where the aqueous conditions are either alkaline, oxic, or both. [6]. High groundwater V concentrations have been documented in oxic as well as suboxic, near-neutral groundwater [6].

Challenges
Despite a sophisticated understanding of the mineralogy associated with Cr, V, and As, there are many challenges associated with understanding the sources of subsurface contamination with these elements. Although elements of interest can be broadly dispersed on aquifer rocks throughout geologic unit or region (as is often the case for As, Cr, and V), the elemental distribution or solid-phase speciation may be quite heterogeneous at a small scale, pointing to "hot spots" where elements are concentrated in the solid phase and become available for release to groundwater. It is worth noting that all Cr, V, and As, though ultimately derived from the sources above, are often associated as secondary deposition with Fe or Al (oxyhydr)oxides [19,49,62], and may be influenced by the distribution of solid oxidants and reductants, like Mn oxides and organic matter [66,67]. As a further complication, trace element release from (oxyhydr)oxides can be dominated by desorption, which is much more rapid than mineral dissolution, resulting in a larger element release rate relative to groundwater flow rate through the aquifer solids. The distribution of these elements in primary and secondary minerals may vary over short horizontal or vertical differences in spatial location. Such a spatially variable distribution results in an uncertain source region that is in the subsurface is the norm, not the exception, for geogenic contaminants. Consequently, most of the required information for developing a quantitative understanding of sources and mobilization mechanisms is initially laborious and site specific, whereas a region-wide synthesis to reduce population-scale exposure should evaluate and encompass a range of possibilities both in terms of solid-phase contaminant occurrence and the hydrogeochemical mechanisms that control contaminant levels in groundwater.

Biogeochemical Processes
Biogeochemical processes can couple with the hydrological and mineralogical factors above to result in the solubilization of geogenic As, Cr, and V. The evolution of groundwater composition (viz. Eh and pH, but also the dissolution other solutes) is driven by a complex array of transport and chemical reactions that together dictate the interactions between the major solid (e.g., Fe (oxyhydr)oxides and silicate minerals) and dissolved phases within the aquifer. This evolution also drives the solubilization and mobilization of lower abundance species, including geogenic contaminants. It is worth noting that contaminants often have contrasting redox controls on their solubility and mobility, meaning that conditions that favor mobilization of one element may favor sequestration of another.

Redox Controls on Speciation
A predominant factor that controls the speciation therefore the solubility, mobility, and bioavailability, of As, Cr, and V is the oxidation state of the elements and pH ( Figure 3). For these elements, a variety of oxidation states is possible under natural aqueous conditions (slightly acidic to slightly alkaline in most settings), and the different states result in species that can have very different solubilities. The similarities in the solubilizing conditions for Cr and V and the opposite solubilizing conditions for As help to explain the concurrences and anti-correlations of certain elements in water supplies [15,16,68].
Inorganic As, under natural conditions, is typically found as either As(III) or As(V). Pentavalent arsenic occurs predominantly as oxyanions at near-neutral to alkaline pH, whereas As(III) commonly occurs as the uncharged As(OH) 3 . In comparison to other oxyanions, including Cr and V, As retains more mobility across a wider range of groundwater pH and redox conditions, particularly under mildly reducing conditions where As(III) predominates as soluble As(OH) 3 [70].
Chromium typically occurs in one of two oxidation states, Cr(III) and Cr(VI). In contrast to As, in which the reduced state is more soluble, Cr(VI) is known to be significantly more mobile than Cr(III), and thus Cr(VI) is more likely to be found in groundwater when conditions are oxidizing and alkaline [71].
Vanadium can occur under typical groundwater redox conditions as two different oxidation states: V(IV), dominant under low pH and suboxic conditions, and V(V), which dominates under high pH or oxic conditions. Both oxidation states can be soluble and can co-occur depending on the pH and co-occurring constituents of the groundwater [21]. However, V(IV) has been much less documented than V(V) in groundwater studies, V(IV) has been inferred based on indirect evidence [6].

Redox Transformations
Fluctuations in redox conditions, or the presence of other redox-active materials, can result in redox reactions that change the oxidation state of As, Cr, or V from a mobile form to an immobile form, or vice versa. The presence and concentration of redox-active elements can have an outsized effect on the redox potential of the groundwater system, which in turn affects the redox state and speciation of trace elements. Geologic or soil systems can have high concentrations of Fe and Mn [9,72], which can act as a redox buffer due to their relatively high abundance compared to other redox-sensitive elements, poising the system to a steady potential, and altering the transformations, and therefore the solubility, of trace elements such as a As, Cr or V [43,73].
Iron and Mn, in dissolved or solid forms, may also directly participate directly in redox reactions with Cr, V, and As ( Figure 4). Manganese oxides, a primary subsurface oxidant, will oxidize As(III) to the less mobile state of As(V), a reaction that occurs more quickly, and therefore more favorably, than oxidation by Fe (oxyhydr)oxides [72]. These reactions can take seconds to days and are generally influenced by concentration of As(III), temperature, and pH. Clays can also act as reaction surfaces for the oxidation of As(III) to As(V), again, strongly influenced by pH, ionic strength, and composition of clay surfaces [49]. Likewise, Cr(III) and V(IV) can also be oxidized to Cr(VI) and V(V) by naturally occurring mineral oxidants, most commonly Mn oxides [74,75]. For Cr(III), the reaction with Mn oxides is dependent on the pH of the environment; where the pH is less than 6, the reaction occurs with dissolved Cr(III) adsorbing onto the surface of solid Mn oxides, whereas for a pH greater than 8 the reaction occurs by the redox cycling of Mn the surfaces of solid Cr(III) surfaces [74].
Reduction reactions are also significant. Cr(VI) can undergo reduction to its more insoluble Cr(III) form via redox reactions with dissolved Fe(II), Fe(II) minerals (pH > 5.5), sulfides (pH < 5.5), and organic matter [74,76]. Similar to reactions with As and Cr, Fe(II)-bearing minerals have also been suggested to actively reduce V(V) to V(III), another oxidation state that typically is only found in anoxic systems, though this reaction occurs over the course of months [77].
Although redox transformations controlling the solubility of elements can occur abiotically, microorganisms often mediate redox reactions because they utilize these transformations for energy synthesis and/or the breakdown of organic material [3,43]. As such, redox transformations tend to occur in a thermodynamically predicted sequence that corresponds to a progressively more reduced redox potential, with microbes preferentially utilizing reactions that will provide the most energy for their metabolism. This preferential order and the associated redox potential of each transformation is known as the redox ladder [78]. In organic-rich aquifers, the predominant electron donor driving microbially-mediated redox reactions is organic matter. Microbially mediated reactions can directly drive the redox environment of the groundwater system by consuming electron acceptors (O 2 , nitrate, Mn oxides, Fe (oxyhydr)oxides, and sulfate), sustaining anoxic conditions in groundwater separated from the water table, and farther along flowpaths promoting reducing conditions characterized by detectable Mn 2+ , Fe 2+ , sulfide, or methane.
Along this redox trajectory, As, Cr, and V also undergo oxidation state changes ( Figure  2). In addition to interactions with abiotic species, specific organisms also oxidize or reduce As, Cr, or V, indicating that trace elements can be converted through microbiallymediated processes (in addition to overall microbially-mediated control of the redox state of groundwater or sediments). In some cases, this can lead to redox-active contaminants being utilized by microbes directly as an energy source. Vanadium, having a high reduction potential and a relatively low bacterial toxicity, has been shown to act as an electron acceptor for the respiration of a number of microbial species, including Pseudomonas vanadiumreductans T-1, Pseudomonas isachenkovii A-1, and Geobacter metallireducens, among others [62,79]. Arsenic-metabolizing species have also been found in both contaminated and non-contaminated sediments [79,80]. In addition to oxidation and reduction, detoxification of As via methylation is a widespread bacteria capability [80,81]. Microbial interactions with Cr tend towards detoxification mechanisms as well; Cr(VI) adsorption by the cell membrane and subsequent reduction to Cr(III) has been shown to occur in both fungi and bacteria [81].
Microbes can also influence the solubility and mobility for redox-active elements by mediating the formation of other materials in the groundwater system. Microbially created amorphous Fe and Mn (oxyhydr)oxides can directly contribute to the reactions that affect the redox potential of the groundwater system. In turn, these metal (oxyhydr)oxides can be utilized by microbes for their metabolic activities, releasing materials that are bound to the (oxyhydr)oxide surfaces, i.e., reductive dissolution. Because of widespread association of As with Fe (oxyhydr)oxides, reductive dissolution is a common mechanism of As release into groundwater, where the decay of organic matter consumes O 2 , creating anoxic conditions that promote reduction of the Fe (oxyhydr)oxide and As to the more soluble As(III) form.
Similarly, V will strongly sorb to Fe (oxyhydr)oxides, and can undergo reductive dissolution in a manner analogous to As; this process can lead to an increase in both dissolved V and dissolved Fe under reduced conditions [64]. The oxidation of sulfide minerals to Fe (oxyhydr)oxides can also concomitantly release As and sulfate into the groundwater environment [13,49]. Microbes will create organic materials that can change the mobility of As, Cr, and V by acting as ligands, or by directly reducing or oxidizing the element; the redox-active element can also be released into groundwater systems when the organic materials are degraded for microbial metabolism.

Sorption and Desorption Reactions
In addition to redox reactions, interfacial interactions such as adsorption and desorption reactions, can affect contaminant mobility. Both As(III) and As(V) will sorb to Fe and Al (oxyhydr)oxides (particularly high surface area or poorly ordered minerals [82]), depending on the pH and competing ions. In addition to oxidizing As(III), as noted above, both Mn oxides and clays will also sorb As [49,72,83], further limiting the mobility of As.
Generally, As(III) is thought to bind weakly as an outer sphere complex whereas As(V) is more likely to form strong inner sphere complexes [84]. Sorption of V is influenced also by its differently charged oxidation state. Although both the anions and cations will sorb to mineral surfaces over a wide range of pH conditions, V(IV) cations will sorb more strongly to the negatively charged (oxyhydr)oxide surfaces than the V(V) anions [5], but this is pH dependent and V(V) is known to form inner sphere complexes [84]. Adsorption and desorption to aquifer materials have been shown to be the primary control on V concentrations where groundwater conditions are oxic [6]. In contrast to As and V, the limiting factor of Cr mobility is the weathering of Cr(III) materials. Weathering enables a small amount to dissolve into the aqueous phase, where it can be adsorbed onto clay minerals or precipitated out of solution with either with (oxyhydr)oxide (Al and Fe) or with free cations (Fe(III) and Mn(II) [74,76]. In acidic environments (pH < 6), Cr(IV) oxyanions can become adsorbed by cationic colloids (usually as an inner-sphere complex [84]), an effect enhanced with increasing total organic carbon content in soil or sediment [74].
In addition to adsorption and desorption, the effect of competition between different ions can affect mobility via surfaces preferentially adsorbing one ion over another. This has been studied with As in particular, as it is ubiquitous and toxic to humans, making its sorption chemistry of particular interest. Of the common ions that appear in environmental aqueous environments, phosphate and carbonate oxyanions stand out as being competitive with As(V) for binding to the surfaces of Fe (oxyhydr)oxides [85][86][87]. When phosphate and carbonate are present in the groundwater system, they may be able to outcompete arsenate for inner-sphere complexation to the surface, resulting in more As in solution [85,86]. Sulfate ions also shows competitive behavior, but to a lesser extent and at least partially dependent on pH, suggesting outer-sphere complexation [85]. Coexisting major ions also modify adsorption sites, in addition to competing with As for adsorption sites. For example, sodium modifies the surface charge of As adsorption sites, resulting in higher As concentrations being found in Na-dominated groundwater rather than Cadominated groundwater. Na-dominated major element chemistry is common to groundwater in volcanic-source basins [4,88,89]. Dissolved organic matter, particularly humic and fulvic acids, have also been found to compete with As for mineral surface binding sites. How effectively these organic acids inhibit As adsorption has been shown to be dependent on the pH of the aqueous phase [90]. However, organic materials can also form ternary complexes with arsenate via an intermediary bridging cations, such as Fe(III), which can favor retention of As rather than release [91]. Sulfate and phosphate have also been shown to compete with Cr for adsorption sites [92], V has been shown to have differing competitive behaviors, specifically coordinating with biogenic ligand and compounds, depending on its oxidation state. Its anionic form resembles and competes with phosphate, whereas its cationic form more closely resembles the behavior of transition metal ions [93].

Challenges
From the sections above, it is clear that a large set of biotic and abiotic reactions control the evolution of groundwater chemistry, and thus the speciation and solubilization of As, Cr, and V. This, along with our understanding of source mineralogy and groundwater flow, allows formulation of a multiscale conceptual model ( Figure 1) in which groundwater becomes contaminated with geogenic contaminants as it travels through the subsurface. In general, we have a reasonable understanding of the reactions that are important to this evolution; however, in many cases, the kinetics of these reactions are unknown or at least poorly constrained. This is made even more difficult by noting that many of the processes are: (1) strongly Eh and pH dependent; (2) involve not only oxyanionic but sometimes cationic (e.g., V(IV)) or uncharged (e.g., As(III)) species; (3) heterogenous reactions involving interactions between aquatic species, minerals, and/or solid-phase organics; (4) microbially mediated processes; and (5) kinetically controlled in the natural setting (that is, slow and may be observed at disequilibrium), which severely limits quantitative approaches to modeling of geogenic contamination of groundwater.

Measurement Methods
Critical to the study of metal(loid) distribution and fate in the natural environment is being able to determine concentration and speciation of the metal(loid) in both the aqueousphase ground water and the solid-phase aquifer substrates. When groundwater samples are analyzed for the total concentration of metal(loid)s, there are number of standard procedures (in the US these are often proscribed by regulatory agencies) that typically involve an instrumental method using inductively coupled plasma (ICP) to vaporize and excite elements coupled to deception by optical emission spectrometry (OES) or mass spectrometry (MS). Detection by ICP-OES utilizes individual elements specific interactions with electromagnetic radiation whereas ICP-MS determines the concentration of an element in a sample by converting elements into gas phase ions and using a series of electromagnets to sort them by their mass-to-charge ratio. Selecting for the appropriate method depends on the sensitively needed, and the occurrence of analytical interferences [94]. Determining the speciation and oxidation state, of metal(loid)s in water can be conducted by coupling one of the previous analytical methods with ion chromatography, either in-line with the analytical instrument [95][96][97] or utilizing field-based separation techniques [98]. It should be noted that for some species (viz. Cr(VI)), specific approaches may be required to obtain adequate sensitivity quantitation in groundwater [99].
Total concentrations of metal(loid)s in aquifer solids can be achieved in several ways, depending on the elements of interest and their abundances. X-ray Fluorescence (XRF) spectroscopy utilizes high-energy X-rays to excite electrons in atoms at the surface of a solid-phase sample, and emitted fluorescence radiation is used to determine the concentrations elements in the sample. XRF has the benefit of being non-destructive and a range of instruments are available, from hand-held field instruments to more precise and sensitive laboratory-based systems; however, XRF may be limited by penetration depth of the incident X-rays, lack of ability to quantify elements with low atomic weights, and challenges in building calibration curves for quantifying elemental concentrations in complex, natural materials. Other advanced, non-destructive analytical techniques for measuring total concentrations of elements in solids, such as neutron activation analysis [100,101], are effective but often have limited availability for routine analysis of environmental samples. In contrast, total concentrations may also be assessed following total dissolution of the solid-phase material, often utilizing concentrated acids, heat and high pressures [101][102][103]; once dissolved, resulting aqueous samples can be analyzed by ICP-OES/MS and concentrations of elements of interest in aquifer solids may be back-calculated. Such destructive dissolution techniques have advantages of high throughputs but may be limited by the risks associated with strong acids (particularly HF) and the frequent inability to completely dissolve the solid-phase sample [104].
Approaches to understanding element-specific metal(loid) speciation in geologic materials are limited, with sequential extractions, spectroscopic methods, or microscale chemical imaging typically employed. Sequential extractions typically involve removing a target analyte (such as As, Cr, or V) using a series of extractant solutions that target increasing recalcitrant pools, often including cation-exchange extractants, carbonate-dissolving extractant, and strong acid extractants that target poorly order and crystalline minerals [105]. However, it is well established that these extracts do not exclusively solubilize materials from their target pool [106,107]. Nonetheless, this approach allows for quantitative assessment of how much of a specific element is associated with the operationally defined extractions, providing an insight on the relative availability [108]. In contrast, X-ray spectroscopic methods (commonly conducted at synchrotrons) utilize a monochromatic X-ray, probing all atoms a single element in a sample. The resulting spectra consist of two regions: an 'edge' region, where the adsorption of the X-rays increase rapidly and is characteristic of the element of analysis, and an extended region, characterized by longer oscillations following the edge. Analysis conducted in the edge region is known as X-ray adsorption near-edge structure spectroscopy (XANES), and is commonly used to determine the valence, and sometime fingerprint the local coordination environment, of an element. In the extended region, extended X-ray absorption fine structure spectroscopy (EXAFS), which is suitable for more concreted elements, is used to determine local coordination environment around an element [109]. Finally, mapping elemental associations and local speciation of elements in solid phases can be accomplished by a wide array of laboratory instrumentation; commonly used approaches include electron microscopy with energy dispersive spectroscopy (SEM-EDS), electron probe micro analysis (EPMA), micro X-ray fluorescence spectroscopy (μ-XRF) and secondary ion mass spectrometry (SIMS). Each of these approaches allows for spatially resolving the distributions of metal(loid)s in relation to other elements (and by association, mineral phases) of interest; however, synchrotron-based μ-XRF elemental mapping has the added benefit that specific metal(loid) oxidation states may be discerned when images are generated across a range of incident X-ray energies [110].

Mechanistic Modeling
Mechanistic modeling is the gold standard for describing environmental systems because it integrates understanding of fundamental processes into a mathematical representation of how the subject system behaves [111][112][113][114][115][116][117]. Such approaches have the promise of improving fundamental understanding of system behavior, predicting unobserved spatial or temporal states, and supporting design, management, and policy approaches to meet a set of constrained objectives. However, formation of incisive models, particularly those rooted in first principles, is complicated even for relatively simple systems. The development of such models in the subsurface is at the forefront of the state of science. In this section, we discuss the steps to formulate mechanistic models seeking to build on conceptual models that incorporate hydrogeochemical heterogeneity, mineralogical diversity, and biogeochemical complexity.

Fundamental Concepts
Subsurface systems have been modeled using mechanistic approaches both as individual systems and as components of coupled systems [111][112][113][114][117][118][119][120][121][122][123]. The well-established nature of mechanistic modeling of subsurface systems suggests a mature state of knowledge that belies the challenges associated with the mechanistic description of certain complex systems, such as those of focus herein. Because mechanistic models mathematically represent our fundamental understanding of a system, realistic models are underpinned by such a mature understanding [112][113][114]116,117,120]. When that mature understanding does not yet exist, limitations with the ability of mechanistic models to describe the system of concern naturally follows.
Even for cases in which fundamental understanding of the operative processes are well understood, challenges with mechanistic modeling can exist. For example, mechanistic models are systems of partial differential algebraic equations that necessarily require initial and boundary conditions, and model parameters (e.g., material properties, reaction rates). When these conditions vary in space and time in such a way that specification of the necessary model inputs are stochastic, uncertainty in the resultant model is a natural consequence [114,117,[123][124][125][126][127][128][129][130][131][132]. Such uncertainty in mechanistic models of subsurface systems is universally the case, but often ignored by purely deterministic modeling approaches.
The issue of scale also introduces an important challenge to mechanistic modeling of certain subsurface systems. Let's presume that at some small scale, which we will call the microscale ( Figure 5), the location and properties of the solid and fluid phase distributions are fully resolved in space and time, and all processes are described accurately using mechanistic approaches. Under such conditions, a mechanistic microscale model can provide an accurate description of the system [113,[117][118][119][120][121]. However, natural subsurface systems of concern have length scales and media complexity that preclude the use of microscale models [111][112][113]116,117,133,134]. To produce mechanistic representations that are tractable, a larger scale representation, called the macroscale, must be used. At the macroscale, the mechanistic representation is based upon an averaged behavior of the underlying microscale representation of the system [111][112][113]116,117,119,133,134]. Such a representation can be approached with varying degrees of rigor ranging from phenomenological to detailed first-principles mathematical approaches for deriving such models. Regardless of the rigor of the chosen approach, macroscale models require approximations of underlying microscale behavior that can introduce additional sources of error [114,116,117,133,[135][136][137]. This error, referred to as modeling error, poses a challenge to macroscale modeling even if perfect knowledge of auxiliary conditions and model parameters exists.

Modeling Metal(loid)s
With this brief introduction into the important fundamentals of mechanistic modeling of subsurface systems, let's consider the specific case of metal(loid)s of focus in this work. Observations of both laboratory and field-scale systems reveal the importance of geogenic sources and biogeochemical reactions as dominant factors influencing human exposure to these compounds in drinking water [83,[138][139][140][141][142]. Mechanistic models can bridge these disparate observational scales to produce insight into system behavior that can be used to protect human health. What follows is a discussion of the contributions and challenges associated with mechanistic modeling of metal(loid)s systems that considers a range of scales that successively increase in complexity. Each successive increase in scale contributes new complexity and challenges.

Batch Models
Batch systems, depicted in Figure 5a, are zero dimensional in space and one-dimensional in time, making them an ideal starting point for our discussion. At the batch scale, advanced spectroscopic techniques can be used to elucidate metal speciation, mineral surface characteristics, and binding mechanisms [114,[143][144][145][146]. Relatively refined descriptions of solution composition are also accessible for such systems. With such data in a wellmixed system, a mechanistic chemical species, reaction, and mass transfer model can be deduced based upon conservation of mass equations [114,141,[143][144][145][146]. Mature publicly available software exists for such purposes, such as PHREEQC, TOUGHREACT, and CrunchFlow, among others [134,147,148]. Parameters and a representative model related to mass transfer to the solid phase are needed, but the homogeneous water phase reactions and corresponding rate coefficients exist in available databases [147][148][149]. In short, at the batch scale, mechanistic modeling of the systems of concern is relatively straightforward. The precise form of such a model will depend upon the characteristics of the solid phase and the total set of species and reactions that are considered. Such approaches lead to a large number of species and component reactions that must be considered [83,118,124,[141][142][143][144]150]. In some cases, the solid particles can be further divided into a set of minerals and organic components, leading formally to a set of solid phases [114,118,141,145,146,151]. The goal of a mechanistic model is thus to describe the species-phase composition in the batch system as a function of time, and the initial conditions must be specified to render such a system solvable.
Batch models are made up of a set ordinary differential-algebraic equations with the number of equations dependent upon the number of unique non-equilibrium species-phase combinations [134,141,147]. Given the totality of subsurface and solid-phase biochemical complexity, the number of equations in such a model can number in the dozens [83,118,123,141,143,146,152]. Despite the system complexity, success in modeling the time evolution of the system and the final equilibrium state can be expected for most systems. Batch modeling is a typical starting point when modeling complicated biogeochemical systems. Agreement to batch system data is the first validation one would seek for a biogeochemical model, and such models can provide additional insight. For example, a systematic study can be used to reduce the complexity of a problem by eliminating relatively unimportant species and reactions. Such model reduction informs the model formulation and mechanistically based correlations for larger scale systems.

Microscale Models
Microscale models, schematically illustrated in Figure 5b, are an extension to batch models that do not require that the system is well-mixed. In addition to not being well mixed, in most microscale systems motion of the fluid phase relative to the solid phase is permissible, such as for the case of water flowing through a porous medium system. Let's assume the common case in which the solid phase is relatively immobile compared to the fluid phase [111,118,[153][154][155]. In this case, the motion of the fluid phase transports dissolved species through an advective process [118,148,[156][157][158]. A model of this system would simulate fluid flow using conservation of mass and momentum equations, along with a state equation relating the fluid density to other state variables (e.g., pressure, composition, and temperature) [111][112][113]155,[158][159][160]. Such a model is embodied in the celebrated Navier-Stokes equations, and readily available solvers exist to solve for a microscale velocity field for complicated fluid distributions typical of porous medium systems [154,155,160,161]. Once a velocity field is available, a set of species transport and reaction equations must be formulated and solved to model the fluid composition [118,124,141,157,162,163]. The set of species composition equations will mimic the form of the batch model and would include diffusive transport mechanisms. The key point here is that extension of the batch model to a porous medium system described at the microscale follows a well-established path that has relatively few open questions.
The firm theoretical footing of microscale models notwithstanding, challenges remain. First, the model form is challenging to solve because of the need to solve both the Navier-Stokes equations and a potentially large system of coupled, nonlinear transport and reaction equations [118,124,141,148,153,157,161,[164][165][166]. Coupling occurs through both homogeneous phase reactions and interphase mass transfer [118,124,141,148,153,162,163,166]. If the computational challenges associated with solving a microscale model can be overcome, a model of roughly equivalent maturity and fidelity with the underlying batch model is expected-extending the mechanistic representation to porous medium systems. Mechanistic microscale modeling can be used to study the evolution of complex geochemical systems in space, for example informing the evolution of a system along a flowpath [126,155,158,163,[167][168][169][170]. A limitation is that the computational burden for solving such a system will limit the domain to a size on the order of typically less than a 20 times the characteristic length of the system in each spatial dimension, and usually much less than this for readily available computational resources [148,154,155,157,171]. For an unconsolidated porous medium system, the characteristic length is often taken as the Sauter mean grain diameter [155,[172][173][174]. Thus, practically speaking, microscale models of complex geochemical systems are restricted to very small domains compared to the length scale of natural systems of concern [129,148,[153][154][155]157,162,167,171]. Because microscale models are based on wellaccepted fundamental understanding, they can provide high-fidelity descriptions of a system that can also be used to advance macroscale models needed to represent field-scale systems [113,119,153,155,157,[162][163][164]175].

Macroscale Models
The limitations on modeling at the microscale means that macroscale models are needed to simulate field-scale systems of concern, which are depicted in Figure 5c,d. The transition from microscale to macroscale models is an area of active research. Macroscale modeling of single-fluid-phase flow is well established and routinely performed in practice with abundant practical benefits well documented [111,113,148,159]. For the geochemical systems of focus in this work, a large system of transport and reaction equations must also be formulated and solved, mirroring the process of developing microscale models based upon batch systems. The increased length scale at which macroscale models are resolved means that field-scale simulations may be possible, though these come with their own computational challenges [157,159,165,176].
Formulation of an accurate macroscale representation of the complex nonlinear reactions that occur at the microscale is one of the biggest challenges facing field-scale mechanistic modeling [124,137,157]. Currently, no one macroscale approach has been found that can produce reliable simulations of complex nonlinear geochemical fate and transport in all porous medium systems of interest [124,136,137,156,157]. What has been well established is that macroscale models do not follow a similar form to their microscale counterparts, and assuming so can introduce large errors [135][136][137]157,177]. Mathematically, the issue arises when attempting to represent averages of nonlinear microscale quantities as a non-linear combination of macroscale averaged quantities that is of the same functional form as the microscale form [114,135,136,177]. Put another way, assuming that each region of a macroscale model can be modeled as a batch system with known inputs and outputs introduces significant errors. A parallel exists with the representation of turbulence, where non-linear velocity deviation effects have been studied for over a century without complete success [178][179][180]. Modeling errors introduced by macroscale models is a serious impediment to the development of useful macroscale simulators for complicated geochemical systems [124,137,157]. Additional obstacles exist when the geochemical system of concern includes both a region above and below the water table. Such twofluid flow problems come with their own model formulation issues [153,158,159,176], although recent advances suggest a new generation of macroscale models may resolve these deficiencies [133,[181][182][183].
Even if accurate macroscale model formulations are developed, further issues related to the establishment of accurate initial conditions for geogenic contaminant sources and complicated redox conditions would remain [126,152,169,171,184]. This would be exacerbated by long-standing issues of the identification of material parameters such as permeability and porosity that vary significantly in space but can only be observed coarsely [123,124,132,152,184]. These issues have been found to have a greater impact on species transport compared to simpler flow problems [118,124,126,129,130,152,155,171,177]. This inevitably leads to models of subsurface systems that are best represented stochastically, reflecting existing uncertainty [123,126,132,152,177,184].

Conceptual Illustration: Bimolecular Reactive Transport
To illustrate the modeling approaches discussed above, as well as their limitations, we will describe the case of bimolecular reactive transport, which has been extensively discussed in the literature [126,155,157,163,170,177]. The significance of this class of problem stems from the recognition that bi-and monomolecular reaction cases have been identified as the building blocks of many significant reactions of interest [170,185,186], including metal(loid)s reactions in the environment [141,187]. With this motivation, batch, microscale, and macroscale model formulations will be discussed to illustrate the general principles described in the previous sections. Modeling errors that accompany the jump from each scale to the next will also be discussed to illustrate the limitations of common modeling approaches and to highlight research opportunities.
A general bimolecular reaction may be described by the following chemical equation: where the concentration distribution of species A, B, and/or AB might be of interest in environmental systems-implying that the concentration of all species should be resolved by modeling at any scale. The first step taken to describe reactive fate and transport in a porous medium is reaction rate characterization using batch experiments. A batch system is fully mixed and may be modeled using a system of reactive conservation of mass equations given by where C iw is the concentration of species i in the fluid phase w, r iw is the reaction rate, subscripts indicate quantities that are not at the macroscale, and i is a species index that may be A, B, or AB for the example considered here. Equation (2) is not solvable without a closure relation for r iw , which will depend upon the specific form of the reaction rate law [157,164,167]. Here, a second-order reaction rate law will be used based upon Equation (1) and written for species AB, given by where k ABw is the rate constant, and r Aw = r Bw = −r ABw . For a batch system, the concentrations of the species of interest as a function of time is the desired solution. A well-posed batch model requires that the initial concentrations of all species, C Aw0 , C Bw0 , and C ABw0 be specified. Although the model involves three species, the batch nature of the system allows for simplifications and only a single differential equation needs to be solved for this simple case, either analytically or numerically, from which the concentration of all species can be retrieved. If k ABw is unknown, experimental observations can be used to solve an inverse problem to determine this rate parameter [157]. For many geochemical reactions, the necessary reaction rate parameters are known and tabulated in available databases [147][148][149].
The next level of complexity in the hierarchy of scales is the microscale. At the microscale, the concentration of dilute species is distributed in both space and time, with spatial contributions due to advective and diffusive fluxes. A common microscale model consists of a system of advective-diffusive-reactive equations (ADRE) [118,157,164,177], given by where D iw is the molecular diffusion, and v w is the velocity of the water flowing through the system. The reaction rate expression in a microscale model corresponds to the batch model for all microscale points in the domain. Additionally, the traditional ADRE relies on the assumption of Fickian diffusion at the microscale [126,157,159,161,164,177]. Although each of these assumptions could be checked, the microscale ADRE has been found to model effectively many single-fluid phase systems of interest [155,157], and would be expected to apply similarly to the fate and transport of metal(loid)s in the environment. To render the microscale model complete, v w must be known, which in general requires a solution of the Navier-Stokes equations, a continuity equation [119,152,157,161], and an equation of state for the fluid density [111][112][113]159,188]. Although the batch model simplified to a one equation solve, the situation is more complicated at the microscale because of the field nature of the quantities, and the need to solve for the fluid motion.
Macroscale modeling within a porous medium is ideally developed based on microscale conservation principles. Macroscale models can be formulated either by posing the model phenomenologically, such that they have a form that is similar to the corresponding microscale model, or they can be derived by direct averaging of microscale equations. Although the former approach is encountered often enough to be considered standard [114,134,137,189], the latter approach can result in models that are consistent across scales. Careful averaging from the microscale to the macroscale provides a verifiable link between scales that can allow one to check any assumptions or approximations that have been used to derive and close the model at the macroscale [113,119,157]. Although such an endeavor is technical, it has been accomplished for many systems of interest [113,115,119,122,183], including for systems experiencing reactive transport [130,190,191].
The use of phenomenological macroscale continuum models has led to well-documented modeling error [114,117,130,135,137,155,157,177,189], and some in the literature to question whether continuum models are viable for reactive transport [114,158]. To illustrate the sources of error, we will describe the commonly used macroscale ADRE, which has been phenomenologically posed, followed by a discussion of modeling error that results from this approach. The standard macroscale ADRE is given by where D w is the assumed Fickian form of hydrodynamic dispersion experienced in the averaging region, and superscripts indicate macroscale quantities. Equation (5) relies on the standard assumption that the averaging region for the macroscale model is well mixed [135,136,155,159] and the porosity is constant in space and time. A result of this assumption is that the reaction rate may be expressed for species AB in our case as r ABw = k ABw C Aw C Bw , (6) where k ABw would be the same as in Equation (3). The Fickian approximation is not applicable in general for natural systems [130,155,157,161,170] and open questions exist regarding the dispersive term for reactive systems even for the limiting case where the Fickian approximation is reasonable [130,155,157,170]. These approximations lead to modeling error, and their individual contributions to error are still an open question [126,130,155,157,177]. Additionally, Equation (5) cannot be used in cases where dissolution/precipitation leads to changes to the porous medium [118,129,164,175,192], which is possible for some systems [141,187].
The phenomenological macroscale ADRE essentially amounts to treating each macroscale averaging domain as a batch system with known inputs and outputs, which is not an accurate treatment in most field cases. A classic experiment that investigated the errors introduced by Equation (5) for bimolecular reactive transport is available in the literature [135]. It was found that the macroscale ADRE regularly predicted a mass production of the bimolecular product species that was 10-20% higher than what was observed. Although this level of error may be acceptable for laboratory experiments that are on the order of centimeters, this error would accumulate rapidly, becoming intractable in field-scale simulations that have length scales on the order of kilometers, which would be expected when modeling metal(loid)s fate and transport in the environment [137,148,155,170,187]. Model solutions that maintain the phenomenological framework of the ADRE have been proposed, but most are based on hysteretic models [124,136,157]. Such models are not ideal for heavy metal reactive transport modeling, where often the history of the system is not known, where the system is heterogeneous, and/or where inverse modeling is the explicit goal of the model [124,170,187]. The next generation of macroscale simulators requires a reactive transport model that only requires knowledge of the current state of the system, as well as accessible medium parameters. Vigorous work on these issues is underway in the geochemical modeling community [116,118,134,155,185,189].
Despite the prevalence of toxic metals such as As in the environment, field-scale mechanistic reactive transport models are rarely reported in the literature, and when they are reported they rarely incorporate all relevant physics. Typically, attempts to simulate As fate and transport in groundwater are based on statistical models [39,40,193,194], or they are phenomenologically based on posited macroscale extensions to pore-scale physics [193,[195][196][197]. Some of the most common methods that are used in the field today rely on coupling reaction modeling using PHREEQC with some other macroscale solver [193,[195][196][197]; however, these macroscale simulators often make use of the traditional ADRE, which does not apply to most reactive transport scenarios. Recent efforts have been made to more rigorously develop mechanistic field-scale models for As adsorption by averaging up from the pore-scale [198], but the resultant models cannot be considered mature.

Challenges
Barriers to the development of reliable mechanistic models include theoretical issues related to upscaling nonlinear reactive transport models from the microscale to the macroscale, managing the complexity associated with biogeochemical systems that may include dozens of species and hundreds of reactions, and the stochastic aspects of subsurface systems, including source conditions for Ar, Cr, and V that can be in a variety of mineral forms that vary markedly in space in a manner that cannot be identified deterministically. Despite the significant challenges, work toward mechanistic modeling should continue to demonstrate the maturity of our level of understanding and to make inroads toward reliable predictive modeling of these complex systems. Potential paths forward include considering wellcharacterized model systems, complexity analysis to determine and represent the leadingorder biogeochemical reaction effects as a starting point, and theoretical advancements to reduce modeling error for nonlinear transport and reaction problems in general.

Outlook and Conclusions
The release of geogenic contaminants to groundwater is controlled by the chemistry of the subsurface environment, which is in turn controlled by the hydrology and composition of the aquifer substrates. Interactions between these two result in a wide variety of potential recation involving the contaminant metal(oid)s, including: acid/base and redox transformation, solubilization, precipitation, sorption and desorption, among others. Here, we explored the common and contrasting release and retention mechanisms for As, Cr, and V. As is known to be more soluble and toxic in its reduced form, As(III), rather than its oxidized form, As(V); conversely, both Cr and V are more mobile and toxic in their oxidized forms, Cr(VI) and V(V), rather than their reduced forms, Cr(III) and V(III). The contrasting geochemistries illustrate how challenging it can be to account for all the factors that may contribute to their release, particularly in environments with multiple potential contaminants and groundwater conditions (i.e., pH, redox potential, etc.) that can alter over both distance and time. However, in many cases interactions with mineral surfaces, such as sorption, dissolution, and oxidation, have commonality among the key drivers of the processes (e.g., Fe and Mn minerals). Modeling these interactions is important for determine how the potential hazardous metal(oid)s mobilize and move in the subsurface environment, which requires connecting the hydrology, geology, and geochemistry at scales ranging from the microscale to the macroscale.
The understanding of the hydrogeochemistry, and thus efforts to model, the subsurface are hindered by the heterogeneous nature of the subsurface. It is difficult to specifically determine the concentrations and distributions of As, Cr, and V in the aquifer substrates, the exact flow of the groundwater, how the groundwater conditions evolve, and what biotic or abiotic interactions are occurring to mobilize or hinder contaminants. The complex nature of the subsurface in turn makes scaling mechanistic models from the microscale to the macroscale difficult. Constructing and validating such models is an ongoing area of research, with improvements being made in technique, approach, and computation power, as well as the expansion of basic kinetic and sorption data from related laboratory experiments.
Moving our understanding and modeling of subsurface metal(oid) contamination forward will likely require a cross-discipline effort in research. Being able to understand the complexity of the subsurface requires a variety of complementary expertise (hydrology, geochemistry, quantitative modeling), so designing projects that bring together members of different fields will facilitate the transfer of information. Drawing from experts across disciplines also reduces the potential for gaps in fundamental knowledge, and address the challenges to developing a comprehensive body of basic knowledge around the hydrogeochemistry and creating robust quantitative models. Furthermore, an important portion of this research will be making the information available to those outside of the scientific community who will have a need for it. This will include not only the general public who may need to be aware of potential issues with groundwater should they utilize it, but also stakeholders in adjacent fields. This may be as varied as public health officials, who would need to know what communities might be impacted by potentially hazardous contaminants, to land managers and environmental regulators, who would need to know where surficial inputs or spills might cause release of geogenic metal(oid)s in the groundwater. Above all, research into the solubilization and transport of metal(oid)s into groundwater, and potentially the drinking water of millions of people, is an urgent need, in order to remediate and prevent human exposure to potentially toxic, geogenic contaminants. Conceptual model showing geochemical evolution of groundwater in the saturated and unsaturated zones, including the hydrogeochemical factors that influence the contaminant chemistry at multiple scales. Dashed lines indicate lithologic boundaries (modified from [27][28][29]).  Stability diagrams for As, Cr, and V based on 10 μg/L for each element, HCO 3 − and SO 4 2− activities of 10 −3 and 10 −4 M using data from [6,69]. Figure constructed using the student version of the Geochemist's workbench. Redox controls on elemental mobility and toxicity that highlights contrasting redox controlled solubilities of As, V, and Cr. Low solubility species tend to form precipitates or sorb strongly to aquifer materials. Red indicates the oxidation states thought to pose more threat to health. OC = organic carbon. Mechanistic modeling approaches: (a) well-mixed batch system; (b) microscale modeling approach where the colored dots indicate reactive species; (c) a macroscale point, which can be resolved using microscale methods or averaged over at the macroscale; and (d) a field-scale system with a macroscale point denoted.