From legacy contamination to watershed systems science: a review of scientific insights and technologies developed through DOE-supported research in water and energy security

The Savannah River Site, located in south-central South Carolina, near Aiken, is an 800 square kilometer area where facilities were constructed in the early 1950s to produce special radioactive isotopes (e.g. Pu and tritium (³H)) for the Department of Defense nuclear weapons stockpile. The major facilities constructed included production reactors, chemical processing plants, and solid and liquid waste storage sites. It is estimated that the Savannah River Site has approximately $172 \times 10^6$ m³ of groundwater, soil, and debris contaminated with metals, radionuclides, and organics (National Research Council 2000). Contamination of the environmental resources is the result of disposal practices conducted on-site. For example, the F-Area Seepage Basins (located in the north-central portion of Savannah River Site) consists of three unlined, earthen surface impoundments that received approximately 7.1 billion liters of acidic low-level radioactive (e.g. U, iodine (I)) waste solutions.

The Oak Ridge Reservation (formerly known as Clinton Engineer Works) was established in the early 1940s as part of the Manhattan Project on 239 square kilometers of land in eastern Tennessee, approximately 40 kilometers west of Knoxville, Tennessee. Its original missions, conducted at three large facilities, included U isotope enrichment, construction, and operation of the world's first continuously operating nuclear reactor to demonstrate Pu production and separation, and radiochemical research and development. Today, the Oak Ridge Reservation covers an area of approximately 134 square kilometers and includes three major DOE complexes: the Oak Ridge National Laboratory (ORNL), the East Tennessee Technology Park, and the Y-12 National Security Complex (Y-12 NSC). Over the past 70 years, historical activities have resulted in significant releases of contaminants into the environment on the Oak Ridge Reservation, contaminating soils, sediments, groundwater, surface water, and biota. Waste disposal areas, unlined infiltration pits, trenches, and spills and leaks have created extensive subsurface contamination. Additionally, from 1950 through 1963, lithium (Li) isotope separation processes at the Y-12 NSC resulted in the release of about 212 000 kg of Hg into the environment, of which about 108 000 kg were estimated to have been lost to the East Fork Poplar Creek (EFPC). Both historical and ongoing Hg releases continue to negatively impact downstream waters, including many kilometers of the river–reservoir system outside the Oak Ridge Reservation.

The Hanford Site, located in southeastern Washington State, was established by the federal government to conduct defense-related nuclear research, development, and weapons production activities. As a result of activities at the Pu production complex and processing facilities (chemical separations plants and solid–liquid disposal and waste storage sites), it is estimated that the 1517 square kilometer area has approximately $1400 \times 10^6$ m³ of groundwater, soil, and debris contaminated with metals, radionuclides, and organics (National Research Council 2000). Contamination of the sites' environmental resources occurred as a result of 1.3 trillion liters of wastewater from chemical processing operations being intentionally discharged into the ground through settling ponds and other subsurface drainage structures (DOE 1997). Additionally, DOE estimates that 67 of the 177 underground high-level radioactive waste storage tanks have leaked 3.8 million liters or more of the highly radioactive waste into the subsurface (DOE 1997). Most of Hanford's subsurface contamination is concentrated in two locations: the Columbia River Corridor (100 Area and 300 Area) and the Central Plateau (200 Area). The 544 square kilometer Columbia River Corridor contained the production reactors, several waste burial sites, and major research facilities (DOE 1996). The 194 square kilometer Central Plateau, near the middle of the Hanford Site, contained the chemical processing facilities for extracting U and Pu from irradiated reactor fuel, associated waste storage facilities (18 tank farms), and waste disposal facilities (i.e. surface settling basins and underground drainage cribs) (DOE 1996). It is estimated that more than 220 square kilometers of groundwater at the Hanford Site is contaminated above current standards, mostly from operations in the 100 and 200 Areas. Section 5 briefly discusses the modeling of U contamination in the Hanford 300 Area and cesium (Cs) migration at the Hanford 200 Area.

The Nevada National Security Site (NNSS), located in southern Nevada, was the primary site for USA's underground nuclear testing. It led to the deposition of substantial quantities of 43 different radionuclides into the environment (Smith et al 2003). While ³H is the most abundant anthropogenic radionuclide deposited in the NNSS subsurface from an activity standpoint, Pu is the most abundant anthropogenic element by mass. Between 1951 and 1992, 828 underground nuclear tests (and 100 aboveground tests) were performed, and approximately 2.8 metric tons of Pu remains in the NNSS subsurface (DOE, Nevada Operations Office 2000). The site is approximately $\sim145$ km northwest of Las Vegas, Nevada, and covers approximately 3500 square kilometers. The DOE-EM program has been developing a monitoring strategy that involves identifying contaminant boundaries, restricting access to contaminated groundwater, and implementing a long-term monitoring program.

The Riverton, Wyoming DOE legacy site is a riparian floodplain that was the location of a U and vanadium (V) ore processing mill that conducted operations between 1958 to 1963 (DOE-LM 1998). It is situated at ∼$1500$ m elevation along the gravel bed of Little Wind River and has a semiarid to arid climate (25 cm mean annual precipitation). Temperatures exceed 0 ^∘C 154 days per year. The floodplain is vegetated with steppe flora, including sagebrush, grasses, and willows. Roots are relatively abundant, and robust upward solute transport through evapotranspiration has resulted in extensive evaporite mineralization of unsaturated soils (Dam et al 2015). The upper 2.5 m of soil is predominantly loam with abundant clay lenses, deposited abruptly over the underlying gravelly–cobbly–sandy alluvial bed. Groundwater contains persistent U, molybdenum (Mo), and sulfate (${\mathrm{SO}_{4}}^{2-}$) contamination from an upgradient legacy U ore processing facility. Soils are loamy, but clay layers and sand lenses are present above and within the unconfined shallow sandy–gravelly aquifer (2.5 m–3 m below the ground surface). The site experiences intense seasonal redox activity triggered by rising and ebbing water tables and frequent flood events. Clay lenses exhibit molecular oxygen (O₂) depletion upon water saturation and the establishment of reducing conditions (including iron sulfide (FeS) precipitation). Significantly, reducing conditions extend into the surrounding sandy aquifer and persist in proximity to clay lenses.

The Rifle, Colorado DOE legacy site is a riparian floodplain that was the location of a former U and V ore processing facility that operated from 1924 through 1958 (DOE-LM 1998). The Rifle Site is approximately 750 m in length along the Colorado River shore and 250 m at the widest point. The former processing facility contained large piles of mill tailings on-site, from which residual U and other associated contaminants leached into the subsurface. All surface structures and contaminated soil were subsequently removed, and the site was capped with fill material and vegetated. However, the potential for infiltration of groundwater contaminants remained until that time (Williams et al 2011, Long et al 2012). Site-specific investigations revealed that the groundwater contained dissolved U at the micromolar concentration level. Other contaminants of concern were identified as As, selenium (Se), and V. Groundwater at the site flows through alluvial flood plain deposits, which are 6 m–7 m thick and underlain by the relatively impermeable Wasatch formation. Groundwater enters the alluvial aquifer from upgradient sources above the floodplain and exits into the Colorado River.

Given this background, earlier research was focused on in situ immobilization of contaminants. More recently, an increased focus was placed on constraining the exchanges and fates of different forms of C and N in river–floodplain settings because of their important roles in driving biogeochemical interactions with contaminants. More broadly, efforts were dedicated to developing a fundamental understanding of microbial communities and how they mediate biogeochemical cycles in the terrestrial subsurface. More than 15 years of focused investigations at the site have created a rich legacy of understanding subsurface contaminant transport and provided transferable insights.

The East River Watershed, located in western Colorado, is a mountainous headwater testbed for exploring how climatic perturbations impact downgradient water availability and quality. The East River Watershed is located approximately 160 kilometers southeast of the Rifle Site and covers approximately 300 square kilometers. The East River is a pristine watershed. The Berkeley Lab's Watershed Function Scientific Focus Area (SFA; https://watershed.lbl.gov/) is examining the contribution to the riverine budgets of C and N that were neglected in previous assessments of contaminated sites like Rifle. The East River testbed began operation in 2014. Since then, it has enhanced process understanding of watershed dynamics and developed new approaches (such as monitoring strategies using sensor networks and scale-adaptive approaches) to improve watershed-dynamics simulation (Hubbard et al 2018, Arora et al 2020).

It is widely understood that microbial communities, not individual microbes, are the drivers of most subsurface biogeochemical cycles (figure 3(A)). Yet, in the early period of biological research at the Rifle Site, the focus was primarily on one family of bacteria, the Geobacteraceae, which have a laboratory-demonstrated capacity for U reduction (Gorby and Lovley 1992). This body of work motivated in situ acetate injection experiments to test the potential of a Geobacter-based bioremediation strategy at the Rifle Site (Anderson et al 2003). Proteomics of subsurface-derived cells (metaproteomics) were attempted to assay for Geobacter-specific proteins in amended aquifer groundwater. The only reference protein sequences available for this analysis were from Geobacter isolate genomes (Wilkins et al 2009). Thus, the effectiveness of protein identification was limited because the aquifer Geobacter genotypes differed from those in pure cultures.

The first proteomics experiment motivated one of the first cultivation-independent genomics (metagenomic) analyses of groundwater to improve the representation of Geobacter protein sequences. DNA from a Geobacter isolate not present at the Rifle Site provided an internal standard that was used to verify the accuracy of the recovered sequences. The experiment uncovered the very substantial microbial diversity of groundwater microbiomes. Draft genomes were reconstructed for over 80 bacteria, most of which were from previously unsampled or unknown lineages. Approximately half were from a group of novel bacteria whose gene content suggested they were symbionts of other microorganisms (Wrighton et al 2012). These groups were subsequently deeply sampled in a larger acetate injection experiment at the site and from unamended groundwater, and are now recognized as essentially ubiquitous in groundwaters worldwide (He et al 2021). The lineages share a common ancestry, so they were described as the Candidate Phyla Radiation (CPR) and were shown to comprise a substantial fraction (15% to ∼$50$%) of all bacterial diversity (Brown et al 2015, Hug et al 2016). The uniformly small genomes for CPR bacteria (Wrighton et al 2012, Kantor et al 2013) predicted their small cell size and thus potential to pass through the filters routinely used for cell recovery. Imaging of cells from the smallest filtrate fraction from acetate-amended groundwater revealed that the cells are around the theoretical minimum size for life (Luef et al 2015).

Also first detected in Rifle groundwater were novel archaea with genomic features reminiscent of CPR bacteria. These diverse archaea, now part of the DPANN group, substantially extended the breadth of known archaeal diversity (Castelle et al 2013). Given that both CPR bacteria and DPANN archaea are generally predicted to grow as surface-attached symbionts of other cells, it has been inferred that they have the potential to impact biogeochemical cycles through their impacts on host cells as well as via their own capacities for C, H, N, and S compound transformations. Subsequent research identified numerous bacteriophages (phages), the virus of bacteria, in subsurface communities (Al-Shayeb et al 2020), bringing them into focus as important contributors to C compound turnover.

Sampling of subsurface microbial communities from pumped pore fluids overlooks particle-associated cells. This motivated the study of sediment-associated consortia and uncovered very significant microbial diversity (Castelle et al 2013). One of the first applications of long-read sequencing to complex microbiomes revealed the presence of thousands of different organisms, most of which are at similar, very low abundance levels in Rifle sediments (Sharon et al 2015). Metagenomic analyses of saturated whole sediment samples recovered from 4, 5, and 6 m depths identified numerous new lineages, one of which is the bacterial phylum Zixibacteria. From a complete, curated genome, Zixibacteria were predicted to have an incredible diversity of metabolic capacities that enable growth despite shifts in conditions and the available aquifer resources (Castelle et al 2013). Another study targeting sediment-associated Chloroflexi revealed, among other traits, an unexpected pathway for carbon dioxide (CO₂) fixation that brings into focus the potential for this process in C cycling in the aquifer (Hug et al 2013).

A study that utilized groundwater and sediment metagenomes sampled from below and above the water table indicated that microbial cohorts consistently establish across the aquifer (Hug et al 2015). This, along with the desire to develop models for integrated microbial community function (Zhuang et al 2011), motivated a major genomes recovery effort. Genomes from over 2000 different microbial community members were analyzed simultaneously (an unprecedented effort). The major conclusion from this study was that biogeochemical cycles are attributed not to single organisms with capacities for specific pathways but to consortia whose metabolisms are interlinked by handoffs of pathway intermediates (Anantharaman et al 2016). This ecosystem structure probably confers resilience and has since been documented by genomics-based studies in other ecosystems, including soil (Diamond et al 2019).

Important limitations of the Rifle research site were its comparatively small size, lack of access to undisturbed soil and riparian zone sediments, and limited vegetation, slope, and aspect types. Thus, methods developed at the Rifle Site were scaled up to tackle the challenge of watershed ecosystem-scale microbial biogeochemical process analyses that incorporate genome-derived information. Initial studies focused on the soils surrounding the meanders of the East River Site in Colorado. Motivated by the need to link micron-scale microbial analyses to a transect that is many kilometers in length, it was hypothesized that soils within the arcs of successive meanders would represent a system representative of meanders along the river length.

Metagenomic datasets were generated from a sampling grid at three sites in the river's upper, middle, and lower reaches. Extensive genome recovery enabled documentation of biogeochemically important capacities across each site. A core floodplain microbiome was identified and found to be enriched in capacities for aerobic respiration, aerobic carbon monoxide oxidation, and thiosulfate oxidation. Metatranscriptomic data revealed that the most highly transcribed genes were amoCAB and nxrAB (for nitrification). Low soil organic C correlates with the high activity of genes involved in methanol, formate, sulfide, hydrogen, and ammonia oxidation, nitrite oxidoreduction, and nitrate and nitrite reduction. Overall, it was concluded that the meander-bound regions serve as scaling motifs that predict aggregate capacities for biogeochemical transformations in floodplain soils (Carnevali et al 2021). Similar research has addressed patterns of organism distribution and function across hillslopes and under four different vegetation types.

DOE's research in water and energy security arose from its historical and present-day need to resolve the many environmental management issues across its vast nuclear facilities, and its responsibility to develop a long-term solution to the USA's growing nuclear waste stockpile associated primarily with civilian nuclear energy production. The range of radiologic contaminants quite literally spans the entire periodic table: from ³H to the short-lived heavy elements (e.g. curium-244 (²⁴⁴Cm)) (Kurosaki et al 2014). More traditional contaminants of concern include Cr, Hg (discussed in section 6), and highly radioactive Cs and radioactive U (discussed in section 5.2), and more recently C and nutrients, including N and phosphorus (P). Below, we briefly describe DOE advancements in understanding the biogeochemical processes impacting contaminant fate and transport, as well as the major elements and nutrients (i.e. C, N, Fe, S, and manganese (Mn)) that drive subsurface biogeochemistry.

4.2.1. Contaminant biogeochemistry in redox-dynamic environments

4.2.2. Biogeochemical critical elements that mediate contaminant behavior

4.2.3. C and N biogeochemistry

4.2.4. Coupling of contaminant and nutrient redox cycles

Hydro-biogeochemistry is essential for tracking and monitoring contaminants and nutrients that affect river water quality through subsurface geochemical fluxes. The subsurface geochemical behavior of contaminants and nutrients differs across sites due to distinct environmental, geological, and hydroclimatic conditions (Zachara et al 2013, Neeway et al 2019). In particular, sediments, mineral composition, geology, and hydroclimatic conditions give rise to contrasting mechanisms controlling plume persistence in the subsurface and riparian corridors. For example, several studies have demonstrated that the effect of mineral surfaces and sediments is important for the transformation and export of geochemical species (Keppler et al 2000, Fox et al 2012, Allard and Gallard 2013, Maher et al 2013, Cumberland et al 2016, Dwivedi et al 2017a, 2019), while hydrology plays a critical role in mobilizing contaminants (Li et al 2010, Flores Orozco et al 2011, Williams et al 2013, Dangelmayr et al 2017, Jemison et al 2020).

In a comparative study, Zachara et al (2013) outlined key differences in the geochemical nature of residual contaminant U at the Rifle and Hanford 300 Area sites. They suggested that the major differences between the sites and the kinetic processes influencing the U distribution in space primarily resulted from the spatial geochemical heterogeneity of detrital organic matter and the magnitude of groundwater hydrologic dynamics controlled by river stage fluctuations.

At DOE sites, including the Rifle and Columbia River Corridor at Hanford, it was found that the interplay between groundwater dynamics and river stage oscillations may result in gradient reversal, changes in flow velocities, redox zonation, and alterations in groundwater composition (Greskowiak et al 2010, Yabusaki et al 2011, Lezama-Pacheco et al 2015, Dwivedi et al 2018a, Noël et al 2019). The relative variations in groundwater levels and river stages can alter the mixing zone's extent, where groundwater meets the surface waters and exchanges solutes (Ellis et al 2007, Fritz and Arntzen 2007). The mixing zone may span over distances of the meter- to kilometer scales (Dwivedi et al 2017b, 2018b, Zachara et al 2020). For example, Arora et al (2016b), Dwivedi et al (2018a), and Yabusaki et al (2017) demonstrated that the subsurface geochemical exports, and thus river water quality, were influenced by the spatial geochemical heterogeneity of detrital organic and temporal variability in groundwater flow direction at the Rifle Site. Similarly, the Columbia River Corridor site at Hanford also experiences frequent hydrologic gradient reversals, primarily due to dam operations and distinct subsurface hydrogeology control (Shuai et al 2019). Overall, the gradient reversal at sites like Rifle and Hanford can result in episodic oxic and anoxic conditions in the subsurface (Yabusaki et al 2011, Lezama-Pacheco et al 2015, Dwivedi et al 2018a, Chen et al 2021). Consequently, changes in the oxic environment due to groundwater mixing and microbially mediated sequential biogeochemical processes can lead to distinct redox gradients through a succession of electron acceptors.

4.3.1. Influence of hyporheic exchange flows and water-table fluctuations on mobilizing contaminants

4.3.2. Influence of hydro-climatic conditions on contaminant and nutrient dynamics

Over the last several decades, there have been significant advances in the characterization of above- and belowground properties and dynamics necessary to parameterize, calibrate, and validate reactive transport models. At the point and single-borehole scales, there have been rapid advances in in situ sensors to measure groundwater and surface qualities, such as soil moisture, pH, water table, and electrical conductivity. These sensors are typically placed inside wells to continuously monitor groundwater or surface water.

Geophysical methods—including electrical resistivity, seismic, and radar—have been increasingly used to characterize the subsurface in a noninvasive manner (e.g. Hubbard and Rubin 2005, Binley et al 2015). Geophysics can bridge the gap in sparse wellbore locations by providing high-resolution and spatially extensive information in a minimally invasive manner. It can image subsurface contaminant plumes (e.g. Johnson et al 2010, 2012, Dafflon et al 2011), as well as map flow and biogeochemical properties that are important for predictive modeling and understanding (e.g. Johnson et al 2010, 2012, Dafflon et al 2011, Wainwright et al 2014). In particular, geophysics has been powerful in identifying biogeochemical hot spots (Wainwright et al 2016a). For the watershed or regional scale, airborne geophysics—particularly airborne electromagnetic surveys—was originally developed for mineral exploration, but is now increasingly used for water-resources applications (e.g. Barfod et al 2018, Ball et al 2020). While a number of different field experiments have demonstrated the feasibility of such subsurface imaging techniques (e.g. Johnson et al 2015), real-time monitoring was still challenging until recently. In particular, autonomous electrical resistivity and phase tomography (ERT) monitoring have the potential to achieve rapid and automated detection and identification of changes in the subsurface (Dafflon et al 2011, Johnson et al 2015).

In parallel, recent advances in remote sensing have revolutionized how we characterize and monitor ecosystem and watershed dynamics. A high-resolution digital elevation model (DEM) derived from the Light Detection and Ranging (LiDAR) method is increasingly available across the USA, which has been applied to understand better the relationship between geomorphology and hydrology (Prancevic and Kirchner 2019), as well as to measure snow depths and snow-water-equivalent (SWE) over a basin scale (Painter et al 2016b). Because LiDAR data has been shown to inform near-surface soil properties (Gillin et al 2015, Patton et al 2018), hydrological connectivity (Jencso et al 2009) and biogeochemical hot spots (Duncan et al 2013) in the literature, it is being increasingly used in several DOE sites, including the East River (Hubbard et al 2018) and Columbia River Corridor Site at Hanford (Chen et al 2021). Similarly, hyperspectral remote sensing is used extensively in DOE- and non-DOE supported research to map plant traits (e.g. Asner et al 2015), leaf water contents (e.g. Colombo et al 2008), leaf chemistry (e.g. Feilhauer et al 2015) and other properties, which are also proxies for soil biogeochemistry (e.g. Madritch et al 2014). The East River is making use of these methods to characterize a range of properties, including leaf area, wet, and dry weights for leaf sample, strait models developed independently for needle and non-needle leaf species (Chadwick et al 2020a, 2020b, 2020c). For capturing the spatiotemporal heterogeneity, time-lapse images from satellites have enabled the monitoring of vegetation dynamics associated with hydrological disturbances (e.g. Wainwright et al 2020). The remote sensing images are available at increasingly high resolution and high frequencies (Falco et al 2018, Devadoss et al 2020).

Representativeness and quality of collected environmental datasets (e.g. groundwater, precipitation) impact accuracy and precision of climate, hydrological, and biogeochemical analyses and predictions. Often, these are incomplete and may contain unspecified or missing entries due to unavoidable reasons, such as unfavorable weather conditions or delays in collecting sensor data. Several statistical and machine learning-based methods have been developed to render Quality Assurance (QA)/Quality Control (QC) and address related issues using datasets at the East River (Dafflon et al 2020, Mital et al 2020, Faybishenko et al 2021, Dwivedi et al 2022).

The need to understand and predict contaminant transport and the export of nutrients and other important constituents has driven the development of much-improved modeling tools for multiphase (variably saturated) hydrological flow, such as TOUGH (Pruess 1987, 2004, Finsterle et al 2014, Jung et al 2018, Sonnenthal et al 2021), STOMP (White and Oostrom 2000, 2003), PFLOTRAN (Hammond et al 2007, 2014, Lichtner et al 2015), and reactive transport such as Crunch (Steefel et al 2005, 2015), ToughReact (Sonnenthal et al 2021), PFLOTRAN (Hammond et al 2007, 2014, Lichtner et al 2015). In turn, the availability of the new and evolving software has made it possible to address many scientific and engineering issues quantitatively at levels that were not possible previously.

5.1.1. Multiphase hydrology software

The multiphase (variably saturated) hydrological flow software developed for the DOE sites has historically been intended primarily for simulation of subsurface processes based on continuum formulations for porous media Darcy flow. Several hydrological codes such as HP1/HPx (Jacques and Šimŭnek 2005, Jacques et al 2006, 2012, 2013), PHT3D (Appelo and Rolle 2010, Greskowiak et al 2010), OpenGeoSys (Kolditz et al 2012, Nagel et al 2013), MIN3P (Mayer et al 2002), and HYDRUS (Šimŭnek et al 2006, 2008, Šimŭnek and van Genuchten 2008) have been developed globally since the late 1980s and early 1990s. However, many of the multiphase codes for variably saturated flow were developed within the DOE Laboratory system, including the TOUGH family of codes (Finsterle et al 2014), PFLOTRAN (Hammond et al 2012), STOMP (White et al 2012), NUFT (Hao et al 2012), FEHM (Zyvoloski et al 1997a, 1997b, Dempsey et al 2012), and Amanzi–Advanced Terrestrial Simulator (ATS) (Garimella et al 2014, Coon et al 2016, 2020). Variably saturated flow has been represented by the Richards equation, which assumes a passive gas phase, as well as with a full multiphase treatment in which multiple phases are subject to gradients and allowed to flow. Enhancements over the years have focused on improving the numerical stability of the Richards or multiphase flow simulations, which are notoriously difficult because of the nonlinearity of the problems they address. In addition, implementation of massively parallel computer architectures has been the focus of much of the effort (e.g. Hammond et al 2012, Garimella et al 2014, Jung et al 2017), although the details of this effort are beyond the scope of this review.

In order to extend beyond the subsurface so as to consider integrated surface and subsurface water, as required for analyzing terrestrial systems that consist of watersheds and river basins, much additional effort has been required. DOE-supported researchers have made important advancements to the state of the art in these integrated surface–subsurface hydrological models, with a focus on highly parallel open-source community codes. For example, extensions (Wu et al 2021) of PFLOTRAN (Hammond et al 2014) to include surface water and recent updates (Kuffour et al 2020) to the integrated surface–subsurface hydrology model PARFLOW (Kollet and Maxwell 2008, Kollet et al 2010) were supported by DOE projects. The Amanzi–ATS code (Coon et al 2016, 2020) was initiated as subsurface flow and reactive transport code (Amanzi) to support work at DOE legacy sites, and then extended as Amanzi–ATS to include surface flow and reactive transport processes to support research across DOE's network of watershed research testbeds.

The state of the art with respect to reactive transport modeling was summarized recently in Steefel et al (2015). In this summary of a special issue on benchmarking of reactive transport software for environmental applications emphasizing continuum formulations, such codes as PHREEQC, TOUGHReact, OpenGeoSys, PFLOTRAN, MIN3P, eStomp, HydroGeochem, and CrunchFlow were described. The first of these benchmarking workshops was supported by DOE-BER, although at a bare-bones level. The codes were shown to be capable of achieving similar or identical results for a range of complex environmental problems that required numerical solutions, i.e. those with nonlinear or coupled processes that are well beyond what can be described with analytical solutions subject to many simplifying assumptions. A significant portion (but not all) of this software development was undertaken with funding from the DOE, while the applications of the software were motivated in many cases by environmental problems on DOE lands and later by the desire to understand and predict biogeochemical cycling in ecosystems.

The development of biogeochemical reaction capabilities that could be used in reactive transport codes has followed a somewhat separate path. These efforts have followed closely the fundamental scientific work on understanding the role of microbes in surface and subsurface environments, particularly as they contribute to contaminant retardation through their control of the redox environment (NABIR 2006). Thus, the focus in reactive transport development here has been on incorporating arbitrarily complex biogeochemical reaction networks mediated by one or more competing or collaborating microbial communities. The microbial communities mediate rates that depend on the local environmental conditions and thermodynamics and contribute to reaction pathways that affect contaminant sorption (and thus retardation) and mineral dissolution + precipitation, as well as gas exchange. One example of such a biogeochemical reaction network with multiple microbially mediated pathways is shown in figure 4 (Arora et al 2016b).

Zoom In Zoom Out Reset image size

As might be expected for complex environmental problems like those facing the DOE Hanford Site, where several processes come into play over a large spatial extent, massively parallel computing for reactive transport is now playing an important role. Arguably the first parallel code for reactive transport was the MCTracker code described and verified in Yabusaki et al (1998). Following this by a few years was the parallel software PFLOTRAN (Hammond et al 2012) that built on many of the capabilities found in FLOTRAN. The power of this new code and its application to the real-world environmental problem of U contamination in the Hanford 300 Area was highlighted in Hammond and Lichtner (2010), discussed below (figure 5). As part of the DOE Advanced Simulation Capability for EM (ASCEM) program, the high-performance computing code Amanzi was developed to simulate subsurface flow and transport. This Amanzi framework was later integrated into the ATS code to bring in surface processes like freeze and thaw (Painter et al 2016a), overland flow, and river flow (Coon et al 2020).

Zoom In Zoom Out Reset image size

5.1.2. Integrated surface–subsurface reactive transport software

Modern reactive transport methods have made significant contributions to the topic of contaminant transport in recent years, and many of these applications are to DOE legacy sites or to sites where DOE contributed funding (e.g. Naturita, Colorado). The contributions can be broadly divided into two important themes:

• (a)  
  The comprehensive and rigorous treatment of geochemistry in contaminant transport models, and
• (b)  
  The role of aquifer heterogeneity in transport and mixing rates.

The second of these topics is not considered further here, since it is primarily the purview of hydrology and transport disciplines and is beyond the scope of this review paper. We focus here on the geochemical and mineralogical aspects of the problem.

The starting point for contaminant hydrogeology has always been the linear distribution coefficient (or $\mathrm{K}_{\mathrm{d}}$) models for sorption, which have the advantage of being easily incorporated into transport equations without rendering them nonlinear. It has also been asserted that they have an advantage in that the data required for the implementation is minimal, or at least less than that required by the multicomponent models discussed below. However, it should be pointed out that the $\mathrm{K}_{\mathrm{d}}$ models, while potentially being as simple as the determination of the sorbed and aqueous concentration under a set of environmental conditions, typically require unique constraints for each and every site to which they are applied. In addition, as will be apparent from the discussion below, such $\mathrm{K}_{\mathrm{d}}$ models may not even apply to a single site where conditions (temperature, salinity, competing ion concentrations, sorption site density) change over time and space. In other words, there is no real generality in the case of linear distribution coefficients, in contrast to the more rigorous ion exchange and surface complexation models (like that of metals on Fe-hydroxides) that can be applied quite widely.

The use of linear $\mathrm{K}_{\mathrm{d}}$ coefficients came under strong attack when it became apparent that ¹³⁷Cs had migrated further than expected as a result of highly saline tank leaks in the Hanford 200 tank farm (Zachara et al 2002, Steefel et al 2003, 2005). To address this issue, multicomponent ion exchange and surface complexation models were then considered in the context of reactive transport. Below, we briefly review them.

5.2.1. Ion exchange and Cs migration at Hanford 200 Area

The problem of understanding and predicting ¹³⁷Cs migration was made more complicated as a result of the multiple exchange sites that occur in sediments below the Hanford 200 tanks (Zachara et al 2002, Steefel et al 2003). In the Hanford case, the effective $\mathrm{K}_{\mathrm{d}}$ depends on the Cs⁺ concentration in addition to the competing sodium ion (Na⁺) concentration (figure 6(A)), since the ion exchange sites with a very high preference for Cs⁺ are present in only low concentrations in the sediment. The low concentration of the high affinity exchange sites means that they could be completely filled, resulting in a reliance on relatively lower affinity sites for the ¹³⁷Cs sorption.

Zoom In Zoom Out Reset image size

An example 2D reactive transport simulation was presented in Steefel et al (2005) to demonstrate the effect of the competing sodium nitrate (NaNO₃) (a major component of the leaking tank fluids) concentration. This simulation was carried out using ion exchange selectivities and site concentrations determined from Hanford bulk sediment (Steefel et al 2003). Figure 6(B) shows the transport of the nonreactive nitrate, with the plume extending over the entire domain size considered. The fractional migration of ¹³⁷Cs at 1 M NaNO₃ is significantly less, in keeping with the expectation that some retardation of this contaminant will occur (although still more than would be the case for ¹³⁷Cs in a typically dilute soil or vadose zone water). At 5 M NaNO₃, the migration of the ¹³⁷Cs is significantly farther, even if still retarded relative to the nitrate. The enhanced ¹³⁷Cs observed below many of the Hanford 200 tanks can thus be explained largely by a classical ion exchange model that accounts for the competing NaNO₃ concentrations in the plume, and by the elevated concentrations of Cs⁺ that exceed the number of high affinity sites that are available for strong sorption. The model also suggests that as leaking of highly concentrated NaNO₃ tank fluids ceases, further migration of the ¹³⁷Cs is unlikely. This is because infiltration by more normal rainwater is not capable of displacing the ¹³⁷Cs presently sorbed on the Hanford sediment. In any case, the use of constant linear $\mathrm{K}_{\mathrm{d}}$ to describe the transport fails completely, so either the more rigorous ion exchange models need to be used, or at minimum, a model that accounts for environmental condition-dependent sorption (i.e. a 'smart $\mathrm{K}_{\mathrm{d}}$').

5.2.2. Surface complexation and U migration

Perhaps an even more convincing demonstration of the inadequacy of the classical $\mathrm{K}_{\mathrm{d}}$ models is provided by reactive transport analyses of metal and radionuclide migration influenced by surface complexation on Fe-hydroxides. In fact, an entire generation of geochemists has investigated the surface complexation behavior of metals on ferric hydroxides over many years (Dzombak and Morel 1990, Davis et al 1998, 2004), but only more recently has the use of surface complexation models been demonstrated convincingly in reactive transport frameworks at the field scale.

Most likely, the first successful demonstration of the use of surface complexation models to describe field-scale reactive transport was presented by Davis and co-workers for the DOE-LM U-contaminated site at Naturita, Colorado (Curtis et al 2004, 2006, Davis et al 2004). Davis et al (2004) used both electrostatic and nonelectrostatic surface complexation models to describe U sorption on the Naturita sediment, but finally settled on the nonelectrostatic models because of their flexibility in treating natural and complex multimineralic sediments (figure 7(A)). In order to match the pH dependence of sorption, in particular, it was necessary to use a three-site SCM consisting of weak, strong, and very strong sites. Figure 7(B) shows the match with the experimental data at a variety of CO₂ partial pressures, an important variable because of the strong competition between U carbonate complexes in solution and the surface complexes developed on the Naturita sediment surfaces (Hsi and Langmuir 1985, Prikryl et al 2001, Davis et al 2004). Calcium uranium carbonate complexes in solution can further reduce the sorption of U, especially Ca₂UO₂(CO₃)$_{3_{aq}}$ (Bernhard et al 2001, Brooks et al 2003, Fox et al 2006).

Zoom In Zoom Out Reset image size

Curtis et al (2006) used the nonelectrostatic surface complexation model presented in Davis et al (2004) to investigate field-scale transport at the Naturita, Colorado Site (figure 7(B)). They demonstrated that the surface complexation model accurately reproduced the available data when combined with the realistic flow and transport parameters for the aquifer. They also demonstrated the inadequacy of a constant $\mathrm{K}_{\mathrm{d}}$ model (figure 7(B)).

Many other studies have been carried out that demonstrate the ability of the surface complexation models to describe field-scale behavior, including those at the U-contaminated Rifle Site (Yabusaki et al 2007, Zachara et al 2013), and at the contaminated Savannah River Site (Arora et al 2018).

5.2.3. Modeling of microbially mediated reaction effects on U migration at the Rifle Site

The fundamental studies of biogeochemistry mediated by microbial communities undertaken as part of the DOE research portfolio over the last 10–15 years have contributed to a much improved ability to model redox processes in the subsurface. This is a key development for predicting the dependence of redox-sensitive contaminants like U and Tc in the subsurface. The redox state of the subsurface depends on the classical sequence of electron acceptor-donor microbially mediated pathways from aerobic respiration under oxidizing conditions to methanogenesis under reducing conditions. The redox state of the subsurface can be manipulated with the injection of a suitable electron donor that drives the system towards reducing conditions, as was demonstrated at the Rifle Integrated Field Research Challenge in western Colorado (Anderson et al 2003, Dullies et al 2010, Williams et al 2011). The engineered biostimulation of indigenous microorganisms like Geobacter (see section 4.1) to catalyze the conversion of U^VI to a reduced immobile form of U, U^IV, was demonstrated convincingly at the field scale, although the longevity of such treatments has remained in question.

In the field experiments undertaken at the Rifle Site, acetate was injected into a series of wells upgradient of a series of monitoring wells (figure 8(A)). However, to simulate the acetate injection and resulting U reduction and immobilization as a result of microbial (primarily Geobacter) activity, it proved useful to begin with 1D column experiments (Li et al 2009). With this approach, it was possible to show that the reactive transport modeling could capture both the sequence and distribution of redox processes (figure 8(B)).

Zoom In Zoom Out Reset image size

To capture the combined behavior of the acetate, U, and other electron acceptors/donors in the 2008 field experiment, it was necessary to use a fully 3D model for the site (Yabusaki et al 2011). The injection of acetate drives the system to reducing conditions, with the result that the U^VI is immobilized as solid U^IV, and the ability to capture this behavior (along with associated redox processes) at the field scale is demonstrated convincingly with the simulations as shown in figure 8(C) (Yabusaki et al 2011).

Similarly, to understand the biogeochemical behavior of U within the groundwater and aquifer solids at the Oak Ridge Field Research Center, and as part of their Integrated Field Research Challenge, a conceptual model for simulating an emulsified vegetable oil injection was developed and incorporated into the geochemistry code PHREEQC (Tang et al 2013a, 2013b). Numerous field experiments and subsequent interpretation of pyrosequencing and qPCR results in conjunction with geochemical data provided the basis for this hybrid microbial-geochemical conceptual model development (Gihring et al 2011). Consistent with the field observations, the model demonstrated the biologically mediated reductive immobilization of U in groundwater with various electron donors and under diverse geochemical conditions (Wu et al 2006a, 2006b, Watson et al 2013).

Long-term (∼$10$ year) observations made at baseflow along the length of EFPC demonstrate consistent opposing patterns of Hg and MeHg concentration. Total and dissolved Hg (Hg_T and Hg_D, respectively) concentrations decrease from the headwaters downstream. Conversely, total and dissolved MeHg (MeHg_T and MeHg_D) concentrations increase with distance downstream. Additionally, the MeHg concentrations show a seasonal dependence, with higher concentrations in spring and summer and lower concentrations in autumn and winter (figure 9).

Zoom In Zoom Out Reset image size

Discharge in EFPC increases with distance downstream, suggesting that one explanation for the decreasing Hg_T and Hg_D concentration is that the additional water is diluting it in the channel. Nevertheless, dilution or particle settling cannot account for the concentration changes. In other words, the fact that the flux of Hg increases with distance downstream suggests that legacy sources of Hg in the watershed contribute additional Hg mass to EFPC (Peterson et al 2018). These mass-balance-based estimates are supported by natural abundance Hg stable isotope signature studies (Donovan et al 2014, Demers et al 2018). By comparing Hg stable isotope signatures from different possible source areas and along the creek, hyporheic water was identified as the likely source of the added Hg, and the reduction of Hg^II to Hg⁰ by both microbially mediated and photochemical mechanisms was indicated.

Two controlling variables on Hg speciation and transport are NOM and dissolved organic matter (DOM), respectively, and total suspended solids (TSS). Much of the Hg transported downstream is associated with suspended particles comprising a mixture of mineral particles, particulate organic detritus, and algae, predominantly diatoms. The nature of the association between Hg and particle surfaces has broader implications with respect to the reactivity and biogeochemical cycling of Hg in the creek ecosystem, including its release from particles and availability for methylation. Synchrotron-based x-ray fluorescence mapping of the suspended particles indicated that the mineral-bound Hg was closely associated with Fe and NOM (Gu et al 2014). In conjunction with complementary Fourier-transform infrared spectroscopy (FTIR) analysis, these results suggested NOM-facilitated Hg bonding to Fe oxides, in which Hg was tightly bound to reduced S groups in the NOM, and the NOM was bound to Fe-oxides via its carboxyl and hydroxide functional groups. The diatom-associated Hg was found on the outer cell surface, suggesting passive sorption onto the cell surface rather than active uptake.

The filter-passing, or dissolved, Hg is strongly complexed by reduced S groups in DOM (Skyllberg et al 2006). Using ultrahigh resolution Fourier transform ion cyclotron resonance mass spectrometry (FTICR-MS) coupled with electrospray ionization, investigators found that molecules containing multiple N and S atoms dominated the complexation of Hg (Chen et al 2017). From both EFPC and Suwanee River, a 2:1 complex ([C₅HNS$_{2}-$Hg−S₂NHC₅]⁺) that dominated Hg complexation was identified in DOM isolates based on both molecular mass and Hg stable isotope distribution. Possible molecular structures of this complex were estimated using density functional theory calculations.

MeHg flux also increases with downstream distance, which is expected given that both concentration and discharge increase. As is the case for EFPC, MeHg is typically not the initial form of Hg entering the ecosystem; rather, it is formed in the environment from an inorganic Hg precursor via a microbially mediated reaction carried out by a diverse group of anaerobic bacteria possessing the hgcAB two-gene cluster (Gilmour et al 2013, Parks et al 2013, Podar et al 2015).

Numerous studies have demonstrated that watershed wetland abundance is positively correlated with MeHg concentrations in the water bodies that drain those watersheds (Wentz et al 2014). Nevertheless, the EFPC watershed has 3% wetland abundance, but EFPC has MeHg concentrations typical of watersheds with 30% wetland abundance. Like inorganic Hg, MeHg in the water column is either particle associated or complexed with NOM. An overarching goal of the research is to understand where and under what biogeochemical conditions Hg is methylated and MeHg demethylated.

Targeted sampling of multiple storm-driven flood events revealed specific linkages between the terrestrial and aquatic environments. For instance, precipitation events delivered DOM and associated Hg to the creek, but not MeHg or particulate phases. Additionally, the recovery of dissolved solute concentrations to pre-flood levels closely followed the flood hydrograph, except for MeHg_D which lagged significantly behind other solutes. Recovery of MeHg_D concentrations closely followed the re-establishment of the diel cycle of dissolved O₂ concentration (figure 10). Taken together, the results of targeted flood sampling implied that in-stream biological mechanisms are primarily responsible for MeHg generation in EFPC (Riscassi et al 2016).

Zoom In Zoom Out Reset image size

Natural cycles and anthropogenic activity affect stream ecology, water quality, and water quantity on timescales as short as months to hours. For example, the progression of the seasons leads to cyclic variation in creek discharge, water temperature, sunlight intensity, and duration. Well-established diel patterns of dissolved oxygen (DO) concentration and pH are driven by photosynthetic activity (Allan and Castillo 2007). Wastewater treatment plants and hydroelectric dam operations are anthropogenic causes of short-time-scale water quality variations and river discharge.

Do Hg and MeHg concentrations have analogous cyclic variations to things like the daily photocycle, pH, and DO, and if so, how do they compare with other water-quality variations? Multiple sampling campaigns have been conducted to monitor changes in water chemistry over 30 to 48 hours to address this question. Previous research has documented the photochemical demethylation of MeHg and the role that various photosensitizers, including NOM, play in enhancing or inhibiting the reaction (Qian et al 2014). Given these earlier results, lower MeHg concentrations were anticipated during the day, driven by photodemethylation. Instead, MeHg concentrations were positively correlated with the daily photocycle increasing during the day with overnight minima. While photodemethylation may have been occurring, the rate of MeHg production outpaced that of all demethylation reactions, resulting in net positive MeHg production (Brooks et al 2021).

Long-term observations demonstrated seasonally dependent watershed sources of MeHg to EFPC. Targeted flood event sampling pointed to in-stream as opposed to out-of-stream sources of MeHg, and that the recovery of MeHg concentrations parallels the recovery of diel patterns in DO. Diel MeHg concentration variations are closely correlated with the daily photocycle. Together, these observations led us to hypothesize that periphyton biofilms are a source of MeHg in EFPC.

Periphyton, or benthic algae, biofilms are complex assemblages comprising biotic (e.g. algae, fungi, bacteria, and archaea, insect larvae), abiotic (mineral grains, entrained sediment particles), and detrital (coarse and fine particulate organic matter) components often embedded in an extracellular gel-like matrix. These biofilms are in-stream, and their photosynthetic inhabitants produce O₂ in a pattern driven by the daily photocycle. During floods, the biofilms are buried and scoured, disrupting the daily DO cycle, which is gradually re-established as the biofilms return. However, bacterial MeHg production is only carried out by anaerobic bacteria. Specifically, the only known methylators are iron-reducing bacteria, ${\textrm{SO}_{4}}^{2-}$-reducing bacteria, fermentative microorganisms, and methanogenic archaea. For Hg methylation to occur in periphyton biofilms, these oxygenic systems must also harbor O₂-free zones where these anaerobic microorganisms can be active.

Chemical gradients within intact periphyton biofilms collected from EFPC were studied in the lab using microelectrode voltammetry. Using this technique, concentrations of DO, Fe^II, and S^−II were measured vertically from the periphyton-water interface down into the biofilm. Over a vertical distance of a few millimeters, conditions progressed from O₂-saturated to iron-reducing and ${\textrm{SO}_{4}}^{2-}$-reducing. The steep biogeochemical gradients within the biofilms were consistent with the biofilms harboring conditions conducive to MeHg production.

To further explore the role of periphyton biofilms in Hg cycling, biofilms were grown in the creek, returned to the lab, and amended with enriched stable isotopes of inorganic ²⁰¹Hg and Me²⁰²Hg (Olsen et al 2016, Schwartz et al 2019). The production of isotopically labeled Me²⁰¹Hg and loss of Me²⁰²Hg was monitored over time as measures of Hg methylation and MeHg demethylation, respectively (figure 11). Experiments were conducted in each of the four seasons. Although the biofilms were a net positive source of MeHg during most incubations, MeHg production was greater during spring and summer and for samples collected downstream versus upstream. Both results are consistent with the long-term field observations that show higher concentrations in those seasons and increased MeHg flux with downstream distance. Biofilms grown or incubated in the dark had substantially lower MeHg production, and in some cases, were net demethylating, which is consistent with the diel sampling in which MeHg concentrations increase during the day and decrease at night. Finally, biofilm structure was important to the measured activity and function of these communities; in biofilms with disrupted structure, methylation potentials decreased by 50%.

Zoom In Zoom Out Reset image size

High-throughput amplicon sequencing has been applied to methylating periphyton communities to understand better the complex interactions that occur within them (Carrell et al 2021). Bacterial and archaeal richness was lower in summer relative to biofilms in the autumn, while fungal richness showed the opposite pattern. Interestingly, Hg methylation potential correlated with numerous bacterial families that do not contain the Hg-methylation genes hgcAB, suggesting that overall microbiome structure and relationships within periphyton communities influence rates of Hg transformation.

Hg and MeHg in contact with periphyton are subject to other reactions in addition to methylation and demethylation, including sorption, desorption, and, in the case of Hg^II, reduction to Hg$^{\textrm{0}}$ and its reoxidation. These ancillary reactions were measured, and a novel model was created—the transient availability model (TAM)—to more accurately describe the periphyton methylation-demethylation results (Olsen et al 2018). Consistent with lab incubations data of field-derived samples from EFPC, the TAM reproduced concentration versus time behavior that had previously been thought to be anomalous. Importantly, the TAM predicts faster Hg methylation rates than previous models, with broad implications for net MeHg generation in lotic ecosystems.

Using Hg biogeochemistry as a use case, we demonstrated that long-term observations in the EFPC, a DOE site (see section 3.2, Oak Ridge Reservation Site), coupled with natural abundance stable isotope fractionation, identified sources and locations of Hg entering the creek. Additionally, it identified seasons of the year that supported higher MeHg concentration. Targeted event sampling demonstrated that discrete and short-duration events connected the terrestrial and aquatic ecosystems in the watershed. In other words, these floods constitute hot moments with major implications for material and energy exchange. Event sampling also eliminated the floodplain as a source of MeHg to the creek and focused attention on in-stream sources. High temporal resolution sampling over daily periods produced the unexpected result of increasing MeHg concentration during the day. Controlled laboratory experiments confirmed that periphyton biofilms are a net positive source of MeHg and constitute hot spots of Hg, Fe, and S cycling. Further, biofilms with a higher degree of connectivity among its populations also have higher Hg methylation potential. This combination of approaches buttressed some ideas, but more importantly, highlighted previous misconceptions, identified HSHMs of governing processes, and opened new avenues of discovery to advance further our understanding of these systems that support diverse and vital ecosystem services for the global economy.

Historically there has been a mismatch between the spatial and temporal scales at which fundamental hydrologic and biogeochemical processes are most readily studied and understood (molecular to laboratory) and those of the ecosystems that they impact (field to watershed/basin). Because of nonlinear process interactions and complex heterogeneity in natural ecosystems, small-scale (e.g. laboratory) studies do not directly translate to field settings. Key reactive transport parameters, including permeability, dispersivity, and reaction-rate coefficients are all well known to be scale dependent and, in some cases, boundary-condition dependent. These challenges have motivated a variety of modeling approaches that are aimed at integrating information across scales and representing ecosystem processes at multiple levels of fidelity and complexity through DOE- and non-DOE supported research.

Multiple approaches have been developed to address this challenge, each with advantages and limitations. A summary of selected approaches used in DOE-supported research is provided here; for broader reviews, see Keyes et al (2013) and Scheibe et al (2015a).

Multiresolution methods use variable discretizations to capture effects of small-scale processes and properties on large-scale phenomena (e.g. Özgen-Xian et al 2020). These methods include adaptive grid refinement (refining grids dynamically to put higher resolution where and when needed); multigrid methods (hierarchical sets of grids used to speed numerical convergence of a fine-grid solution), and nested or spatially varying grids that enhance resolution in specific subdomains (e.g. figure 12).

Zoom In Zoom Out Reset image size

Upscaling methods use assumptions about the nature of subgrid variability to derive effective parameters valid at the model resolution. These methods include numerical upscaling (direct solution of subgrid domains to derive approximate grid-scale parameters), various forms of homogenization and volume averaging, and perturbation solutions to stochastic partial differential equations.

Hybrid multiscale methods are similar to multiresolution methods. However, they use fundamentally different models in selected subdomains (for example, pore-scale simulations in portions of a domain where mixing-limited precipitation occurs; Scheibe et al 2015b) and are sometimes referred to as 'adaptive physics' models. These include hierarchical methods (e.g. the heterogeneous multiscale method; E et al 2003) in which overlapping domains exchange information and concurrent methods (e.g. mortar methods; Balhoff et al 2008), within which distinct subdomains pass information at domain boundaries.

The HSHM construct provides the basis for developing a robust predictive understanding of how contaminant behavior in subsurface heterogeneous media varies across spatiotemporal scales. Subsurface heterogeneities are sites of intense biogeochemical activity and are critical to watershed function. Recent studies have highlighted in greater detail the roles of the naturally reduced zones (or 'redox heterogeneities') on biogeochemical function. For example, heterogeneity caused by clay lenses in high permeability soils provides opportunities to sorb contaminants in place (Campbell et al 2012, Janot et al 2016), while fracture networks can lead to rapid infiltration of contaminants to groundwater (e.g. Berkowitz 2002, Arora et al 2019a). Thus, high contaminant concentrations can be spatially diffused across large scales or concentrated at higher rates in specific areas, thereby forming contaminant hot spots. Similarly, these concentrations can be temporally diffused or more concentrated during specific times, thereby exhibiting contaminant hot moments. In addition, hydrological perturbations in transiently saturated zones drive changes in soil microbial communities, triggering a range of geochemical responses, thereby forming dynamic HSHMs. Figure 12 illustrates static hot spots (i.e. fixed in space) and dynamic hot spots (i.e. migrating in space).

Although the HSHM construct simplifies the characterization of complex heterogeneous systems and time scales at which biogeochemical processes occur—by subdividing the region and time into a finite number of relatively homogeneous units—it lacks a quantitative definition. Moreover, the HSHM construct is underutilized because of its dynamic nature. In response to that, Bernhardt et al (2017) argued that the temporal dynamics of a putative hot spot are a fundamental trait that should be used in their description and, therefore, stressed the merger of HSHMs into a single descriptive term 'ecosystem control points.' This definition indicates that knowledge of the rate distributions has more relevance than knowledge of maximum rates. Partly supported by DOE, Bernhardt et al (2017) suggested four distinct categories of ecosystem control points: (a) permanent control points, experiencing high rates of biogeochemical activity relative to their proximity (e.g. floodplain), (b) activated control points, where higher biogeochemical activity is set off when certain substrates are delivered (e.g. low-lying topographic positions), (c) export control points that accumulate reactants until some hydraulic or diffusion threshold is not surpassed (e.g. geochemical exports through meanders), and (d) transport control points having extremely high transport capability without their direct intervention in exhibiting high biogeochemical activity (e.g. preferential flow paths). Table 2 presents some examples of HSHMs from DOE research in the backdrop of ecosystem control points.

7.2.1. Impact of HSHMs on biogeochemistry of legacy contaminants and nutrients

Modeling eco-hydro-biogeochemical processes requires the characterization of the bedrock-to-canopy properties across watersheds. The tremendous heterogeneity of these properties is often a significant obstacle toward scaling modeling efforts from the hillslope and floodplain scales to a larger scale. Although various spatially extensive characterization techniques such as remote sensing and geophysics are available, it is still challenging to estimate each property distribution individually over the watershed scale, particularly soil and subsurface properties in the subsurface.

A functional zonation approach holds tremendous potential for tractably estimating individual property distributions over the watershed scale. In this approach, zones or spatial units within a landscape with unique distributions of multiple properties relative to neighboring regions can be identified using machine-learning-based spatial clustering approaches (e.g. figure 12). The zonation approach essentially builds on a hydrofacies approach, which has been used in groundwater hydrology to define hydrological parameters in each facies (e.g. Fogg et al 1998, Klingbeil et al 1999). The reactive facies (Sassen et al 2012, Wainwright et al 2014) defined both hydrological and geochemical properties in each facies and stochastically mapped the reactive facies based on multiple geophysical data. Bea et al (2013) documented the value of incorporating estimated reactive facies information within a reactive transport model to improve predictions of long-term U plume transport. These concepts have proved valuable in identifying zones and hot spots of biogeochemical activity (Flores Orozco et al 2014, 2018, Wainwright et al 2015).

More recently, the availability of enhanced remote-sensing (including hyperspectral) products and advanced machine-learning approaches has led to a paradigm shift in the zonation approach and demonstrated possibilities for characterizing the bedrock-to-canopy properties across watersheds. With this new paradigm shift, the zonation approach harnesses a suite of remote sensing datasets for identifying and distributing zones relevant to biogeochemistry and/or ecosystem functioning. In addition, the zonation approach takes advantage of the current scientific understanding of critical zone processes and insights into water and nutrient cycling and their exports within a watershed and its different compartments, particularly through the Critical Zone Observatory (CZO) network (Brantley et al 2017).

Next, we briefly describe some prominent applications of the zonation approach at DOE sites. Hubbard et al (2013) and Wainwright et al (2015) used geophysical and remote sensing datasets to identify zones representing the heterogeneity of key properties and ecosystem functions—such as soil properties and C fluxes—in Arctic tundra. Devadoss et al (2020) defined the ecosystem functioning zones (each of which has distinct time series signatures of soil moisture and plant dynamics) and investigated the relationship among soil, plant, and snow dynamics. Wainwright et al (2016b) subsequently incorporated U reactive facies (e.g. Bea et al 2013) into the zonation approach to identify biogeochemical hot spots (e.g. naturally reduced sediments) based on electrical geophysical methods at the Rifle Site. Several additional modeling studies at the Rifle Site represented naturally reduced sediments in their models to investigate the HSHMs of N, C, Fe, and DO and quantified their subsurface geochemical exports (e.g. Arora et al 2016b, Yabusaki et al 2017, Dwivedi et al 2018a).

Recent investigations have been focused on examining how properties associated with specific functional zones exert controls on water and N exports in response to snow dynamics and contribute to the aggregated watershed concentration–discharge signature in mountainous regions like the East River (Hubbard et al 2018). A suite of remotely sensed bedrock-through-canopy watershed data layers and machine-learning approaches were utilized at the East River to identify watershed zonation. The zonation approach explored the sensitivity of different identified zones to foresummer drought (Wainwright et al 2022). Zonation is also powerful for understanding the co-variability of above- and belowground processes, as well as identifying relevant field experimental sites throughout the domain. For example, research at the Hanford Site is using river channel morphology zonation (hydromorphic units, Ren et al (2021)) as a framework for a transferable understanding of river–groundwater exchange fluxes and transit times.