Site Geologic Units

Within the study area, there are several mapped geologic units that comprise the surficial fractured rock and karstic groundwater aquifer system. These units can be grouped into structural features that form the Blue Ridge, Frederick Valley synclinorium, and Triassic basins.

Boring Log Reanalysis

The records produced during well drilling activities—including drillers’ boring logs and descriptions produced by onsite geologists—provide a valuable primary documentation of the subsurface conditions. At Fort Detrick Area B, a total of 76 drilling log records were reported in appendix A of the draft final remedial investigation report (Arcadis Inc., written commun., December 2019). These records were summarized to identify (1) the thickness of the reported soil overburden by recording the first occurrence of competent bedrock, (2) the number of voids or cavities encountered during drilling, (3) the reported depth of the deepest void or cavity, and (4) the degree of reported weathering and depth to which this weathering occurred. Using this information, the presence or absence of an epikarst zone was interpreted and the depth of the bottom of the zone was estimated. The determination of the epikarstic depth was not conducted for 16 wells whose drilling logs were too short or not complete enough to perform an analysis. The bottom of the epikarstic zone was identified as the deepest reported void zone and (or) the deepest zone of substantial weathering before a transition into unweathered fracture-dominated bedrock. The results of these analyses were mapped to identify spatial patterns within the subsurface data. Due to the heterogeneity of the subsurface within karst systems and the limited spatial extent of an individual borehole, interpolation between boreholes was not performed.

Geophysical Log Reanalysis

Geophysical log data were collected in 23 geophysical logs from 13 wells by ARM Geophysics of Hershey, Pennsylvania while a series of new wells were constructed in 2011 and 2012. The locations of the geophysical logs, which are named with “BMW” followed by a 2-digit identifier are shown in figure 9. The names provided within appendix D of Arcadis Inc. (2014) are retained in this report. Individual logs were sometimes collected at the same location but at different depth intervals; the letters A–E following the location name denote the relative depth interval of the data collection. The data were reviewed and the optical and acoustic televiewer (OTV and ATV, respectively) data were independently analyzed using WellCad version 5.2 to identify the orientation of bedding, fractures, and voids. Imagery data collected during logging included OTV and ATV. Line data collected during logging included 3-arm caliper, gamma radiation, fluid resistivity, fluid temperature, and formation resistivity. These data files and well logs are available in appendix D of Arcadis Inc. (2014). The raw data logs were reinterpreted by the USGS to evaluate the orientation of features that commonly act as preferential flow conduits, including bedding, voids, and fractures.

The following steps were undertaken in WellCad to prepare the data and to collect the strike and dip of bedding features and fractures within the geophysical logs. In several wells, multiple geophysical logging events occurred in the same borehole at different depth intervals; these data were assimilated for interpretation. Linear data collected at the site were inspected to ensure that variations in gamma and resistivity were within natural instrumental bounds. In some cases, the resistivity or gamma values were spurious and were not interpreted. Where multiple readings were collected in the same borehole at the same interval, the data were overlain to evaluate the consistency of the observations.

The imagery data was primarily used to make designations of strike and dip of fractures and bedding. Several processing steps were necessary to ensure the accuracy of strike and dip measurements. The instruments used to collect the data use Earth’s magnetic field to orient to magnetic north. Therefore, a rotation is required to transform to true north. A 10.5° clockwise rotation was applied to all image logs to ensure their true north orientation, as the declination at Fort Detrick is 10.5° W. In some cases, the OTV or ATV logs appeared offset from each other or the depth interval collected during one instrument run differed from the other. In this case, common features of both logs were used to apply a stretch or shrink to one log to ensure a consistent depth observation. The gamma and caliper logs assisted in making this determination.

When the imagery logs were properly oriented and aligned, linear features in the logs were identified and classified. As the geophysical logs represent the imagery of the interior of a cylinder, linear features appeared as sinusoids in geophysical logs. Fractures and bedding were interpreted by considering variations in grain size, rock color, variations in the caliper log, and acoustic travel time and intensity. There were many intervals where cloudy sediment impeded the view of the optical log and only ATV logs were used to identify fractures. Large void spaces were encountered in many logs. The orientation of these features was estimated by fitting sinusoids to the sides and middles of the voids; however, the orientation of these features is less certain than the linear features.

Once these features were classified, several corrections were applied to ensure the strike and dip represented the true rock units. First, a correction for borehole deviation was applied by using the tilt and azimuth provided in the ATV logs. This correction is required because drilled boreholes are typically not perfectly vertical. Second, the strikes and dips were corrected for variations in borehole diameter. The diameter may vary from the drilled diameter, which will affect the true dip. Once these corrections were applied, the strikes and dips were determined. Polar plots indicating the poles to the plane of these features projected onto the southern hemisphere of a stereonet plot allow the features to be examined on a well-by-well basis.

Base-Flow Separation

Hydrograph separation was used to differentiate the observed Carroll Creek stream flow into two components—base flow and quickflow. Stream base flow represents discharge from the groundwater systems into the stream-channel network. Stream quickflow represents short residence time water—including surface runoff or shallow lateral subsurface flow—that reaches the stream network following precipitation events. Seasonal and interannual variability in stream base flow reflects changes in groundwater storage. There are many techniques that are designed to separate the base flow component from the quickflow component of stream discharge; many of these are graphical in nature and use the shape of discharge time series to distinguish between two or more components, often base flow and quickflow discharge. These graphical hydrograph separation techniques, although reproducible, do not have a physical basis and are not amenable to estimation of uncertainty.

Chemical hydrograph separation allows for a more physically based determination of the relative influence of base flow and quickflow and incorporates stream chemistry into the algorithm. Specific conductance (SC) is often used as a proxy for groundwater-derived total dissolved solutes in stream, as precipitation typically has low SC whereas groundwater typically has elevated SC resulting from its longer contact time with aquifer materials. In this study, we use optimal hydrograph separation (OHS; Raffensperger and others, 2017; Foks and others, 2019), a method of chemical hydrograph separation that employs the use of a two-parameter recursive digital filter. The two dimensionless parameters are a recession constant (α) and the maximum value of the base-flow index (β). The value of α is calculated using streamflow recession analysis and the value β is optimized using a mass balance constraint. OHS has the following key features

Two-component hydrograph separation, or separation of base flow from quick flow, for the chosen sites was performed using a recursive digital filter (RDF). Eckhardt (2005) proposed the following RDF to estimate the base flow component of streamflow:

where Base flow at time step j is restricted to $Q_{B_j}\le Q_j$. The parameter β is identical to the parameter BFI_max in Eckhardt (2005).

The recession constant parameter, α, was derived for each study site using the log-linear slope of the recession hydrograph for storm events within the period of record, using the methods of Rutledge (1998). An α value was calculated for each recession period of sufficient length and then the median α from all storm recessions was taken as the α used in the OHS calculations. Only α values computed from significant (p>0.05) log-linear regressions were retained. The minimum number of days for α to be calculated was determined using the formula of Linsley and others (1982), which states that the number of days (N) is proportional to the drainage basin area in square miles (A):

The second adjustable parameter in the two-parameter RDF, β, is a fitted, unmeasurable parameter that ranges from 0 to 1 (Eckhardt, 2005; Raffensperger and others, 2017). For sites without available SC data, the model parameter β was determined using a method developed by Collischonn and Fan (2013), referred to as “CaF.” This method for estimating β is based entirely on discharge records. The method applies a backward-moving filter using a previously calculated recession parameter. Collischonn and Fan (2013) tested the filter using data from 15 streamflow sites in Brazil with differing physical characteristics and showed that the β values obtained were comparable to the predefined values suggested by Eckhardt (2005) based on the class of aquifers encountered in each basin.

For sites with SC data, OHS was applied, in which model parameter β is optimized by applying a SC mass balance constraint (Rimmer and Hartmann, 2014; Raffensperger and others, 2017). Streamflow is initially separated into two components, base flow (Q_B) and quick flow (Q_S), while values of SC are proposed for the base flow (C_B) and quick flow (C_S) components. These variables, as well as streamflow (Q), are used to estimate stream SC (C_sep) by

The root-mean-square error, E, between the estimated stream SC and observed stream SC (C_obs) is calculated as This process is repeated until there is minimized error, E, implying that the parameter value for β is also optimized (Rimmer and Hartmann, 2014; Raffensperger and others, 2017). Optimization was performed in version 3.4.3 of R (R Core Team, 2018) using the algorithm “Bound Optimization by Quadratic Approximation” (BOBYQA) in the R package nloptr, version 1.0.4 (Powell, 2009; Ypma, 2018).

Two OHS models were used to estimate base flow SC (C_B) (Raffensperger and others, 2017). The first OHS model estimated SC by way of a sine-cosine function over time (sin-cos model) to emulate seasonal variation and is defined as

in which the base flow SC on day j ($C_{B_j}$$C_B{}_j$) is a sum of sine and cosine functions of time (t_j) with an annual period, described by a mean value (${\overline{C}}_B$), amplitudes ($C_B^{*s}$,$C_B^{*c}$), and a beginning time/day (t₀). In the optimization, the following values are estimated that minimize the root-mean-square error (RMSE; eqs. 1–3): β, ${\overline{C}}_B$, $C_B^{*s}$, t₀, $C_B^{*c}$, and C_S.

The second OHS model estimates C_B by using a peak-identification algorithm and linear interpolation (SCfit model). Peaks in the observed SC data were identified using the function “findpeaks” in the R package pracma (Borchers, 2021), whereas C_B values were estimated with linear interpolation between the identified peaks. The SCfit model has only two variables to optimize (β, ). SCfit and sin-cos models were accepted if the Nash-Sutcliffe efficiency coefficient (NSE) (Nash and Sutcliffe, 1970) was greater than 0.3 and did not converge to a user-defined optimization bound (lower bound=0.00001; upper bound=1.0) (Raffensperger and others, 2017; Foks and others, 2019). The NSE cutoff value of 0.3 is the same as that used in Raffensperger and others (2017).

Comparing the base flow separation results to other regional watersheds provides insights on the groundwater-surface water exchange processes occurring at Fort Detrick Area B and the applicability of conceptual models developed in other hydrogeologic settings. The intercomparison of base flow separation will focus on two metrics calculated for all available streamgages for the period of October 1, 2017, to September 30, 2020, a period of 3 water years. These metrics are (1) groundwater base flow per unit area (in in/yr) and (2) the base flow index, which is the proportion of base flow to total streamflow and ranges from 0 to 1. A base flow index closer to 0 is expected in basins with flashy storm responses (flashiness describes short-lived and rapid storm hydrograph peaks) and a base flow index closer to 1 is expected in basins with muted storm responses and steady inter-storm flow sustained by groundwater.

The two primary base flow separation techniques discussed in this report are OHS and the CaF method. An additional widely applied method for base flow separation is a graphical method called PART, which is described by Rutledge (1998). This third method is included for comparison because it is more widely used within regional studies of karst hydrogeology and it is readily applied to all available streamgages using only a time series of discharge. The part() function of the USGS DVStats package version 0.3.4 was used to perform this analysis on the daily mean streamflow values stored in the USGS NWIS database (Lorentz and DeCicco, 2021). To apply the method, gaps of 3 days or less in the streamflow record were filled using linear interpolation.

To interpret the base flow separation results, a binary classification of lithology was used to identify watersheds with a potential for influence by karst processes. Using geologic data reported by the USGS state geologic map compilation (Horton and others, 2017), watersheds were classified as “karst” if more than 50 percent of the watershed area was underlain by carbonate geologic units. Although a simple classification scheme, this designation captured many of the USGS streamgages reported in the regional literature to be located in karst environments in Virginia, West Virginia, and Maryland.

Groundwater Gradients

In karst environments, hydraulic gradients, the number of active conduits, flow routes, and directions of groundwater flow may change rapidly with changing hydrologic conditions (Taylor and Greene, 2008). At Area B, previous studies have produced potentiometric surface maps (fig. 4). Although in karst environments potentiometric mapping may not fully characterize preferential flowpaths arising from anisotropy and highly heterogenous subsurface processes, it is still a useful tool for guiding other types of analyses, including the design and interpretation of groundwater tracer tests (Taylor and Greene, 2008). However, potentiometric mapping is limited in its ability to characterize short-lived responses to hydrologic events. The groundwater gradient, along with hydraulic conductivity and cross-sectional area, determines the groundwater flow magnitude along a groundwater flow path. In karst systems, conduits and preferential flow paths may have extremely low hydraulic conductivity. In this case, very low gradients may still convey large flow volumes.

In this analysis we focus on temporal changes in the observed groundwater gradient between wells in the hydrologic monitoring network to characterize the short-lived responses to hydrologic events. Changes to the flow direction and hydraulic conductivity resulting from the complexity of the subsurface conduit system are possible, though the approach we take assumes that the changes in flow direction and hydraulic conductivity through time are negligible. We examine temporal changes to hydraulic gradients between individual well pairs. The selection of these pairs is guided by the potentiometric surface mapping. The magnitude of the gradients is not intended to guide estimates of total groundwater discharge for the reasons previously discussed.

Groundwater gradient time series were computed from the daily mean water levels. Horizontal groundwater gradients between well pairs were computed by dividing the difference in water level elevation by the Euclidean distance, neglecting the influence of topography and observation depth. Vertical groundwater gradients between well pairs were computed by dividing the water level elevation difference by the difference in observation depth below land surface.

Soil-Water Balance Modeling

Accurate estimates of the temporal and spatial variations in water budget components require accounting of water inputs and removal from the hydrologic system. Input to the groundwater system in Pennsylvania and western Maryland is dominated by precipitation infiltrating the soil. Water output from the Pennsylvania and Maryland aquifer systems include processes such as runoff, ET from the ground and vegetation, and withdrawals for human uses such as irrigation or domestic water supply. The soil-water balance (SWB) code developed by Westenbroek and others (2010) was used to compute the water budget components for a regional area spanning portions of western and central Maryland and Pennsylvania using readily available gridded datasets. The ET outputs from this model are used in the water budget calculation described in a following section. The SWB code utilizes a modified version of the Thornthwaite-Mather soil-water balance approach at a daily time step (Thornthwaite, 1948; Thornthwaite and Mather, 1957) and incorporates spatially distributed soil, meteorological, and land cover data (Westenbroek and others, 2010).

In the water budget analysis in this study, the groundwater recharge and ET components of the budget were estimated using the USGS SWB 2.0 model (Westenbroek and others, 2018). The Pennsylvania-Maryland (PAMD) SWB model area is considerably larger than the area of interest around Fort Detrick. The size of the entire model domain was designed to include sufficient watersheds with streamflow data to fit the modeled budget components to observed streamflow and base flow. The entire model domain contained 55 watersheds varying from 0.37 to 817 mi² in total area (mean watershed size was 96 mi² and the median was 22 mi²). The model domain considered for the Fort Detrick water balance calculations was 14,960 mi². The PAMD SWB model grid consisted of 614,760 cells arranged in 705 columns by 872 rows, where each cell was a square with sides 250 meters (m) (approximately 820 ft) long.

Water Budget Analysis

Determining the water budget components in the vicinity of Fort Detrick Area B is critical to evaluate key aspects of the previously described CSM. The focus of the study design is to determine whether there is evidence that interbasin groundwater flow is a significant process. Interbasin groundwater flow, sometimes termed “underflow,” occurs when the groundwatershed is not coincident with the surface watershed and can be the result of topographic effects, structural controls, or anthropogenic influences. Because assessing this process is a primary goal, the water budget analysis makes no assumption that the water budget is a closed system. Many water budget analyses assume a closed system to estimate difficult-to-measure water budget components, most commonly ET. However, imposing this assumption would necessarily prevent an evaluation of interbasin groundwater flow. Instead, in this study each individual water budget term is calculated using available data or through modeling, and a water budget residual term is calculated. If the water budget residual term is small relative to the uncertainty of the underlying data, then an additional import or export of water (in other words, interbasin transfer) is not needed to fully describe the hydrologic system. Refer to Kampf and others (2020) for further discussion or water budget analyses that do not assume a closed system, including examples of this approach to identify groundwater interbasin transfer.

A water budget analysis is used to describe the inputs and outputs of water within an area. In general, the components of a water budget for a given amount of time are summarized by equation 6.

On sufficiently long time scales, the changes in storage may be set to zero on the assumption that over the time period of interest the change in storage is negligible when compared to the magnitude of the inputs and outputs. In this study, this assumption is likely not valid for the relatively short (approximately 3 year) monitoring period. We therefore use the hydrologic monitoring network to create an independent estimate of the storage change during the water budget period.

In this analysis, the water budget components defined in equation 7 are evaluated within the two nested topographic watersheds at the Carroll Creek streamgages. The budget components are evaluated within each watershed to compute the water budget residual. The following sections describe the methods for quantifying each component. A schematic representation of the water-budget calculation is shown in figure 10.

where

Uncertainty is an important consideration for evaluating the residual term from a water budget analysis. Though the values for individual water-budget components can be subject to substantial uncertainty, many water budget studies neglect this completely (Healy and others, 2007). The water budget residual includes all the errors associated with the individual water budget components, as well as any unquantified budget components. To assess if a residual indicates the presence of an unquantified budget component, it must be compared to the propagated errors from all of the components in the calculation. In this analysis, we estimate the uncertainty of each water-budget component and combine them into an overall standard deviation error of the residual using addition in quadrature, as described in equation 8. Addition in quadrature is appropriate under the assumption that the errors for the components are normally distributed and independent.

where

If E_tot from equation 8 is greater than or equal to ε_tot from equation 7, then this analysis will consider the water budget residual to be insignificant. An insignificant water budget residual suggests that the hydrologic system can be reasonably conceptualized as a closed system and that there is no evidence to support significant interbasin transfer. Due to the uncertainties inherent in a water budget analysis and the coarse spatial representation of the analysis, small amounts of interbasin transfer cannot be ruled out. However, the water budget analysis would be able to evaluate if large magnitude interbasin transfer is occurring.

Groundwater Recharge

Groundwater recharge is defined as the water that enters the groundwater system vertically through the unsaturated soil zone. In karst systems, preferential vertical flowpaths including sinkholes, swallets, and solutionally enlarged conduits can enable rapid, focused recharge. The heterogeneity of groundwater recharge in karst systems makes this a difficult component of the hydrologic system to estimate. In this study, several methods were used to estimate groundwater recharge at Fort Detrick within the Carroll Creek watersheds, each with independent assumptions and limitations. However, due to the difficulty of representing the heterogeneous karst system using these methods, the estimated groundwater recharge is not used in a budget calculation; instead, the information is provided to evaluate the temporal changes during the water budget period and to place the hydrologic system observed at Fort Detrick into a regional context. The following sections describe the multiple methods used to estimate groundwater recharge.

Oxidation-Reduction

Oxidation-reduction (redox) reaction describe the transfer of electrons and chemical potential energy between ions and is a crucial determinant of metal, nutrient, and microbial activity in groundwater. The oxidized or reduced state observed in a sample may only be a recent event in the geochemical evolution along the flow path. Conservative estimates of redox state are made in this study to guide geochemical characterization of the sample and system. In groundwater, dissolved oxygen (DO) is the primary oxidizer, after which NO₃, manganese, iron, and then SO₄, in succession, begin to donate electrons as the water is in an increasingly reduced state.

In this study, a categorical assignment of redox condition was made for each sample following established criteria (McMahon and Chapelle, 2008; Jurgens and others, 2009). The redox condition of each sample was categorized as “oxic” if concentrations of DO were greater than or equal to 0.5 milligram per liter (mg/L) and iron was less than or equal to 100 micrograms per liter (µg/L), or “anoxic” if concentrations of DO were less than 0.5 mg/L and iron was greater than 100 μg/L (McMahon and Chapelle, 2008; Jurgens and others, 2009). Samples that met the criteria for both oxic and anoxic conditions were categorized as “mixed.” During dissolved gas analysis, anoxic or mixed conditions were used as indicators for the possibility of excess N₂ gas resulting from denitrification; anoxic or mixed conditions were also used as indicators during ¹⁴C geochemical corrections for the possible influence of organic carbonate geochemistry.

Stable Isotopes of Water

Stable isotopes of water are useful indicators of the source of recharging water to groundwater. In addition to information about the precipitation sources, the methods can indicate if substantial evapotranspiration occurred between precipitation and groundwater recharge. Typical sources of postprecipitation evaporation include evaporation from surface water bodies and plant transpiration. Understanding the isotopic composition of precipitation at the study area is essential in stable isotope analysis. For this study, the isotopic composition of precipitation was compiled from Rice and others (1993), Plummer and others (2001), and Good and others (2014) to include observations within 100 mi of the study area. A local meteoric water line (LMWL) was fit by least-squares to the observed precipitation isotope data. The line-conditioned excess (LC-excess; Landwehr and Coplen, 2006) was calculated for groundwater samples to indicate additional evaporative fractionation relative to the LMWL. Negative LC-excess values indicate that postprecipitation evaporation likely occurred, whereas positive LC-excess values indicate an input from distinct recharge source(s), including higher elevations or stronger seasonal influence, with stable isotope compositions different from the assumed regional modern precipitation.

Carbon Isotopes

Carbon isotopes provide information on biologic and abiotic processes within the carbon cycle. Additionally, carbon-14 (¹⁴C) is a useful tracer for old (hundreds to thousands of years old) groundwater. Use of ¹⁴C for groundwater age dating often requires geochemical correction for nonatmospheric dissolved inorganic carbon (DIC) contributions. Analytical solutions describing the combined effects of differing DIC isotope sources on the ¹⁴C and δ¹³C of groundwater provide one method of estimating the nonatmospheric contributions. The graphical method of Han and others (2012) facilitates estimation of the isotopic composition of the various carbon exchange reservoirs, grouping of groundwater samples, and identification of dominant reactions and reaction pathways through carbon isotopic space. In combination with observed carbonate lithology and groundwater chemistry, the graphical method indicated the Revised Fontes and Garnier (RFG; Han and Plummer, 2013) analytical correction model is appropriate for most samples in accounting for DIC from soil zone gas and saturated aquifer carbonates. Carbon-14 and δ¹³C values for the soil zone gas and the aquifer carbonates were estimated from the plot and further refined by checking corrected final ¹⁴C values against other tracers. Sites appropriately captured by the RFG model were corrected using the open-system model, which accounts for continued gas exchange between the atmosphere and groundwater.

Dissolved Gases

Dissolved gases were used in this study to estimate groundwater recharge conditions and dissolved gas exchange/production along groundwater flow paths. Dissolved gases were collected with two methods and analyzed at separate USGS laboratories.

Copper tube samples were analyzed at the Denver Noble Gas Laboratory for isotopes of He, Ne, Ar, Kr, and Xe. Dissolved gases O₂, CO₂, N₂, and CH₄ were analyzed at the USGS Reston Groundwater Dating Laboratory. Measured helium isotopes are of particular interest and reported as the ratio of helium-3 to helium-4 in the sample (R) divided by the ratio of helium-3 to helium-4 in the atmosphere (R_a; 1.38 × 10⁻⁶). By definition, air equilibrated water (AEW) R/R_a values close to 1 indicate an atmospheric source. Values less than 1 indicate a crustal terrigenic source and greater than 1 indicate tritium decay or a mantle source (Solomon and Cook, 2000). Dissolved noble gas concentrations of Ne, Ar, Kr, and Xe were interpreted using the closed-equilibrium (CE) model (Aeschbach-Hertig and others, 2000) to determine noble gas recharge temperature (NGT), a proxy for elevation and seasonal timing of recharge, excess air (EA; dissolved gas at time of recharge, which is related to the entrapped air [Ae] CE model parameter by the equation EA=Ae–F×Ae), and the fractionation factor (F) of the gases during recharge. Inverse CE models were optimized by minimizing the sum of chi-squared statistic (χ²) between measured and modeled gas concentrations. Noble gas samples were modeled both including and excluding N₂ as a fitting parameter in the CE model. Helium isotope mass balance following Solomon and Cook (2000) was used to calculate the helium contribution from the atmosphere, the decay of ³H (³He_trit), and terrigenic sources from radioactive decay (⁴He_terr). The atmospheric contribution of AEW+EA is differentiated from nonatmospheric radiogenic sources based on the helium isotopes. In cases where a large proportion of helium comes from radiogenic sources, the helium balance is too uncertain to accurately estimate the tritiogenic component.

Septum stopper glass bottles were collected and analyzed for Ar, N₂, O₂, CO₂, and CH₄ in accordance with USGS Reston Groundwater Dating Laboratory sampling procedures. No attempt was made to halt biological activity in the copper tube sample after collection, whereas stoppered N₂/Ar bottles were chilled to slow such activity. CE-modeled concentrations of DO, CO₂, and CH₄ were compared to measured concentrations from stoppered bottles to determine the relative change in reactive gases from estimated recharge conditions to sampled groundwater conditions.

Age Tracers

Groundwater age dating is the technique of inferring the time since groundwater recharge using the concentration of specific constituents in the water as a proxy for the age (Cook and Herczeg, 2000). One class of constituents used for age dating are atmospheric tracers. These constituents are typically gasses that have had a varying concentration in the atmosphere (typically as a result of human activity). When precipitation enters the groundwater system, it carries information about the concentration of these atmospheric tracers, which are used to estimate the time of groundwater recharge. In some cases, these atmospheric tracers are radioactive, which can provide additional information on the time since groundwater recharge. Another class of constituents used for age dating are accumulating tracers. These constituents enter groundwater through natural processes as it flows through natural aquifer materials. If the rate of accumulation can be estimated, the time since recharge can be calculated.

To estimate the age (time since recharge) of the water samples collected in this study, samples were analyzed for a suite of age tracers, an approach commonly taken because of the varying time sensitivities of age tracers. Under ideal conditions, multiple tracers can be used together to increase confidence in age estimates. In this study, we evaluate ³H (tritium), ³He (helium-3), ⁴He (helium-4), SF₆ (sulfur hexafluoride), and ¹⁴C (carbon-14). Chlorofluorocarbons (CFCs) are also commonly used in groundwater age dating, but these constituents are reported as widespread contaminants within Fort Detrick Area B and are therefore not suitable for groundwater age interpretation in this study. There are multiple ways to estimate groundwater age; the methods applied in this study reflect the best available methods consistent with the conceptual understanding of the hydrologic system and the available tracers analyzed in this study.

The past concentration of ³H in the atmosphere at the study location determines if ³H may be used as a stand-alone age tracer. Nonatmospheric tritium inputs must be considered at Fort Detrick Area B since laboratory and medical wastes are sources of groundwater contamination. The atmospheric background concentration of ³H at Fort Detrick was approximately 10.4 tritium units (TU) (Michel and others, 2018) in 1950 and would have decayed to 0.22 TU in the year 2019. Because Fort Detrick was established as a research center in 1943, nonatmospheric contributions of ³H to groundwater likely occur near to or after the 1950 cut-off. In this study, we define “modern” water as recharged after 1950 and “pre-modern” water as recharged prior to 1950. The presence of modern water can be identified in a sample by a ³H concentration greater than the decayed background level of 0.22 TU.

Tritium/helium-3 dating (Solomon and Cook, 2000) is based on the ratio of parent (³H) to daughter product (³He_trit) and is not influenced by the variable concentration of ³H in the atmosphere. If the components of the helium budget are reasonably constrained, ³H/³He ages are insensitive to initial ³H concentration and to mixing with groundwater containing no tritium. Comparing the estimated initial ³H (³H₀; sum of ³H and ³He_trit) to the atmospheric ³H concentration for a given ³H/³He_trit estimated recharge validates the helium budget calculations and identifies nonatmospheric ³H. Nonatmospheric ³H affects the parent-daughter ratio; therefore, the ³H/³He_trit ages provide a minimum age estimate when nonatmospheric ³H is present. The ³H/³He_trit method provides an estimate of groundwater age assuming no mixing along or between flow paths, approximating a piston flow system. It is recognized that piston flow conditions are infrequent in natural systems and the resultant age is often referred to as the “apparent” tracer age. Apparent ³H/³He_trit ages are quantitative and provide a good relative measure of age of the modern component, but do not represent the full distribution of flow paths captured at a sample location with ³H-dead (negligible ³H) premodern groundwater being specifically excluded.

As a quantitative estimate of premodern groundwater contributions, ³H can be used to estimate the fraction of modern recharge in the sample (Jasechko, 2016). The model assumes mixing of two well mixed reservoirs, approximating exponential mixing, of modern recharge post-1950 and pre-modern recharge. That representation potentially is not representative of physical hydrology in Area B but provides a simplified estimate for the contribution of premodern recharge in each sample. The ³H modern fraction is calculated using a cut-off year of 1950 and the ³H atmospheric history corrected to the respective sample date. The estimated values of ³H modern fraction for samples with likely ³H contamination from nonatmospheric sources are greater than 1 and are not meaningful.

Terrigenic ⁴He (He_terr) accumulates in groundwater from radioactive decay along flow paths such that it can, in some instances, be used to estimate groundwater age (Solomon and others, 1996; Solomon, 2000). The challenge in a quantitative ⁴He age analysis is identification of any potential discrete sources of helium and determining an appropriate effective He accumulation rate along flow paths (Torgersen and Ivey, 1985). Helium-4 accumulation rates can be estimated through comparison to other age tracers (for example, Plummer and others [2012]), but more often the elevated ⁴He is used to qualitatively identify the presence of old groundwater.

Sulphur hexafluoride (SF₆) is an anthropogenically sourced tracer released to the atmosphere resulting from industrial processes starting in the 1960s, making it an effective tracer for young groundwater. Noble gas closed-equilibrium solubility-modeled recharge conditions were used to convert measured concentrations to air-equilibrium values that can then be related to historical atmospheric concentrations to estimate groundwater age (Busenberg and Plummer, 2000).

Groundwater Age Lumped Parameter Models

Age tracer concentrations were analyzed using TracerLPM (Jurgens and others, 2012), which applies lumped parameter models (LPMs) to understand the source and relative contribution of varying flow paths. In most cases, a groundwater sample is a mixture of flow paths with varying ages, and the age distribution represents the probability of a water parcel with a given estimated age occurring in the sample. Age distributions are estimated by LPMs constrained to represent the physical hydrologic system and optimized to minimize the misfit between multiple measured and modeled tracer concentrations. In groundwater studies, fitted age distributions are interpreted to represent the aquifer dimensions and the extent of flow path mixing. Groundwater age metrics, including the mean age and the susceptibility index (described in a following section), provide a useful summary of the fitted age distributions. The four groundwater age metrics used in this study, along with the conceptualization of groundwater age as a distribution of mixed young and old components, are shown in figure 11. The conceptualization of groundwater age observed at a well as a mixture of younger water that reaches the well more directly through fractures and conduits, and older water that has a less direct or slower flowpath is shown in figure 11A; an example age distribution of this conceptualization is shown in figure 11B. The modern fraction metric used in this study, where water recharged prior to 1950 is considered to be pre-modern and its fractional contribution to the water age distribution is estimated using a binary mixing approach, is shown in figure 11C. This date is the cutoff date because the tritium isotopic system used in this study relies on atmospheric concentrations following this date. The mean LPM age metric used in this study is shown in figure 11D. The age distribution is modelled using the tritium isotopic system, so the mean age of the distribution is not sensitive to pre-1950 water. The apparent tritium-helium age, which represents the age as a very narrow distribution (minimal mixing is assumed), is shown in figure 11E. The susceptibility index used in this study, where the LPM modeled modern water age distribution is compared to a reference distribution (in this study a piston flow age distribution with an age of 1 year), is shown in figure 11F. One metric (apparent tritium-helium age) assumes no mixing and the other two metrics (modern fraction mean age and susceptibility index) incorporate mixing and reflect the distribution of estimated post-1950 ages. A modified version of the USGS program TracerLPM (Jurgens and others, 2009) was used for calibration of LPMs from which these two age dating metrics are derived.

The approach for LPM selection and interpretation was dictated by the utility of the various tracers for age dating, comparison of tracer concentrations to expected atmospheric inputs, and the conceptual hydrogeology. At Fort Detrick Area B, the dispersion model (DM) that describes the output tracer concentrations according to the advection-dispersion equation was selected for LPMs. The fitted dispersion parameter (DP) is directly proportional to the amount of mixing along and between flow paths, with low DP values (DP=0.001) approximating piston flow and high values (DP=3) approximating exponential mixing of flow paths. The limited number of tracers useful for input to LPMs in this study reduce how well the models are constrained and are representative of all flow paths captured in a sample. As a result of tracer selection, LPM-derived age metrics of this study provide information on the young portion of the true age distribution.

Statistical Methods

To assess the geochemical similarity among samples, multivariate statistical methods were applied. Statistical analyses provide a formal measure of the multivariate relation of groundwater chemistry across sample sites. The statistical similarity between sample chemistries were interpreted here to indicate similarity of recharge source(s) and flow-paths. Variables for analysis included select tracer data and the major ion chemistry. Variables were removed prior to analysis based on (1) a combination of correlation coefficients (Spearman’s ρ) greater than 0.7 (p-value <0.05), (2) having similar geochemical signals based on working knowledge of system chemistry and ion sources, and (3) a high number of censored values, in other words, greater than 10 percent of values below the laboratory reporting limit. Censored values for the remaining censored variables (Br, Ba, B, Sr, V, Li, Se, and U) were set to the reporting limit. For some kinds of statistical analyses, reporting limit substitution can be problematic due to sample-by-sample variations in reporting limit. In this analysis the laboratory analytical methods and reporting limits were the same across samples and this substitution creates equal values that do not inject false distinctions between values during clustering. Furthermore, the conservative use of substitution for constituents with 10 percent or fewer censored values likely has a relatively small impact on clustering results. Variable values were standardized prior to analysis to the sample population mean and a standard deviation of one. Hierarchical cluster analysis (HCA) was applied to 30 selected constituents. HCA provides an objective method of evaluating the similarity in groundwater chemistry between sample sites for a suite of chemical constituents. Constituent values were centered and scaled, a distance matrix was computed using Euclidean distance, and Ward’s method of minimum variance (Ward, 1963) was applied to cluster samples with similar geochemical characteristics. The heatmaply package (Galili and others, 2018) in R statistical software (R Core Team, 2018) was used to apply the HCA and visualize the results using these methods. HCA produces a dendrogram from which groups can be defined. In this study we use the HCA-produced dendrogram along with professional judgement to define and interpret the groupings.

Correlations between groundwater chemistry and sample site characteristics (including depth of sample interval and site location) were assessed to determine the spatial relation of groundwater chemistry. Correlation coefficients for Pearson’s r, indicating the deviation of sample values from a linear relation, and Spearman’s ρ, indicating a monotonic relation between the variables, were calculated at the 95-percent confidence level. Only statistically significant correlations are reported with a preference toward Spearman’s ρ if both correlation metrics are meaningful. Correlation coefficients and respective statistical significance (p-values) were determined using cor.test() function in R (R Core Team, 2018).

Boring Log Reanalysis

The 76 available logs recorded by drillers and onsite geologists during the installation of monitoring wells across Area B were interpreted to identify (1) the depth of the soil-rock contact, (2) the number of voids or cavities reported during well drilling, and (3) the interpreted thickness of the epikarst, a conceptualized layer of increased transmissivity resulting from solutionally enlarged fractures and weathering. Although the soil-rock contact is readily identified in the boring logs, the number of cavities/voids and the epikarstic thickness require professional judgement to identify. In some cases, a zone of numerous voids was reported; in this analysis we provide the number of distinct voids or void zones reported even though these may vary in thickness. Where these were not clearly enumerated, we indicated the number with a greater than symbol. The transition from epikarst to the fracture-dominated competent rock is rarely abrupt and the degree of confidence in identifying this layer is dependent to some extent on the level of detail within the drillers’ or geologists’ logs. In some cases, this layer is evidently missing, as there are no reported weathered zones matching the conceptualized features of an epikarst.

The depth to the soil-rock contact for individual reported boring logs generally decreases from west to east across Area B, with this depth commonly reported to be 20 to 40 ft on the west edge of Area B and 5–20 ft in the eastern third of Area B and to the east of Carroll Creek (fig. 12). The greatest depths to the soil-rock contact are reported south of the west Area B boundary, where depths of 75 ft (WVLY-2), 65 ft (WVLY-5), and 51 ft (WVLY-1) are reported. The average depth across all 73 reported soil-rock contact depths is 24.5 ft.

Of the 75 boring logs from which the number of distinct voids could be evaluated, 41 (54.6 percent) had at least one void reported during drilling and 14 (18.6 percent) had three or more distinct voids or void zones. The greatest number of distinct voids were reported along the southwest boundary of Area B, with three wells (BMW58D, BMW59D, and BMW47D) having six or more distinct, solutionally enlarged cavities large enough to be detected during drilling. This region also contains a greater number of voids encountered during drilling (fig. 13). The majority of wells drilled in the western third of Area B encountered at least one void large enough to be noticed during drilling. Voids are also encountered elsewhere across Area B, including in the north, which is underlain by the New Oxford Formation, and to the east of Carroll Creek.

Of the 60 boring logs from which the depth to the bottom of the epikarst could be evaluated, 19 (31.6 percent) had no identifiable zones of weathering or no reported voids. No epikarst was identified in these wells and are labeled with “None” in figure 14. The remaining 41 boring logs had an average depth to the bottom of the epikarstic zone of 86 ft below ground surface, with a median value of 59 ft. Six boring logs had voids and weathering to a depth greater than 100 ft below ground surface. One boring, BMW67E, had a reported water-bearing void at 326 ft below ground surface. The thickness to the bottom of the epikarst is variable and, although depths exceeding 75 ft below ground surface were most commonly observed in wells along or south of the southwest Area B boundary, the available data do not have readily apparent trends. This likely reflects the subsurface heterogeneity, the variable detail in reported site boring logs, the varying depths of the logs, and the challenge in interpreting this feature. However, the presence of sites throughout Area B that have limited reported weathering and solution-enlarged cavities suggests that the epikarst is not a laterally continuous layer and that is observed in the boring logs. This interpretation agrees with the description in Arcadis Inc. (2014) of the epikarstic layer as laterally discontinuous.

Geophysical Log Reanalysis

The orientation of fractures and linear features interpreted from geophysical log data reported in appendix D of Arcadis Inc. (2014) are shown in figure 15. Individual logs were sometimes collected at the same location but at different depth intervals, with the relative depth interval of the data collection denoted by letters A–E following the location name. Each point on figure 15 represents the projection of a point orthogonally intersecting a plane, as projected into the lower half of a sphere. The geophysical logs are grouped by lithology and ordered alphanumerically within each group.

The bedding orientation interpretations generally agree with what is presented in Arcadis Inc. (2014), such that the New Oxford Formation units in the northern portion of the site dip gently to the west (BMW29, BMW71A, BMW71BC, BMW72, BMW75, BMW76), whereas the bedding direction in the Frederick Formation is less consistent. Some geophysical logs (BMW70BC, BMW69A, and BMW53BC) show the limestone of the Frederick Formation dipping to the southeast, consistent with regional mapping of the Cambrian units (fig. 2). Some logs (BMW68AB, BMW78BC [New], BMW78BC [Old]) show bedding dipping to the north, and yet other logs show no clear directionality to the bedding (BMW68CD, BMW69BC, BMW74, BMW67A, BMW67B, BMW67D). In all logs, fracture and void orientations were also recorded. A dominant orientation of fractures is not observed, except in the case when small fractures tend to follow the bedding orientation (BMW69A, BMW70BC).

Streamflow Monitoring

Streamflow monitoring at the two streamgage locations revealed important dynamics in the total streamflow and stream base flow during the change from dry to historically wet conditions observed during the monitoring period. Total streamflow peaked during a flood event on May 17, 2018, and stream base flow remained elevated from May 2018 until June 2019, when decreased precipitation accumulation led to a decline in streamflow through the late summer and early fall of 2019 (fig. 16C). In 2020, a more typical seasonal pattern of higher spring streamflow and lower fall streamflow was present. The downstream Carroll Creek streamgage record of stream base flow began on October 11, 2017. In order to estimate the total stream base flow observed for the full water budget period, the October 2017 stream base flow was estimated using a linear regression model (y=1.01x+0.348; R²=0.984) with the upstream Carroll Creek gage monthly totals.

Base flow at the two streamgage locations revealed the seasonal patterns of groundwater discharge (base flow) and surface runoff (quickflow) to Carroll Creek. The two streamgage locations bound the Carroll Creek stream reach that includes numerous springs. By evaluating the differences in stream base flow and quickflow between the streamgages, the aggregate net changes in groundwater discharge and surface runoff to the stream were evaluated. To allow direct comparison between quickflow and base flow, values were normalized by the watershed area and reported in units of a length per unit watershed area (fig. 17). To compute the watershed difference values, the difference in monthly total base flow or quickflow volume between the two streamgages was normalized by the area shown in yellow in figure 17B. Throughout the monitoring period the area-normalized quickflow for the two streamgages was generally similar in timing and magnitude, although the downstream gage usually exhibited slightly more quickflow per unit area, with the exception of May 2018—a period of intense precipitation and flooding. The downstream watershed encompasses more urbanized area and impervious surface in the vicinity of Area A, a likely cause of greater quickflow per unit area within the watershed difference area.

Both streamgages typically have more monthly mean base flow than quickflow. This is especially clear in late 2018 and early 2019, when stream quickflow declines but base flow remains elevated. Base flow increases once again in the spring of 2020. The downstream gage consistently has more base flow per unit area than the upstream gage. That requires high values of base flow per unit area within the watershed difference area, which encompasses the region with identified springs on the east edge of Area B.

The observed stream temperatures (fig. 16D) exhibited a strong sinusoidal annual variation, with the annual daily mean temperature ranging from 5 to 20 degrees Celsius (°C). In this area, stream temperatures typically reach their maximum value in July and their minimum value in January. The temperatures observed at both streamgages were similar, except for a period around the annual maximum and minimum, when the upstream gage had a greater maximum or lesser minimum. Observed specific conductance, also shown in figure 16D, exhibited a lagged annual sinusoidal variation during 2019 and 2020, reaching its maximum value of 600–700 microsiemens per centimeter (µS/cm) in October and its minimum value of 350–450 µS/cm in May of those years. This pattern was disrupted from the beginning of the monitoring period through early 2019, with no evident annual maximum in the summer months of 2018. The summer of 2018 was a period of elevated stream base flow (fig. 16C). During individual storm events, the specific conductance exhibited short-term variation representing different sources of water input to the stream. Most often, a storm event causes the specific conductance to abruptly decrease as relatively low specific conductance water reached the stream through precipitation and overland runoff. The specific conductance increased as base flow became the predominant water source. During winter months, several storm events exhibited the opposite behavior, with pronounced spikes in specific conductance exceeding 800 µS/cm. These spikes represent the influence of road salt applied to roadways that reached the stream in snowmelt or other runoff periods. For most of the observation period, the specific conductance in the downstream gage exceeded the upstream gage, which is consistent with increases in downstream specific conductance observed in synoptic USGS monitoring shown in appendix 1.

Stream Base-Flow Results

The results of the OHS and CaF base-flow separations are shown in table 4 for the 52 streamgages within the region surrounding and including Fort Detrick Area B that are evaluated in the SWB domain. Only five streamgages have the specific conductance time series required for OHS analysis; the remaining gages are evaluated with the CaF method. The computed alpha values range within the 52 watersheds from 0.91 to 0.99, with a mean of 0.96. A mean of 0.95 was reported by Raffensperger and others (2017) for selected Chesapeake Bay watersheds.

The two streamgages installed for this study on Carroll Creek (USGS station numbers 01642198 and 01642199) were both successfully evaluated with OHS, with Nash-Sutcliffe efficiencies of 0.71 and 0.69, respectively. The results of the OHS evaluation for each streamgage are presented in figures 18 and 19. Values of β are suggested to be related to aquifer material, with values of 0.8 for perennial streams with porous aquifers, 0.5 for ephemeral streams with porous aquifers, and 0.25 for perennial streams with hard rock aquifers (Eckhardt, 2005; Raffensperger and others, 2017). The values calculated for the two Carroll Creek gages (0.737 and 0.763) are consistent with a karstic environment with a relatively porous aquifer material.

A comparison between three different hydrograph separation methods (fig. 20) shows that the characteristics of the Carroll Creek watersheds relative to 52 regional streamgages are consistent regardless of which method is used, even if the values differ slightly among the methods. The results of comparisons between PART, OHS, and CaF for two base flow separation metrics include (1) the annual average (water years 2017–20) base flow per unit area (fig. 20A–C) and (2) the base flow index (fig. 20D–F). During the monitoring period, numerous large storm events resulting from localized intense precipitation for the Carroll Creek gages led to flood events in water years 2018 and 2019 (fig. 16) that were not experienced by all regional gages. However, the results for water year 2020 were generally consistent with 3 water-year values shown in figure 20. For the five streamgages for which OHS could be applied, there is close agreement in the estimated recharge per year, though PART slightly overestimates the base flow per unit area when compared to OHS (fig. 20B) and CaF (fig. 20C). The karst-designated watersheds tend to have a slightly higher base flow per unit area (fig. 20C) than those designated non-karst, though there is considerable overlap. The downstream Carroll Creek gage (01642199) had the highest base flow per unit area of all the gages examined over this 3-year period, regardless of which separation method was used. For water year 2020 alone, this watershed had the third highest base flow per unit area (PART) or the sixth highest base flow per unit area (CaF), depending on the algorithm. The upstream Carroll Creek gage (01642198) had a base flow per unit area more typical of other regional karst watersheds. This result suggests that focused groundwater recharge and discharge processes within the lower watershed compose a greater proportion of the total streamflow than many regional karst streamgages.

The base flow index values computed using OHS and CaF are generally in close agreement for the five USGS gages for which CaF and OHS could be computed (fig. 20D). PART tends to overestimate the base flow index relative to both OHS (fig. 20E) and CaF (fig. 20F). This is somewhat in contrast to the results of Sanford and others (2012) for Virginia watersheds, but more consistent with Chesapeake Bay watersheds (Raffensperger and others, 2017). The karst-designated watersheds have a higher base flow index (0.75+) than non-karst watersheds (0.38–0.75), regardless of which method is used. The two Carroll Creek streamgages have intermediate values between these two groups (fig. 20F), indicating that the flashiness of the system is greater than many other regional karst watersheds. Contributing factors to this relative flashiness in the Carroll Creek watersheds may include geology (those watersheds are only composed of approximately 50 percent carbonate units) and development (the lower watershed area [Area A] is developed and impervious surfaces efficiently route stormwater to the stream).

Groundwater Monitoring

Groundwater levels from the monitoring network were interpreted to evaluate several aspects of the hydrologic system. The timeseries of groundwater levels were used to (1) identify the responsiveness of the groundwater system, (2) examine changes in potentiometric gradients during varying hydrologic conditions, (3) estimate a specific yield required for the water budget calculations reported in a subsequent section, and (4) estimate recharge using the EMR method.

Precipitation Monitoring

Monthly precipitation totals observed at the upstream Carroll Creek streamgage are shown in figure 16B. As previously discussed, historically high precipitation in the spring and summer of 2018 was observed. For the water budget analysis, gridded precipitation estimates spanning the watersheds of interest are needed. In this study, PRISM data are used as an input to the SWB model. To ensure the precipitation estimates, which are developed at a national scale, are suitable for use at the watershed scale, a comparison was made between the gridded estimates and the observed precipitation. This comparison was made at two different time intervals; at a monthly interval and for the full duration for which both measures are available (May 1, 2018–December 31, 2020). The total USGS-reported total precipitation and PRISM data agree to within 3.3 percent (157.1 in and 162.3 in, respectively). The monthly observed and PRISM precipitation sums (fig. 25) are more variable but generally show a good agreement. Overall, there is a root mean squared error of 0.81 in. There is a slight (−3.2 percent) bias where low USGS-observed precipitation amounts are overestimated by PRISM. During wet months (greater than 5 in of observed precipitation), observations and PRISM agree within 15 percent. Because longer periods of record will have less relative uncertainty than shorter periods of record (as random noise cancels; see Winter [1981]), a 5-percent estimate of precipitation uncertainty is used, which is more representative of the difference between measured and PRISM totals over the water budget period than individual monthly totals.

Fluorometric Monitoring

Fluorometric monitoring evaluated if quantifiable concentrations of dye are present within the groundwater system from the previous dye tracer tests in 1995 and 2013. Discrete split samples were collected in December 2020 and continuous fluorometric monitoring at key spring sites occurred for a duration of approximately 3 months each.

Comparisons with Regional Streamflow

The SWB model was compared against regional streamflow records to evaluate the suitability of the ET, interception, and net infiltration outputs at the watershed scale and annual time period. Since the Carroll Creek streamgages had a short period of stream gaging (2017–present), a model domain much greater than the Carroll Creek watershed was used to capture streams with a longer record of streamflow data. Of the total monitored watersheds in the region (n=55) the authors selected 34 reference gages that had less than 1,000 mi² of area, weren’t dominated by urban land-use areas, and didn’t have flow control structures such as hydroelectric facilities. A total of 36 watersheds were selected for comparison with the SWB model (34 reference gages and the 2 Carroll Creek gages; fig. 30A). The model was run for the entire domain for the period from January 1, 2000, to December 31, 2020. A simple steady-state budget that did not consider changes of storage was calculated for each gage for each year where streamflow data was available (fig. 30B); observed total streamflow per unit area is shown on the x axis and the net water entering the watershed as precipitation exceeding ET and interception is shown on the y axis. The SWB model in this study tends to slightly underpredict the total streamflow per unit area within the total watershed domain, with a –7.1 percent negative bias. The model consistently underpredicts the downstream Carroll Creek total streamflow, potentially due to local factors not represented in the model. Although there is scatter in the annual total streamflow per unit area, potentially resulting from interannual storage changes, focused groundwater recharge, or other processes not represented by SWB, the fit was deemed acceptable for the purposes of using ET and infiltration component estimates from SWB in the water budget calculation.

The annual results for the SWB infiltration are compared against base flow estimates computed using the CaF method (fig. 30C) and the PART method (fig. 30D). Again, for this evaluation changes in storage are not considered. The SWB model in this study tends to overestimate the annual base flow per unit area calculated with the CaF method, though the model consistently underpredicts the base flow per unit area for the Carroll Creek watersheds. The SWB model does not overestimate the base flow values derived from PART as much (1.1 percent bias), though a greater amount of scatter is observed. Taken together, the fit between observed stream base flow and simulated net infiltration was not deemed to be suitable enough for the net infiltration output of SWB to be used in the water budget calculation.

Carroll Creek Watersheds Results

Two SWB model outputs (ET and interception) and one SWB input (precipitation) are used directly in the water budget calculations in this study; these gridded values are presented in figure 31 for the area including the two watersheds evaluated. The period shown is the total water budget period from October 1, 2017, to December 31, 2020; the annualized values are reported as inches per year.

Sensitivity Analysis

The results of the sensitivity analysis indicate modeled actual ET and net infiltration are most sensitive to soil curve number. Sensitivity analysis was completed on the full model domain, the results for the lower Carroll Creek watershed (which encompasses the upper Carroll Creek watershed) are shown in figure 32. Changes in soil curve number inputs to as much as ±20 percent result in net infiltration output changes of approximately –5 percent or +30 percent, respectively. Soil groups 1, 3, and 4 were found to exhibit the most sensitivity in lower Carroll Creek watershed. Additional sensitivity to root zone is observed but results in less than 5 percent changes in net infiltration output for changes of 20 percent in root zone depth. Very little sensitivity to changes in readily evaporable water or total evaporable water is observed in the three examined SWB outputs.

Actual Evapotranspiration

In the SWB model, the actual ET represents the water evaporated or transpired from the soil. This amount is a function of the available water to be evapotranspired as well as the climate conditions and season, which determine the amount that could potentially be evapotranspired if available. In this study, the actual ET is computed using the SWB model. For comparison, the annual total actual ET estimated in this study was compared to four other sources of actual ET estimates. These products are inclusive of water transpired from the soil as well as evaporated from intercepted vegetation; therefore, the products are most comparable to the sum of the SWB actual ET and interception outputs. The SSEBOP, MOD16A2, and TerraClimate gridded evaporation products are remote-sensing based estimates of ET. A fourth method relies on a simplistic water budget equation (ET = P − Q / A − ΔS, where P is precipitation, Q is total streamflow, A is the watershed area, and ΔS is the change in storage) that neglects interbasin transfer and uncertainty; it is provided as a reference.

A comparison between the estimates of ET is shown in figure 33. The results are presented as the spatial watershed mean for the two study watersheds at the annual to 3-year timescale. In general, the estimates of actual ET at this temporal and spatial scale have a wide range, with 3 water-year mean values ranging from 17.6 to 35.7 in. The difference between the products reflects varying data sources, methodologies, and spatial resolutions. The TerraClimate dataset shows little interannual variability and consistently exceeds the other estimates of actual ET. The remaining products generally show increased actual ET during water year 2018, the wettest year monitored. Conversely, the lowest rates of actual evapotranspiration are reported during water year 2018, the driest year monitored. The P – Q / A – ΔS estimate of ET differs between the two gages, resulting from greater streamflow per unit area in the lower watershed. During all four intervals shown in figure 33, the SWB actual ET and interception sum is within 20 percent of one or both P – Q / A – ΔS estimates.

Interception

Interception, the process by which vegetation surfaces prevent precipitation from reaching the ground surface, is sometimes combined with the overall flux of ET though it is represented as a distinct process in the SWB model. The importance of interception as a water budget term varies with rainfall amount and intensity, as well as the properties of the vegetation, including species and season (Kampf and others, 2020). There are few readily available global or national datasets with a temporal scale that aligns with the study time period. However, Miralles and others (2010) estimated the canopy interception using a multisatellite approach with data from the period of 2003–07. Globally, they report a canopy interception of 19 percent of precipitation for broadleaf deciduous forests. In forested regions of the mid-Atlantic United States, 10 to 20 percent of precipitation is estimated as interception. The interception estimated by the SWB model accounts for approximately 5 percent of the total precipitation during the water budget period. According to the 2016 National Land Cover Database, the upper and lower watersheds are 25.6 percent and 22.3 percent deciduous forest by land area (Dewitz, 2019). If in these forested areas of the watersheds approximately 20 percent of precipitation were intercepted, that would account for an approximate watershed-mean interception value of 5 percent of overall precipitation.

Net Infiltration

In the SWB model, the water that recharges the groundwater system is termed net infiltration. It is a representation of the water estimated to bypass the root zone and approximates the groundwater recharge if unsaturated zone flow and focused groundwater recharge is assumed negligible. In karst systems, this assumption may not be valid. The SWB estimated net infiltration is compared against other estimates of groundwater recharge in figure 34 for the 3 water-year period and for individual years. The SWB net infiltration is compared to (1) stream base flow per unit area, as calculated using OHS, which would equal groundwater recharge if there were no change in storage, no interbasin transfer, and no riparian zone ET; (2) the recharge estimated using the RORA displacement curve method; and (3) the recharge estimated for individual wells using the EMR analysis. The EMR results are only available for water years 2019 and 2020. The EMR method estimates recharge at individual well locations rather than the entire watershed; these wells are all located in the lower watershed area. The median value of estimated recharge from 10 groundwater wells, with the range of reported recharge shown by a vertical black error bar, is shown in figure 34.

Over the 3 water-year period, the recharge estimated from SWB for the lower watershed is consistently less than the stream base flow, recharge estimated using the RORA method, and recharge from EMR analysis. This difference likely reflects the effect of focused groundwater recharge, which is not simulated in SWB but would be characterized by the methods using streamflow (OHS base flow and RORA) and methods using groundwater time series data (EMR). The lower watershed contains a greater proportion of the karstic Rocky Springs Station Member of the Frederick Formation and identified springs and sinkholes than the upper watershed. The SWB recharge equals or exceeds the stream base flow in water year 2018 and water year 2020 and approximately equals the RORA recharge in water year 2020. The closer correspondence between the SWB outputs and other methods for the upper watershed likely reflects less focused groundwater recharge in this watershed, which contains the less karstic Triassic units. Water year 2019 contained a transition from historically wet conditions to more typical hydrologic conditions. Water was released from storage over this period (table 8), which may account for the discrepancy as SWB recharge is simulated to occur instantaneously following precipitation, whereas in the real system stream base flow remains high following a wet period.

Major Ions

The major ion data describe the primary dissolved constituents in groundwater (table 11). For the most part, dissolved ions in groundwater are derived from water-rock interactions with some additional sources from anthropogenic activities.

Oxidation and Reduction

The majority of samples indicate oxic conditions with high dissolved oxygen (DO) and nitrate (NO₃) with the exception of three sites (FR Dd 222 [BMW65B], FR Dd 231 [BMW70A], and FR Dd 239 [AMW 568-13B]) that are anoxic, denitrified (except FR Dd 231 [BMW70A]), and have elevated concentrations of metals (table 11). An oxic categorization does not exclude some reductive consumption of electron donors in the groundwater sample. Chemical reduction is possible if the oxic conditions occurred in a more limited capacity, such that the threshold criteria of Jurgens and others (2009) is not met, or that the groundwater has since become more oxidized along the flow path closer to the sample point after being in a reduced state. Overall, the high concentrations of DO in samples across the study area suggest that chemically reducing conditions do not affect groundwater chemistry along the flow paths intercepted by the majority of samples collected in this study. However, the sample locations in this study were selected to evaluate the conceptual site model and not to characterize contamination source areas.

Tracers

Groundwater tracer chemistry was interpreted to improve understanding of groundwater source and movement across Area B. The tracers discussed in this section include stable isotopes of water, carbon isotopes, dissolved noble and reactive gases, as well as tritium and tritogenic helium-3.

Groundwater Age

Groundwater age was evaluated through multiple techniques using the ³H and ³He_trit results to understand contributing flow paths and time since groundwater recharge at sample sites across Area B. Observed SF₆ samples corrected for excess air exceeded the expected atmospheric inputs, indicating regional SF₆ contamination, likely from industrial sources, and unreliability of SF₆ for groundwater age dating at Area B. As detailed in other sections of the report, this study does not fully describe the carbon isotope evolution and ¹⁴C could not be used for estimating groundwater ages. Limited to ³H and ³He_trit, groundwater age interpretations in this study quantify age metrics for groundwater recharged since approximately 1950 (modern) and the relative proportion of recharge prior to 1950 (premodern). A more complete description of the full distribution of groundwater ages, including quantified ages of the premodern component, was not possible. The metrics used in this study are shown conceptually in figure 11.

Results of the ³H fraction modern analysis indicate the relative contribution of premodern (pre-1950) groundwater in samples. Although the mean age of the premodern component cannot be estimated, the ³H fraction modern provides valuable information on the portion of the age distribution not captured in the ³H/³He ages and LPM results. Results show 16 of the 20 sites are estimated to have equal to or greater than 50 percent modern water (fig. 43, table 14) indicating that although the majority of sampled Area B groundwater is modern (recharged after 1950), substantial fractions of premodern water exist, particularly deeper within the aquifer. The correlation between ³H fraction modern and depth of sample interval (springs set to 0-ft deep) is statistically significant (Pearson’s r=0.48, p-value=0.042), indicating increasing amounts of premodern groundwater with depth. Decreasing values of ³H fraction modern with depth indicate longer flow paths, potentially from fracture flow or diffusion-dominated portions of the aquifer and (or) regional groundwater flow, being captured at the sample sites. Without useful tracers for determining premodern groundwater ages, it is not possible to determine the most likely source of this premodern groundwater with any confidence.

Apparent ³H/³He_trit groundwater ages describe the modern component of the sample age and are insensitive to contributions of ³H-dead premodern water and nonatmospheric ³H. Apparent ³H/³He_trit ages range from about 3 to about 33 years (table 12) at Area B. Elevated He_terr increased uncertainty in the helium balance and precluded calculation of ³H/³He_trit ages for samples from FR Dd 222 (BMW65B) and FR Dd 243 (PW). Calculated ³H₀ values are in relatively good agreement with the historical atmospheric record of ³H in precipitation (fig. 42), providing further validation of helium mass balance calculations and increasing confidence in ³H/³He_trit ages. Apparent ³H/³He_trit ages increase with depth of the sample interval below land surface (Pearson’s r=0.61, p-value=0.007). Ages of modern recharge are oldest along the south boundary of the study area (fig. 44). These sample locations are deeper and located further downgradient in the Carroll Creek watersheds.

Lumped-parameter models (LPMs) estimate age distributions for each sample by describing the mixing behavior of groundwater along and between flow paths. In this study dispersion model (DM) LPMs were fit to ³H and ³He_trit. The selected tracers and model provided the best fit and estimates of the groundwater age distribution for the data. Three measures of the groundwater age were derived from the LPM solutions: mean age, the dispersion parameter, and the susceptibility index.

Lumped parameter model (LPM) mean ages range from about 3 to about 72 years (table 14; fig. 45) and are in general good agreement with apparent ³H/³He_trit ages. Mean ages could not be estimated by LPMs—which rely on atmospheric tracer histories—for three sites with indications of nonatmospheric ³H (FR Dd 239 [AMW 568-13B], FR Dd 241 [BMW5], and FR Dd 242 [BMW5D]). The largest differences in LPM and ³H/³He_trit ages are for samples with young ages, where relative age differences approach method resolution, and for samples with large LPM dispersion parameters consistent with increasing mixing that is not accounted for in the apparent age analysis. The two sample sites with the oldest LPM mean ages are located within the fractured crystalline bedrock (FR Dd 243 [PW]) or within a deep sample interval along the valley margin (FR Dd 222 [BMW65B]), suggestive of a flow system disconnected from the shallow groundwater in Area B. LPM mean ages increase with sample interval depth (Spearman’s ρ=0.66, p-value=0.004) and from east to west (Spearman’s ρ= −0.57, p-value=0.01), supporting increased capture of regional groundwater flow paths at increasing depth.

The LPM-calculated dispersion parameter (DP) ranges from about 0.001, which approximates piston flow and limited mixing, to about 0.9 (table 14). A DP of 0.001 was assumed for two sites (FR Dd 222 [BMW65B] and FR Dd 243 [PW]) for which ³He_trit could not be reliably calculated and the DM was underparameterized. Dispersion parameters are highest for samples collected from the middle of the Area B study area (fig. 46), indicating greater mixing along (longitudinal dispersion) and between (wider range of flow path lengths) flow paths captured at those respective sample sites. In this study, increased dispersion might be related to capture of additional recharge sources from the center of Area B. Conversely, decreased dispersion with depth (fig. 46) suggests more focused recharge and (or) advection dominated flow at greater depths consistent with recharge at sink hole locations and groundwater flow through solution enlarged cavities in the deeper carbonate aquifer. It should be reiterated that estimates of dispersion in this study are based on young groundwater tracers, and as such have limited utility for identifying the most likely premodern groundwater sources.

The susceptibility index (SI) is a comparative measure of the deviation of the full LPM age distribution from a reference piston flow age distribution with mean age of 1 year (in other words, highly susceptible groundwater). As with the other age metric presented in this report, the SI does not account for potential differences in transport properties between the age tracers used to derive the SI and contaminants of concern. Values of SI range from about 0.11 to about 0.61 (table 14), indicating a wide range of relative susceptibility to land surface contamination across Area B. As a combined metric of LPM mean age and DP, values of SI show a similar pattern across Area B (fig. 47) with groundwater being most susceptible to land surface contamination across the center of Area B. Values of SI decrease with depth (Spearman’s ρ = −0.64, p-value=0.005) supportive of local recharge source across Area B and potential for connectivity between local recharge and the deeper groundwater flow system. Sample sites with high values of SI have higher concentrations of NO₃ (Spearman’s ρ=0.48, p-value=0.054), illustrating an observable connection of interpreted susceptibility to land surface activity. Statistically meaningful relations between the SI and nitrate in this study and various land-surface contaminants in another study (Solder and others, 2020) suggest the relative values of SI are useful for understanding susceptibility across Area B.

Hierarchical Cluster Analysis

Hierarchical cluster analysis (HCA) provides an objective method of evaluating the similarity in groundwater chemistry between sample sites for a suite of chemical constituents. The HCA combines the set of 20 samples into seven groups based on geochemical similarity (figs. 48–50; table 15).

Group 1 (table 15) represents shallow groundwater with a substantial component of recharge from within or near Area B. Group 1 is composed of springs and wells BMW35 (FR Dd 229) and BMW35D (FR Dd 230) in the New Oxford Formation (fig. 49). This group is characterized by relatively high NO₃ and PO₄—a possible anthropogenic impact—and δ¹³C, consistent with open-system carbon exchange. The group 1 samples generally have line- conditioned (LC)-excess centered on zero, indicative of similarity to regional precipitation, and Ca:Mg mole ratios that narrowly range from about 2.3 to about 2.4 (fig. 48; table 11). The group 1 samples have a major ion composition that is tightly clustered when plotted on a Piper diagram (fig. 50).

Group 2 samples represent deep groundwater recharged from precipitation from higher elevations and farther from Area B. Group 2 wells generally have deep screened intervals in the Fredrick Formation, with the exception of FR Dd 240 (BMW71C), which is screened in the New Oxford Formation, and chemistry that is characterized by low CO₂ and low ¹⁴C consistent with groundwater interaction with carbonate aquifer units. Group 2 samples have high LC-excess and cooler noble gas temperatures, indicative of higher elevation recharge. The spatial pattern of group 2 samples along the south study area boundary are somewhat suggestive of preferential flow (fig. 49). Group 2 samples generally have a lower ³H relative to the remaining samples, with the exception of FR Dd 222 (BMW65B), indicating older groundwater and lack of connection with shallower groundwater flowpaths with elevated ³H across Area B. On a Piper diagram (fig. 50), group 2 samples have a similar major ion chemistry to group 1, though with increased variability.

Group 3 samples represent a possible mixture of heavily evaporated and high Cl:Br recharge from shallow flowpaths and deeper groundwater from a high-elevation recharge source (group 2). Group 3 wells (1) have a wide range of screened interval depths, from 36 to 231 ft bls; (2) are all located in the Frederick Formation; (3) are tightly grouped in the southeastern part of the study area (fig. 49), and (4) have geochemistry that is characterized by high Cl, Ca:Mg mole ratios, Sr, δ¹⁸O, Ca, and K (table 15). The increasing LC-excess with depth and average-to-warm noble gas temperatures of group 3 samples indicate a mixture of shallow and deeper groundwater flow paths at these wells. Elevated calcium and Ca:Mg ratios in combination with high CO₂, HCO₃, and average pH (table 11) suggest interaction between recent recharge with high CO₂ from unsaturated zone processes and carbonate aquifer units, although very low δ¹³C at two of the wells (FR Dd 231 [BMW70A] and FR Dd 238 [AMW 568-13A]) indicate carbonate chemistry may be impacted by organic carbon sources. Group 3 samples have relatively high sodium and SO₄ (table 11, fig. 48), possibly a result of interaction with more geochemically evolved waters, as represented by samples from FR Dd 222 (BMW65B) and FR Dd 239 (AMW 568-13B), or as a result of the complex carbonate geochemistry previously mentioned. The sodium and SO₄ values result in a grouping distinct from groups 1 and 2 on the Piper diagram, especially for the cation composition (fig. 50).

Group 4 consists of two wells (FR Dd 241 and FR Dd 241), completed in the sandstone of the New Oxford Formation in the northern part of the study area near the active landfill area. Group 4 has the highest ³H values, possibly from an anthropogenic source. It has high magnesium, moderate SO₄, and low NO₃ and Cl concentrations and plots similarly to group 2 samples on the Piper diagram (fig. 50). Finally, three samples are distinct from the other groups (table 15), indicating a unique geochemical composition at those sites relative to the other groundwater samples. FR Dd 243 (PW) has relatively high silica and low HCO₃, indicative of the hydrologic position of the well in the crystalline bedrock upgradient of the carbonate aquifer units. Although it plots similarly to group 3 samples on the Piper diagram (fig. 50), the magnitude of the ion concentrations is much lower.

Both FR Dd 222 (BMW65B) and FR Dd 239 (AMW 568-13B) have relatively high trace metals, low NO₃, and high excess N₂ in the case of FR Dd 239 (AMW 568-13B), likely associated with the anoxic conditions at the two sites (table 11). FR Dd 222 (BMW65B) has relatively high sodium and SO₄, possibly related to the well location and depth along the western basin margin, and low ³H and δ¹⁸O that indicate old groundwater recharged from winter precipitation. Water level monitoring shows this well is slow to recover from pumping and to respond to recharge compared to other wells in the study area. This suggests the well is completed in part of the Frederick Formation with a weak hydraulic connection to the rest of the groundwater system in the study area. The well may be representative of the fracture- and matrix-diffuse flow part of the aquifer rather than conduit-flow dominated parts. FR Dd 239 (AMW 568-13B) has relatively high chloride and Cl:Br mass ratio, likely as a result of land surface contamination, high sodium possibly related to well depth, and relatively low LC-excess indicative of a more evaporated recharge source. The elevated sodium, chloride, and SO₄ in these wells results in a distinct water type (fig. 50).