Reactive transport modeling for supporting climate resilience at groundwater contamination sites

Climate resilience is an emerging issue at contaminated sites and hazardous waste sites, since projected climate shifts (e.g., increased/decreased precipitation) and extreme events (e.g., flooding, drought) could affect ongoing remediation or closure strategies. In this study, we develop a reactive transport model (Amanzi) for radionuclides (uranium, tritium, and others) and evaluate how different scenarios under climate change will influence the contaminant plume conditions and groundwater well concentrations. We demonstrate our approach using a two-dimensional (2D) reactive transport model for the Savannah River Site F-Area, including mineral reaction and sorption processes. Different recharge scenarios are considered by perturbing the infiltration rate from the base case as well as considering cap-failure and climate projection scenarios. We also evaluate the uranium and nitrate concentration ratios between scenarios and the base case to isolate the sorption effects with changing recharge rates. The modeling results indicate that the competing effects of dilution and remobilization significantly influence pH, thus changing the sorption of uranium. At the maximum concentration on the breakthrough curve, higher aqueous uranium concentration implies that sorption is reduced with lower pH due to remobilization. To better evaluate the climate change impacts in the future, we develop the workflow to include the downscaled CMIP5 (Coupled Model Intercomparison Project) climate projection data in the reactive transport model and evaluate how residual contamination evolves through 2100 under four climate Representative Concentration Pathway (RCP) scenarios. The integration of climate modeling data and hydrogeochemistry models enables us to quantify the climate change impacts, assess which impacts need to be planned for, and therefore assist climate resiliency efforts and help guide site management.

1 Introduction

Changing climate may pose a major risk in environmental remediation, especially with regard to fate and transport, including both hydrologic and reactive processes (Maco et al., 2018). In particular, many sites are managed with monitored natural attenuation strategies where an expanded contamination plume with high concentrations of tritium, uranium, and other chemical species remains in the subsurface (Denham et al., 2020). Hydrological shift has been identified as one of the key drivers influencing such risk and uncertainty. In a changing climate, precipitation and evapotranspiration (ET) regimes can change in both magnitude and timing, significantly affecting infiltration. Precipitation regimes are expected to change depending on where the site is located (e.g., Lambert et al., 2008). Increasing ET is usually predicted in climate model projection due to increasing temperatures under global warming (e.g., Abtew and Melesse, 2013; Milly and Dunne, 2016). Extreme events, such as heavy rain and prolonged droughts, are expected to become more frequent and thus may result in faster plume remobilization (e.g., Rahmstorf and Coumou, 2011).

We may define climate resilience at contaminated sites as the capacity of individual waste disposal sites to return to the system's original condition when affected by climate trends, climate variability, extreme events, and other climate-change-related impacts. A critical need exists for understanding climate change impacts on contaminated sites (e.g., U.S. Department of Energy, 2014; U.S. Environmental Protection Agency, 2017); however, a quantitative estimation with climate change projection is still missing. Evaluating the effect of climate change on the abundance of water resources has been widely studied (e.g., Gellens and Roulin, 1998; Green et al., 2011; Middelkoop et al., 2001; Pfister et al., 2004); however, water quality and contamination issues were less investigated (Visser et al., 2012). Most previous researches study surface water (Wilby et al., 2006; Van Vliet and Zwolsman, 2008; Van Bokhoven, 2006; Futter et al., 2009; Schiedek et al., 2007) because of the accessibility and data availability (Green et al., 2011). In the limited studies for climate change impacts on groundwater in the subsurface domain, agricultural effluents at the regional scale are the research focus (Bloomfield et al., 2006; Futter et al., 2009; Li and Merchant, 2013; Olesen et al., 2007; Sjøeng et al., 2009; Whitehead et al., 2009; Wilby et al., 2006; Darracq et al., 2005; Destouni and Darracq, 2009; Park et al., 2010).

Recently, Libera et al. (2019) investigated the potential impact of climate change on residual contaminants in vadose zones and groundwater, using a groundwater flow and transport model. They investigated the complex effect of precipitation and recharge shifts, leading to either dilution and remobilization of residual contaminants. Libera et al. (2019) showed the effects of dilution and remobilization on contaminant concentrations before and after changing precipitation, depending on the well locations, and that surface barrier and source-zone monitoring are critical for mitigating the impact. However, Libera et al. (2019) only simulated tritium with decay but did not couple with a reactive transport model to simulate other chemical species, sorption, and mineral reactions. In this study, we hypothesized that increasing recharge would decrease reactive species concentration further, since increasing the volume of water in the domain would increase pH, which limits the mobility of uranium. The impact of hydrological shifts on reactive contaminants is expected to be more complex, especially redox and pH-sensitive heavy metals. Remobilization would also be affected by additional clean infiltration water. To test those hypotheses and evaluate the impacts, process-based flow and reactive transport models that can characterize sorption and ion-exchange processes are essential for quantitatively analyzing the contaminant plume and understanding climate resilience.

This study aims to evaluate the effects of climate-driven hydrological shifts on (1) reactive contaminants and (2) mineral reactions in vadose zones and groundwater. We assume that the effect of changing precipitation and temperature can be represented by perturbations/shifts in natural recharge through the aquifer system. In our case, climate resilience is evaluated by the concentrations of monitoring wells after climate events in comparison to the background baseline concentrations. We demonstrate our approach at the Department of Energy (DOE)'s Savannah River Site (SRS) F-Area Seepage Basins, South Carolina (SC), USA, where soil and groundwater have been contaminated by various metals and radioactive contaminants. This is the same site studied by Libera et al. (2019). The SRS F-Area seepage basin was chosen because of the historic contamination monitoring activities and extensive characterization and model development (e.g., Flach, 2004; Bea et al., 2013; Sassen et al., 2012; Wainwright et al., 2014, 2019; Schmidt et al., 2018; Denham and Eddy-Dilek, 2017; Libera et al., 2019). More importantly, the contamination in F-Area is indicative that our study may provide broader insights into other contamination sites. The vadose zone residual contaminants are quite common (e.g., Stubbs et al., 2009; Zachara et al., 2005), and the contaminant export through the wetland region is a fairly common feature across many contaminated sites as well (e.g., Mansoor et al., 2006; Li et al., 2014; Chang et al., 2014). Hence, the SRS has become a unique study site for investigating the potential consequences of climate change for contamination remobilization and mineral reactions/interactions. Our research focuses on the effect of climate change on progress toward return to natural conditions of the plume between the basins and the funnel and gate. This is important for the timing of the transition of the site from enhanced to monitored natural attenuation and hence important for the overall effectiveness of remediation. More importantly, this work will support the risk management under changing climate conditions.

2 Site description

The Savannah River site F-Area in South Carolina is approximately 100 miles (i.e., 161 km) away from the Atlantic Ocean and occupies an area of about 800 km² (Fig. 1). The site was used to produce special radioactive isotopes, plutonium, and tritium for nuclear weapons during the Cold War era. The F-Area is located in the northern–central part of the SRS. There are three hydrostratigraphic units within the Upper Three Runs Aquifer, shown in Fig. 1b: an upper aquifer zone (UUTRA), a tan clay confining zone (TCCZ), and a lower aquifer zone (LUTRA). The UUTRA and LUTRA are mainly composed of clean sand, while the TCCZ is a low-permeability mixed sand-and-clay layer. The piezometric head measurements indicate that the UUTRA and LUTRA units are hydrologically connected. The bottom of the LUTRA consists of a competent clay layer confining unit that separates the deeper aquifer (Gordon Aquifer) from the upper two aquifer units (Fig. 1). The historical monitoring data collected at the SRS have shown that the F-Area contaminant plume migrates within the UUTRA and LUTRA (Fig. 2), discharging into a local stream called Fourmile Branch Creek (FMB).

Figure 1 (a) Location of seepage basins in the F-Area of the Savannah River site (SRS). (b) Hydrostratigraphic units defined for the F-Area; (c) 2D cross-section model domain. Modified from Bea et al. (2013).

Low-level radioactive acidic waste was disposed of in three separate unlined seepage basins (F-1, F-2, and F-3) and leached into the groundwater. The basins received approximately 7.1 billion liters of waste solutions from processing irradiated uranium between 1955 and 1988. After the waste discharge operation was terminated in 1988, the F-Area basins were closed and capped with a low-permeability material. Currently, an acidic contaminant plume extends from the basins approximately 600 m downgradient to the groundwater seepage near the FMB. Several measures have been taken to reduce the environmental impacts at the F-Area site, including capping the basins and pump-and-treat remediation of contaminated groundwater. Since 2004, the site has been undergoing enhanced natural attenuation using a funnel-and-gate system, which consists of groundwater flow barriers to decrease the groundwater gradient and base injection to neutralize pH and in turn immobilize uranium. The funnel-and-gate system is operating and requires the periodic injection of a base solution to increase pH and immobilize contaminants. Quantitative estimation from the modeling study will provide insights for site management and stakeholders on when it is appropriate to transition the site to natural attenuation status without any treatments. Despite the many active remediation measures, the groundwater remains unnaturally acidic upgradient of the funnel and gate and contaminated with various radionuclides.

One characteristic of the SRS F-Area is the high acidity of the plume, making U(VI) highly mobile. The natural groundwater pH is slightly acidic, between 5.0 and 6.0, and decreases to values approaching 3.2 in the most contaminated locations. Despite many years of active remediation, contaminated groundwater still remains highly acidic, and the concentrations of U(VI) and other radionuclides are still significant (Seaman et al., 2007; Savannah River Nuclear Solution, 2021). It should be noted that in the acidic pH range in the SRS F-Area, solid/liquid partition coefficients (Kd) values for U(VI) could change between 10² and 10⁶ (Davis et al., 2004; Dong et al., 2012). In addition, competing sorption between U(VI) and H⁺ is important in remediation and has been well studied in the F-Area site (Davis et al., 2004; Bea et al., 2013; Arora et al., 2018). Because of the difficulty and apparent uncertainty in assessing the adsorption properties and mobility of U(VI) under complex geochemical conditions in groundwater, several researchers have performed quantification (UQ) related to U(VI) and H⁺ competing sorption in the F-Area site (e.g., Curtis et al., 2006; Hammond et al., 2011).

3 Modeling methods

3.1 Reactive transport modeling with Amanzi and PFLOTRAN

Groundwater flow and contaminant transport are simulated by the numerical code Amanzi (Moulton et al., 2014, https://github.com/amanzi/amanzi, last access: 31 January 2022), which provides a flexible and extendable parallel flow and reactive transport simulation capability for environmental applications. Amanzi has the capabilities to solve coupled unsaturated and saturated groundwater flow as well as advection–dispersion transport equations. It includes a general polyhedral mesh infrastructure and provides multiple advanced nonlinear solvers with open source libraries. The reaction of contaminants and minerals carried by flow through the surrounding rock and soil is modeled by coupling with the geochemistry engine of PFLOTRAN (Lichtner et al., 2015) via the generic interface Alquimia (https://github.com/LBL-EESA/alquimia-dev). PFLOTRAN simulates the mineral reactions, adsorption, and ion exchange, while groundwater flow and transport are simulated by Amanzi.

3.2 Model setup and boundary conditions

The two-dimensional (2D) flow and transport Amanzi model developed in Libera et al. (2019) was employed and expanded with the coupling of a reactive transport model in this study. Our Amanzi simulation used the same conditions of mineral composition and kinetic reactions as the TOUGHREACT model in Bea et al. (2013). The 2D domain is approximately 2600 m long and 100 m deep along the groundwater flow line, passing through the middle of the F-3 basin of the SRS. Bea et al. (2013) calibrated the model and verified it using observational geochemical concentration data from several monitoring wells and also evaluated the sensitivities of key parameters in the modeling. The model includes the vadose zone and three hydrostratigraphic units (i.e., UUTRA, LUTRA and TCCZ) defined in the previous section. We assume homogeneous average hydrogeological properties within each unit (see Table 2), whose values are compiled from available site investigation reports. Table 1 specifies porosity, permeability, and capillary pressure/saturation data for the vadose zone (Flach et al., 2004; Bea et al., 2013). The system is considered to be advection dominated. Based on the study of scale-dependent advection and dispersion processes in Molz (2015) that states contaminant transport will be typically dominated by advection at a scale of 1000 m, the system is considered to be advection dominated, and mechanical dispersion and molecular diffusion transport processes are neglected.

Table 1 Physical model parameters used in the simulations. α, n, and Θr are the parameters of inverse air entry suction, a measure of the pore-size distribution, and residual water content, respectively, in the van Genuchten water retention curve.

Table 2 Chemical composition for the background (initial), recharge and seepage solutions (modified from Bea et al., 2013). Unit is mol kgw⁻¹, except pH and CO₂ (aq).

^a Calculated as electric charge balance. ^b Equilibrium with kaolinite. ^c Fixed by atmosphere pressure of CO₂.

No-flow boundary conditions are assigned along the two vertical sides of the 2D cross section (see Fig. 2) according to the groundwater divides, modified from previous studies (Flach, 2004; Bea et al., 2013). An impervious flow boundary condition (i.e., no-flow) is set at the bottom of the computational domain, since the confining unit at this location is highly clay-rich and continuous (Bea et al., 2013). Recharge rate is computed by the difference of climatological average precipitation and ET. This is appropriate for this domain and most groundwater models in which the groundwater domain is deep compared to the root zone depths.

Figure 2 Left: illustration of the hydrological conceptual model under investigation, representing a vertical 2D cross section driven along the middle line of the contaminant source zone. Right: schematic diagram of the impact of increased recharge on the concentration breakthrough curve (BTC) at an observation well located downgradient from the source zone. Modified from Libera et al. (2019).

The geochemical initial and boundary conditions in Table 2 are set to be the same as Bea et al. (2013), with a small modification of the nitrate-concentration initial condition for better matching with the observation. The pCO₂ concentration is based on previous publications (Bea et al., 2013) and assumed constant over the simulation, as increasing pCO₂ concentration has limited impacts on pH than the acidic species in the rain. Based on previous studies and field investigations, eight minerals are simulated in the reactive transport model in the F-Area. The dissolution and precipitation of initial minerals (i.e., quartz, kaolinite, and goethite) were modeled using kinetic-rate expressions derived from the literature and listed in Table 3. Gibbsite, jurbanite, basaluminite, opal, and schoepite are the species that can form when the plume interacts with the solids.

Table 3 Initial mineral volumetric fraction distribution in the simulation (Bea et al., 2013).

While Bea et al. (2013) implemented an electrostatic sorption model developed previously by Dong et al. (2012), which is less numerically efficient and requires additional parameterization, Arora et al. (2018) developed a non-electrostatic sorption model (NEM) at the F-Area site and demonstrated that the NEM achieved the same predictive performance as a surface complexation model (SCM) with electrostatic correction terms. The SCM approach is computationally expensive and requires the estimation of additional parameters that describe mineral surface characteristics. On the other hand, NEM does not consider the effects of the development of surface charge on the formation of surface complexes, and it also simplifies the parameters needed in the reactive transport modeling. In Arora et al. (2018), three mineral surface sites with different site densities and acidity constants are developed for modeling H⁺ sorption and transport and then further extended to noncompetitive and competitive H⁺ and U(VI). In this paper, we use the competitive H⁺ and U(VI) sorption NEM parameters (including site density and surface complexation constant listed in Table 4), which are derived from an inverse analysis and calibration by Arora et al. (2018), and implement them in the model.

Table 4 NEM model parameters for H⁺ and U(VI) competitive sorption (Arora et al., 2018).

3.3 CMIP5 climate scenarios

CMIP5 (Coupled Model Intercomparison Project, Taylor et al., 2012) is an experimental protocol with an ensemble of global climate model outputs to improve understanding of climate and to provide estimates of future climate change that will be useful to those considering its possible consequences. The climate forcing in our study used the $\mathrm{1}/\mathrm{8}$^∘ downscaled CMIP5 outputs at the F-Area study site from January 1950 to December 2100. The ensemble outputs include 28 models with four climate scenarios (RCP2.6, 4.5, 6.0 and 8.5) in the future climate projection. The top soil at the F-Area study site is sandy (Wainwright et al., 2014), so we assume that surface runoff is negligible. In other words, infiltration is calculated by subtracting ET from precipitation, which is simulated by the atmospheric and land surface models, respectively, from the coupled climate models. The figure below shows that the 10-year moving average of selected variables demonstrates that both precipitation and ET have increased by approximately 6 % since the 1950s to the present and will keep increasing up to an additional 6 % by the end of this century. The differences among climate scenarios are not statistically significant, but the highest greenhouse gas concentration (i.e., RCP8.5) ensemble simulates higher precipitation and ET than others. Although total recharge only slightly increases as both precipitation and ET are increasing (hence the difference offsets), our simulations focus on gaining understanding and quantitative estimation of changing climate impacts on the long-term robustness of contamination remediation in the F-Area.

3.4 Modeling scenarios

The modeling scenarios were developed based on Libera et al. (2019) and are only briefly described here. The modeling scenarios cover two stages of the F-Area historical operation and one additional stage in the future projection. The waste disposal was active during the period 1955–1988, and the basins were capped in 1988, when seepage from the basins into the vadose zone is assumed to have stopped. This study evaluates the effect of climate-change-induced variations in recharge on contaminant transport after 2020. A base case was developed with a constant recharge rate throughout the simulation period for assessing climate change impacts. The uniform recharge rate is 4.743 × 10⁻⁶ kg-water m⁻² s⁻¹ (0.15 m yr⁻¹ infiltration rate) based on the estimation in Bea et al. (2013). Furthermore, we developed three perturbed recharge scenarios with respect to the baseline recharge conditions. The three scenarios are (1) constant positive recharge shift from 2020, i.e., increasing precipitation scenarios, (2) constant negative recharge shift from 2020, i.e., decreasing precipitation scenarios, and (3) cap failure and a constant positive case from 2020. In both increasing and decreasing scenarios, recharge changes by 10 %–50 % after 2020. We hypothesize a complete failure of the containment structure scenario, which is represented by increased source-zone recharge of 10 %–50 % to the level of the surrounding region. Multiple studies have demonstrated that increased vegetation or other mechanisms could threaten or completely damage the integrity of the source-zone capping structure (Worthy et al., 2013, 2015). In addition to the perturbation scenarios, the contaminant transport and plume remobilization simulated by Amanzi are also forced by the four projection scenarios of CMIP5 ensemble climate model data, i.e., climate model scenarios. Instead of the constant recharge rate in those scenarios with changing precipitation, the annual recharge rate is used in CMIP5 climate scenarios from 1950 to 2100.

4 Results

4.1 Base case

The plume migration is depicted in Fig. 4 for the base case results described in Bea et al. (2013). The plume migrates through the vadose zone and then infiltrates vertically downward until it reaches the groundwater table (Fig. 4a). The plume then migrates vertically through the TCCZ into the LUTRA and also horizontally downstream closer to the FMB (Fig. 4b). Despite the low permeability of the TCCZ, leakage from the UUTRA to the LUTRA is observed over most of the flow domain. After basin closure and capping, the seepage from the basin is assumed to stop. The uncontaminated groundwater arriving from upgradient increased pH and reduced the U(VI) concentration (Fig. 4c). After the basin closure, because the vadose zone flow stops, pH remains low and U(VI) concentrations high in the vadose zone. In addition, the uranium concentration is higher in the TCCZ, where the permeability is low. The vadose zone below the basin appears to act as a long-term contaminant source for groundwater in the deeper layers (Fig. 4d). Although aqueous uranium concentration decreased by several orders of magnitude after the basin was capped, it is still higher than background concentration.

Figure 3 Simulated precipitation, evapotranspiration, and net infiltration (precipitation–evapotranspiration) at different climate projection scenarios from the CMIP5 datasets.

Figure 4 Plume profile of aqueous uranium concentration in the downstream of the F-Area study site from 1954 to 2020 in the base case simulation (the vertical exaggeration is ×5).

Figure 5 shows the base case breakthrough curves of pH, aqueous uranium, tritium, and nitrate at the source-zone well (FSB-95DR) and the downgradient well (FSB-110D) for the full simulation period (1956–2100). Both wells are located in the UUTRA layer. The simulated pH values rapidly decrease to 3.3 at both the source-zone well (Fig. 5a) and the downgradient well (Fig. 5c). In general, tritium concentrations (Fig. 5c) decrease more quickly and more dramatically than aqueous uranium and nitrate, owing to their radioactive decay. The uranium concentrations (Fig. 5b) increased from the background level 1.25 × 10⁻¹⁰ to 3.0 × 10⁻⁵ mol kgw⁻¹ (kilogram of water) at both wells in less than a few years and remained constant until basin closure in 1988. After the basin closure, pH rebounds to 4.0 in 2000 and gradually increases throughout the end of simulation. Similarly, uranium concentration (Fig. 5b) decreases by 2 orders of magnitude in 20 years and keeps decreasing to approximately 1.0 × 10⁻⁷ mol kgw⁻¹ by the end of the simulation period. Compared to the downgradient well, the source-zone well consistently has lower pH (Fig. 5a) and higher aqueous uranium (Fig. 5c) concentrations throughout the simulation period. By the end of 2100, pH (Fig. 5a) is higher than 5.0 and aqueous uranium concentration lower than 2 × 10⁻⁷ mol kgw⁻¹ (Fig. 5b) in most of the vadose zone at the source-zone well. Overall, our simulation has a similar fit to the observations by Bea et al. (2013), as the same parameters in Bea et al. (2013) are used.

Figure 5 Breakthrough curves of pH, aqueous uranium, tritium, and nitrate at the source-zone well and downgradient well in the base case over the simulation period (1954–2100).

4.2 Increasing recharge scenarios

The breakthrough curves under the increasing recharge scenarios are shown in Fig. 6. When recharge is increased, pH at the source-zone well (Fig. 6a) is significantly lower compared to the base case scenario. pH values are changed with different recharge rates as relatively high-pH-infiltrated rainwater dilutes the low-pH-contaminated environment in the subsurface system. However, the relationship between recharge and pH is nonlinear, with thresholds such that pH is lowest at +20 % recharge while pH is higher in the cases with +30 % to +50 % recharge. Nitrate concentrations at the source-zone well (Fig. 6b) increase immediately after 2020 and spike 5 years after perturbation, with the highest concentration in the greater recharge (+50 %) case. After 2050, nitrate concentration is highest with +20 % recharge and decreases from +30 % to +50 % recharge (Fig. 6b). The tritium concentrations (Fig. 6c) peak similarly to nitrate, although tritium decreases significantly after 2040 due to radioactive decay. The uranium concentrations (Fig. 6d) are also similar to the breakthrough curves for the nitrate concentrations and show negative correlation with pH oscillation. At downgradient locations, pH (Fig. 6e) is not influenced by the recharge increase up to +30 %. Above the 40 % increase, pH decreases significantly after 2040. Nitrate concentrations at the downgradient well (Fig. 6f) decrease immediately after 2020 due to dilution but increase afterwards, with peaks around 2040. Similarly to the source-zone well (Fig. 6b), the concentration peaks are higher with greater recharge rates and remain higher than the base case throughout the end of the simulations. The tritium concentrations (Fig. 6g) keep decreasing after 2020 with the peaks in 2040, with the similar behavior of sudden increase showing in nitrate concentration (Fig. 6f) and higher concentrations in the high recharge scenarios. The uranium concentrations (Fig. 6h) also exhibit patterns similar to those of nitrate (Fig. 6f) in that both the peak and remaining concentrations are higher in the greater recharge scenarios.

Figure 6 Breakthrough curves of pH, nitrate, tritium, and aqueous uranium at the source-zone well (a–d) and downgradient well (e–h) in the base case and increasing precipitation scenarios from 2020 to 2100.

The reactive (uranium) and non-reactive (nitrate) species are compared in Fig. 7. K_d values are computed by sorbed uranium concentration in the solid phase with the aqueous uranium concentration from the model outputs. Figure 7a shows that K_d values at the source-zone well are lower in the increasing recharge cases than the base case, which is consistent with pH breakthrough curves (Fig. 6a). The +20 % case has the lowest K_d at the source-zone well, while the K_d values are higher in the smaller recharge case (+10 %) and greater recharge cases beyond +30 %. In contrast, at the downgradient well (Fig. 7d), the K_d values are lower in the +40 % to +50 % scenarios, echoing the downgradient pH breakthrough curves in Fig. 6e. In addition, we compare uranium and nitrate concentrations with respect to the maximum concentration (i.e., the peak concentrations that occur after a few years in the increasing recharge scenarios) as well as the average concentration from 2040 to 2100, which illustrates the long-term contamination trend. Figure 7b–c present the ratio of uranium and nitrate, defined as the concentration in each scenario compared to the baseline case. In the maximum concentration at the source-zone well (Fig. 7b), the ratios are mostly higher than 1, demonstrating that the maximum concentration is higher in the greater increasing recharge scenarios. The uranium maximum concentration ratio is higher than the nitrate; therefore, the increasing recharge affects the uranium concentrations more than the nitrate concentrations at the peaks (Fig. 7b). For average concentrations at the source-zone well (Fig. 7c), the ratio increases in the +20 % recharge case but decreases at greater recharge values. Different from the maximum concentrations, the mean uranium ratio becomes lower than the mean nitrate ratio and falls below 1.0 in the greater recharge scenarios (Fig. 7c). At the downgradient well, the maximum concentration ratios are less than 1.0 in the (+10 % to +30 %) recharge scenarios but higher than the base case in greater recharge (+40 % to 50 %) scenarios, while nitrate and uranium ratios are similar (Fig. 7e). The average concentration ratios at the downgradient well after 2040 are generally higher with increasing recharge and reach their highest values at the +40 % scenario (Fig. 7f). The nitrate concentration ratios are lower than uranium in the smaller (+10 % to +30 %) recharge scenarios but are higher in those scenarios of above +40 % recharge.

Figure 7 Breakthrough curves for K_d at the source-zone well (a–c) and the downgradient well (d–f) for the increasing recharge scenario from 2020 to 2100. Maximum and average ratios of base case to increased recharge case for uranium and nitrate concentrations at both well locations.

4.3 Decreasing recharge scenarios

Although decreasing recharge has little impact on pH at the source-zone well up to −30 % (Fig. 8), pH increases significantly in the −40 % to −50 % recharge scenarios. The nitrate concentrations (Fig. 8b) increase immediately after the perturbation of recharge and then decrease throughout the end of the simulation. Similarly to pH, nitrate concentrations (Fig. 8b) do not change significantly in smaller decreasing recharge scenarios but decrease by 2 orders of magnitude in the greater (−40 % to −50 %) decreasing recharge scenario. Tritium concentrations (Fig. 8c) also increase immediately after 2020 and then decrease; the rate of decrease is more rapid than the nitrate concentrations due to radioactive decay and exhibits few differences among decreasing recharge scenarios. The uranium concentration (Fig. 8d) breakthrough curves are similar to the nitrate curves. At the downgradient well, the pH values have a similar trend to the source-zone well in all decreasing recharge scenarios before 2040. However, the breakthrough curves diverge after 2040 and increase more in the greater decreasing recharge scenarios. The pH values are higher than the source-zone well and reach as high as 7.0 in the −50 % recharge scenario in 2100 (Fig. 8e). The nitrate concentrations in the downgradient well (Fig. 8f) keep decreasing in the first 10–15 years after 2020. Concentrations peak around 2025–2035; the decrease is more significant in all the decreasing recharge scenarios than the base case. In general, the peak concentrations occur earlier and higher in the greater decreasing recharge scenarios, and the breakthrough curves decrease more quickly and lower in the long-term projection to 2100. Spikes were observed in the tritium concentration breakthrough curves (Fig. 8g) as well as in smaller magnitudes at the downgradient well 10–15 years after the perturbation. The uranium concentration breakthrough curves (Fig. 8h) are similar to the nitrate but decrease more rapidly in all cases.

Figure 8 Breakthrough curves of pH, nitrate, tritium, and aqueous uranium at the source-zone well (a–d) and downgradient well (e–h) in the base case and decreasing precipitation scenarios from 2020 to 2100.

K_d breakthrough curves generally reflect the pH breakthrough curves in Fig. 8 and are higher in the decreasing recharge scenarios at both well locations. Figure 9a shows that at the source-zone well, the base case ∼ −30 % cases have relatively similar K_d values throughout the simulation period. After 2060, the −40 % and −50 % recharge scenarios both significantly increase. At the downgradient, K_d values are generally higher than the source-zone well, and the differences in K_d values among cases are more pronounced (Fig. 9d). The K_d value is lowest in the higher recharge (base case and −10 %) scenarios and largest in the significantly decreasing recharge (−50 %) case. The −10 % and −20 % recharge scenarios significantly diverge at 2040 and converge at 2100. Similarly to the increasing recharge scenarios, maximum uranium and nitrate concentrations at the source-zone well occur immediately a few years after the perturbation (Fig. 8). With decreasing recharge from −10 % to −50 %, maximum concentration ratios are higher than 1.0 and increase with decreasing recharge, while average concentration ratios are generally lower than 1.0. In Fig. 9b, the uranium maximum concentration ratios are slightly higher than nitrate, with a greater difference in the −50 % recharge case. In Fig. 9c, the ratios of long-term average concentrations show that both uranium and nitrate concentrations are nearly the same as the base case in smaller decreasing recharge scenarios (−10 % to −20 %) but decrease quickly and are significantly lower in the greater decreasing recharge scenarios (−30 % to −50 %). Compared to the results at the source-zone well, the maximum and average concentration ratios at the downgradient well (Fig. 9e and f) have similar trends. Nitrate maximum concentration ratios are higher than nitrate (Fig. 9e), and their differences are greatest in the −20 % and −30 % recharge cases. The average concentration ratios (Fig. 9f) decrease with decreasing recharges, and uranium ratios are consistently higher than nitrate.

Figure 9 Breakthrough curves for K_d at the source-zone well (a–c) and the downgradient well (d–f) for the increasing recharge scenario from 2020 to 2100. Maximum and average ratios of base case to increased recharge case for uranium and nitrate concentrations at both well locations.

4.4 Cap-failure scenarios

In the cap-failure scenarios, pH is always lower than the base case across +10 % to +50 % recharge rates (Fig. 10a). At the source-zone well, these pH values dip below 3.5 in 2030, rebound to 4.0 after 2045, and then slightly increase to 4.3 by the end of the simulation. The +50 % cap-failure scenario has the highest pH value compared to the +10 % to +40 % cap-failure cases. Nitrate concentrations spike and increase by 1 order of magnitude in 2030 and then decrease to the same level as the base case in 2050 (Fig. 10b). The patterns of tritium and uranium breakthrough curves (Fig. 10c–d) look very similar to nitrate. Among the breakthrough curves of nitrate, tritium, and uranium across all cap-failure scenarios, the +50 % cap-failure scenario simulates the earliest peak, while the +20 % scenario simulates the highest peak. At the downgradient well, pH values at all cap-failure scenarios increase with the base case in the first 10 years and then decrease around 2035 and remain lower than the base case (Fig. 10e). pH values only decrease from 5.5 to 5.0 in the smaller +10 % to +20 % recharge rates. However, they decrease significantly with greater recharge rates (+30 % to +50 %) in those cap-failure scenarios. The breakthrough curves of pH increase in the first 10 years after 2020, dip to 3.6 around 2040, and then slightly increase but require several decades to rebound to the same pH level as in 2020. Compared to the nitrate concentration breakthrough curves at the source-zone well (Fig. 10b), the peaks at the downgradient well are simulated in 2040 with a 10-year delay (Fig. 10f). The nitrate concentrations in those greater (+30 % to 50 %) recharge rates occur earlier and are higher than in the smaller (+10 % to 20 %) recharge rates. The tritium concentration shows similar peaks to nitrate, but the earliest peak with +50 % cap failure has the highest values, and the later peaks with smaller (+10 % to 20 %) recharge rates will be lower because of tritium radioactive decay (Fig. 10g). Uranium concentration breakthrough curves show similar behaviors to nitrate in both wells (Fig. 10b, d and f, h).

Figure 10 Breakthrough curves of pH, aqueous uranium, tritium, and nitrate at the source-zone well (a–d) and downgradient well (e–h) in the base case and cap-failure scenarios from 2020 to 2100.

Kd breakthrough curves are highly correlated with pH at both monitoring wells (Fig. 10a, e), and K_d decreases by 2035 and 2040 at both wells, respectively, returns to the 2020 level around 2050, and then keeps increasing until 2100. At the source-zone well, K_d values decrease and rebound most quickly in the +50 % recharge case, and the smallest +10 % recharge rate case shows a similar trend but is delayed by nearly 10 years (Fig. 11a). At the downgradient well, the K_d breakthrough curves at higher recharge cases (+30 % to +50 %) are more closely correlated with pH and decrease around 2040, while smaller recharge cases (+10 % to +20 %) are more similar to the base case (Fig. 11d). A turning point occurs in 2040, when the +30 % case switches places with the +50 % case and has the lowest K_d value until 2100, similarly to the behavior of aqueous uranium breakthrough curves in 2040 (Fig. 10h). When comparing Figs. 11d and 10e, it is clear that although pH is not the highest in the +20 % cap-failure scenario, after 2070, that scenario has the highest K_d value and more adsorption. In cap-failure scenarios, the maximum and average uranium concentration ratios are consistently greater than nitrate in both wells and follow the same trend with increasing recharge rates (Fig. 11b–c, e–f). Both ratios of uranium and nitrate maximum and average concentration are 1 order of magnitude greater than the base case. The maximum concentrations of uranium and nitrate are observed in 2030 and 2040 at the source-zone well and downgradient well, respectively (Fig. 10b, f), although it is difficult to tell the difference from the breakthrough curves because of the magnitude of peak concentrations. The uranium and nitrate maximum concentration ratios are highest in the 20 % cap-failure scenarios (Fig. 11b) and decrease with a greater increasing recharge rate. The ratios of uranium average concentrations against the base case are also persistently higher than nitrate in the long term throughout 2100 and decrease with greater recharge rate (Fig. 11c). At the downgradient well, the maximum concentration ratio against the base case generally increases with greater recharge rate and is largest in the +40 % recharge case (Fig. 11e). The average concentration ratio increases with the smaller (+10 % to +30 %) recharge rates and then decreases with the greater (+40 % to +50 %) recharge rates (Fig. 11f).

Figure 11 Breakthrough curves for K_d at the source-zone well (a–c) and the downgradient well (d–f) for the increasing recharge scenario from 2020 to 2100. Maximum and average ratios of base case to increased recharge case for uranium and nitrate concentrations at both well locations.

4.5 Climate model scenarios

Recharge rates are calculated by subtracting ET from precipitation in the four CMIP5 climate projection scenarios. The highest greenhouse gas concentration pathway RCP8.5 scenario has the maximum simulated precipitation and evapotranspiration. However, the differences in recharge rate are small across those four scenarios as both precipitation and ET increase in the projection (Fig. 3). Therefore, the concentration breakthrough curves are similar under those climate scenarios. The average recharge rate in those scenarios is around 8.0 × 10⁻⁶ kg-water m⁻² s⁻¹ (0.253 m yr⁻¹) or approximately 1.68 times higher than the base case. In general, simulated contaminant concentrations in those climate scenarios are lower than the base case due to dilution effects with a greater recharge rate, except that pH values are also lower than the base case (Fig. 12). The breakthrough curves decrease more quickly before 2020 (not shown in Fig. 12) and reach background concentration sooner than the base case.

At the source-zone well, the pH breakthrough curve gradually rebounds from 4.0 to 4.5 by the end of the simulation (Fig. 12a). Both nitrate and uranium concentrations show annual variability after 2020, as recharge rates are changing annually (Fig. 12b, d). Specifically, nitrate breakthrough curves (Fig. 12b) become steady state sooner than the uranium, as the nitrate background concentration is higher. The oscillation is hardly observed in tritium concentration breakthrough curves, as it decreases more quickly due to decay. At the downgradient well, pH values across climate scenarios are consistently lower than the base case with annual variability (Fig. 12e). Compared to the results at the source-zone well, the nitrate concentrations at the downgradient well (Fig. 12f) are lower than the background level with greater annual variability and become steady state a few years later. The tritium concentration becomes extremely low below 1.0 × 10⁻¹⁵ mol kgw⁻¹ (Fig. 12g), while uranium concentrations return to the background level after 2030 (Fig. 12h).

Figure 12 Breakthrough curves of pH, nitrate, tritium, and aqueous uranium at the source-zone well (a–d) and downgradient well (e–h) in the climate scenarios from 2020 to 2100.

5 Discussion

A balance between dilution and remobilization is a key factor determining the contaminant concentration depending on recharge rates, as discussed in Libera et al. (2019). In the increasing recharge scenarios, contaminant concentrations decrease first due to dilution and then increase because the mobilized contaminants migrate from the source zones to the wells. The highest recharge scenario has the earliest and highest peak in contaminant concentrations due to a stronger remobilization effect, but it has the lowest concentrations and highest pH later due to dilution. In the later period, the increasing recharge again causes dilution due to flushing, resulting in a concentration level below the base case. Because of long-term dilution, the aqueous uranium concentration in greater increasing recharge scenarios is even lower than the base case at the source-zone well after 2035. The relationships between concentrations and recharge are nonlinear and nonmonotonic, depending on different times and locations. Changing recharge rate has less impact at the downgradient well where the spikes are delayed for approximately 10 years since its location is further from the seepage basin, and it takes time for the remobilized plume to reach it. The breakthrough curves of smaller (+10 % to +30 %) increasing recharge scenarios are similar to those of the base case, with slight dilution effects throughout 2100, while concentration spikes due to remobilization in 2040 are observed with the larger (+40 % to +50 %) increasing recharge scenarios (Fig. 6e, h).

In the early stage of decreasing recharge scenarios, contaminant concentrations increase because of diminished flushing and a low flow rate of clean groundwater. Later in the simulation period, contaminant concentrations decrease significantly in the greater decreasing recharge scenarios, when the groundwater table declines and isolates the residual contaminants in the vadose zone. This was not observed in the previous tritium simulation (Libera et al., 2019) because of tritium's radioactive decay. In general, this means that decreasing precipitation and droughts are effective in sequestering contaminants in the vadose zone. At the same time, it implies less flushing and an increase in residence time of the contaminants at the site. The uncontaminated groundwater from the upgradient also migrates more slowly in the aquifer. The larger volume of residual contaminants could potentially increase risk, particularly considering extreme precipitation events, which are projected to happen more frequently in many climate models (USGCRP, 2017). Also, there is more interest in groundwater resources during a drought, which leads to increased pumping in the contaminated aquifer. Although such pumping activities are strictly regulated at our study site, such trade-offs require attention at other sites.

To investigate the impact on reactive species such as uranium, we compared reactive (uranium) and nonreactive (nitrate) concentration ratios to assess the impacts of reaction and sorption.

We originally hypothesized that increasing recharge would decrease reactive species concentration further, since increasing the volume of water in the domain would increase pH, which limits the mobility of uranium. However, Fig. 7 shows that the uranium-concentration ratios compared to the base case increase more significantly than the nitrate concentrations. This is because the remobilization occurs when the pH is still low and also because remobilization happens to both uranium and protons (Fig. 6). In addition, the amount of residual contaminants is larger for uranium than nitrate due to sorption. Later in the simulation period, the uranium average concentrations are lower than for nitrate and decrease with greater recharge scenarios, because increasing pH, due to long-term dilution by additional recharge, immobilizes uranium.

In cap-failure scenarios, sorption of uranium is reduced with increasing infiltration, because K_d is sensitive to lower pH due to remobilization through the basin. At the downgradient well, the greater recharge cases (+30 % to +50 %) have a more closely correlated K_d and pH and have a higher aqueous uranium concentration than the lower recharge scenarios. In our scenarios, there is a clear change in the balance of aqueous and sorbed uranium concentration in the transition from +20 % to +30 % recharge, where the system's sorption in the downgradient fundamentally changes. The cap-failure cases indicate that changing recharge and cap-failure levels can trigger dramatic changes in pH and sorption. Similarly to Libera et al. (2019), this study confirms the importance of cap or surface barriers for limiting the impacts of cap failure under extreme climate regimes. Uniquely, the uncertainty of K_d value was constrained to 10² to 10³ in our study compared to the greater range (10²–10⁶) in the previous studies (Bea et al., 2013).

At this site, pH and uranium concentration fronts are retarded, because they are affected by the adsorption and ion-exchange processes onto kaolinite and goethite (Bea et al., 2013). Limited sites for adsorption/exchange are saturated by the elevated H⁺ and uranium loading, and their concentrations eventually reach steady state at the end of the projection period. Overall in our scenarios, the change in recharge has a similar impact on uranium and nonreactive species, which is largely attributed to pH buffering due to mineral precipitation. The increase in pH due to dilution encourages the precipitation of kaolinite, but the precipitation reaction of kaolinite produces H⁺ ions, which then decreases pH. At low pH, the hydroxyl groups on the octahedral structures of aluminosilicates like kaolinite become protonated, effectively creating a net positive charge on the mineral. This means that uranium cannot sorb to the clay and is therefore mobile in the system. Previous experimental (Dong et al., 2012) and modeling (Arora et al., 2018) studies also reaffirmed that percent U(VI) sorption is greater with a higher, neutralized pH, because U(VI) and H⁺ are competing in sorption. This is the process of dissolution and precipitation of kaolinite.

$$\begin{array}{}\text{(1)}& {\mathrm{Al}}_{\mathrm{2}}{\mathrm{Si}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}(\mathrm{OH}{)}_{\mathrm{4}}+\mathrm{6}{\mathrm{H}}^{+}\leftrightarrow \mathrm{2}{\mathrm{Al}}^{\mathrm{3}+}+\mathrm{2}{\mathrm{SiO}}_{\mathrm{2}}+\mathrm{5}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\end{array}$$

A similar reaction occurs with gibbsite. Dong et al. (2012) showed that there was an insignificant weight percent or volume fraction of gibbsite at F-Area, since it only forms at pH > 5.4. However, in the decreasing recharge scenarios, all the recharge cases at the downgradient well have a pH between 5.4 and 7 after 2070, and pH at the two greater recharge cases at the source-zone well also surpasses 5.4. Decreasing recharge would likely trigger the formation of gibbsite, which could increase pH buffering. Additionally, according to Bea et al. (2013), this mechanism as well as cation-exchange and adsorption processes on kaolinite and goethite explain some buffering of pH. The pH buffering effect is the major mechanism for pH remaining low for an extended period of time in climate resilience studies with reactive transport modeling. Nevertheless, our model is built upon more than 10 years of site characterization, sorption experiments, and reactive transport models (Dong et al., 2012; Bea et al., 2013; Sassen et al., 2012; Wainwright et al., 2019). The NEM sorption model used in this study is based on Arora et al. (2018), which is calibrated with long-term monitoring datasets and considered competitive H⁺ and uranium sorption.

In addition to understanding the impact of a range of recharge scenarios, this study has established a pipeline to use the CMIP5 climate model projections as input to the hydrology and reactive transport modeling simulations. Although increasing precipitation is projected over time, we found that the increasing ET associated with temperature can reduce the recharge rates. We found that, compared to the base case and hypothetical scenarios, the CMIP5 climate data project a small increase or no change in recharge rate over time, indicating that the changing climate has minor effects on the contamination plume and breakthrough curves in our study site. This is similar to the behaviors observed in the increasing precipitation scenarios in Fig. 6, that smaller recharge increases have little impact on the concentration breakthrough curves, because the increasing recharge is below the threshold that may cause significant remobilization. Contaminant migration is more controlled by the transport process. Our reactive transport modeling with CMIP5 projection recharge shows that contaminant migration is sensitive to recharge rate. In our study, ET is prescribed from the ensemble average of CMIP5 datasets and is not computed in our simulations. The annual variability of precipitation and ET, in other words net infiltration, is more significant in RCP8.5 than other scenarios. The uncertainty in the estimation of ET as well as the annual variability in CMIP5 scenarios could significantly affect the assessment of waste disposal and contaminant transport.

6 Conclusion

The climate resilience of residual contamination at the SRS F-Area waste disposal site throughout the projection period from 2020 to 2100 is investigated in this study. Groundwater flow, mineral reactions, surface complexation sorption, and ion-exchange processes are simulated by the Amanzi and PFLOTRAN flow and reactive transport model. We illustrate four scenarios characterized by a range of variable recharge values: (1) increasing recharge after 2020, (2) decreasing recharge after 2020, (3) cap failure and constant positive recharge shift, and (4) recharge rate under different RCP scenarios from the CMIP5 climate model projection. Although exaggerated in the first three cases, this systematic study using changing recharge rates was useful in identifying the phased impacts of increasing or decreasing recharge rates as well as the difference between the reactive and nonreactive species. Plume distribution and breakthrough curves of chemical species are evaluated to assess the impacts of changing recharge rate and flow conditions. The ratios of maximum and average reactive and nonreactive species concentrations between scenarios and base case are used to understand how climate change affects the adsorption and ion exchange of residual contaminants in the subsurface domain. Furthermore, K_d breakthrough curves are evaluated to understand the pH effects on sorption with different recharge rates in those scenarios.

With increasing recharge rates, pH decreases and residual contaminant concentrations increase because of the remobilization of protons and reactive species. The impact on uranium or pH-dependent species is the same as nonreactive contaminants. K_d values are correlated with pH and behave differently when changing recharge rates beyond certain thresholds. In most cases, uranium-maximum concentration ratios against the base case are higher than the nitrate concentration ratios, owing to remobilization, while the uranium concentration breakthrough curves in the later period depend on long-term flow conditions. The results of cap-failure scenarios suggest that reactive transport modeling and analysis of pH effects on reactive species are important for the risk assessment of such engineering failures.

Our results highlight that climate change impacts may not be intuitive and must be analyzed quantitatively by models. ET projection has great uncertainty but is particularly important in determining the recharge rates in reactive transport modeling for climate resilience studies. Reactive transport models which consider pH dependency for reactive species are essential for analyzing the impacts of pH with changing recharge rates. Although this study is focused on one site, we developed the pipeline to use climate projection datasets in reactive transport modeling and thereby demonstrated the capability of assessing climate change impacts on waste disposal sites. We expect that our approach and insights are transferable to other sites that have large amounts of residual contaminants in the vadose zones or in the groundwater.