Hydrogeologic Characterization of Area B, Fort Detrick, Maryland

Scientific Investigations Report 2022-5054
Prepared in cooperation with U.S. Army Environmental Command and U.S. Army Garrison, Fort Detrick
By: , and 



The authors thank the following agencies for their coordination and assistance: Watermark-ECC for providing split-well samples for fluorometric analysis, Arcadis Inc. for logistical coordination regarding site access.

The authors would like to thank the following U.S. Geological Survey colleagues for data collection efforts in this study: Roberto Cruz, Nicole Bellmyer, Brian Banks, Taylor Naglieri, Chance Bowersox, Ellie Foss, Caitlyn Dugan, Shane Mizelle, Andrew Greise, Logan Jeffries, Rich Janke, Cory Wright, and Erica Warner. The authors thank Carole Johnson for assistance in geophysical log interpretation and Steve Westenbroek and Martha Nielson for assistance in applying the soil-water balance model. Mitchell McAdoo and Mark Kozar are thanked for reviews which greatly improved the manuscript.


Groundwater in the karst groundwater system at Area B of Fort Detrick in Frederick County, Maryland, is contaminated with chlorinated solvents from the past disposal of laboratory wastes. In cooperation with U.S. Army Environmental Command and U.S. Army Garrison Fort Detrick, the U.S. Geological Survey performed a 3-year study to refine the conceptual model of groundwater flow in and around Area B of Fort Detrick at the site- to regional-scale. The investigation was designed to review the geologic setting, assess the temporal variability of the hydrologic system, evaluate the potential for interbasin groundwater flow, determine the degree of vertical connectivity of the aquifer, characterize the sources and timing of groundwater recharge, and identify if dyes from previous tracer tests continue to drain from the aquifer. This study established a continuous hydrologic monitoring network of 12 water level gages, 2 streamgages, a precipitation gage, and in situ fluorometric monitoring. A water budget analysis was performed using hydrologic monitoring data and a soil-water balance model constructed for the study. 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. Groundwater and spring samples from 20 locations were collected in a 2019 synoptic geochemical sampling event and analyzed for a suite of analytes that included groundwater age tracer constituents.

The karst groundwater system was found to be highly responsive to hydrologic events, with strong water level and stream base flow responses to individual storm events and a historic wet period in 2017 and 2018. The water budget analysis included historic flooding in May 2018, though more typical hydrologic patterns were observed in 2019 and 2020. During most evaluated intervals, the water budget residual was less than the estimated uncertainty on the residual for the two Carroll Creek watersheds, which suggested no substantial net interbasin flow occurs from these watersheds. The watershed difference area, a region that includes Area B, had a significant negative water budget residual, which may be the result of a net interbasin import of groundwater or the result of focused groundwater recharge not simulated by the soil-water balance model. Geochemical analysis and groundwater age dating reveals shallow groundwater (approximately less than [<] 150 feet deep) appears to be relatively young (approximately <30 years) and to be recharged in the vicinity of Area B. In the deep groundwater sampled in this study (approximately greater than [>] 150 feet deep), older groundwater from a differing recharge source, based on stable isotopes and noble gas analyses, is observed and interpreted to represent less direct connectivity to the surface and increased proportions of water recharged to the north and (or) west of Area B. A clustering analysis to reveal groupings within the suite of geochemical data was used to define seven groups. The groupings generally show that wells in similar depths and lateral aquifer positions generally cluster together, with some exceptions. Although limited by suspended sediments, the in situ fluorometric monitoring at springs did not detect any dye leaving the system above the limit of detection for the method. Dye was only detected above the limit of detection in one well, which was used as an injection well during a previous dye tracer test.

The results of this study support and refine the conceptual site model of groundwater hydrology at Area B. The geologic and geophysical log review in this study agrees with prior assessments of physical controls on groundwater flow. A literature review of mid-Atlantic karst studies identified similar controls reported in these environments. The additional characterization of hydrologic responsiveness in this study suggests that hydrologic conditions and events are important considerations when interpreting potentiometric surfaces and contaminant trends over time and highlights the importance of continuous hydrologic monitoring. There is evidence to suggest that either intense focused groundwater recharge occurs in the vicinity of Area B or net along-valley groundwater interbasin flow from the upper study watershed enters the lower watershed and discharges to Carroll Creek. Geochemical analyses also suggest that water recharged from Catoctin Mountain and the elevated areas to the north and (or) west of the site may be present in the older and deeper Area B groundwater.


The U.S. Geological Survey (USGS), in cooperation with the U.S. Army Environmental Command and U.S. Army Garrison Fort Detrick, collected and interpreted hydrologic information to refine the groundwater conceptual site model in the vicinity of the active U.S. Army installation of Fort Detrick, Frederick, Maryland (fig. 1). The installation consists of three separate areas known as Area A, Area B, and Area C. Areas A and B are the site of groundwater contamination resulting from the past disposal of industrial and laboratory wastes. Administratively, the groundwater contamination is addressed separately for these two areas. The groundwater system at Area B is the focus of this report.

Fort Detrick is located in western Maryland and is divided into two primary areas.
                     Area B located to the left/west of Area A
Figure 1.

Locations of Areas A and B, Fort Detrick, Maryland.

Groundwater Contamination at Area B of Fort Detrick

Groundwater contamination at Area B of Fort Detrick primarily consists of volatile organic compounds (VOCs). A total of 13 VOCs have been detected in groundwater above the U.S. Environmental Protection Agency’s (EPA) maximum contaminant level, including the chlorinated solvents trichloroethane (TCE) and tetrachloroethane (PCE), the chlorinated fluorocarbon CFC-11, benzene, and other compounds (Arcadis Inc., written commun., December 20191 ). Biological and radiological materials, herbicides, and defoliants were also disposed of in Area B. Eight sites within Area B have been identified as former disposal sites that received laboratory and demolition wastes from the 1950s through the 1970s (Arcadis Inc., written commun., December 2019). These sites are primarily within the southwestern portion of Area B, though additional waste disposal sites that are sources of VOCs in groundwater are present in north-central Area B. Laboratory wastes including 55-gallon drums of TCE were disposed of in southwest Area B in a region designated “B-11.” High groundwater concentrations of TCE and PCE observed in the mid-1990s suggested TCE and PCE may be present as dense nonaqueous phase liquids. An interim removal action of the laboratory wastes and associated soils was performed in B-11 from 2001 to 2004, though the scope of the excavation was affected when live biological wastes were encountered (Arcadis Inc., written commun., December 2019). Although approximately 4,000 tons of soil, chemical containers, biological and medical waste, and construction debris were removed during the interim removal action, contaminated soils remained at the base of the excavation pits and the area was capped with impervious cover in 2010. Due to the presence of VOCs within the surficial karst limestone aquifer at Area B, in 2009 Area B Groundwater (site FTD-72) was placed on the National Priorities List with the EPA’s Comprehensive Environmental Response, Compensation, and Liability Act (CERCLA; also known as Superfund) Program (Arcadis Inc., written commun., December 2019).


At the time of publication, data were not available from Arcadis Inc.

Previous Investigations of Area B Groundwater

Numerous environmental investigations at Area B of Fort Detrick have been conducted from the early 1970s through the present. Over the course of the past five decades, more than 160 wells have been drilled and sampled as part of a series of investigations to understand the hydrogeological system well enough to fully characterize the nature and extent of contamination and design a suitable remedial strategy. This section summarizes the key previous investigations to understand the hydrogeologic system underlying Area B. A complete review of investigations to date is summarized by the draft final remedial investigation report (Arcadis Inc., written commun., December 2019).

Groundwater well installation at Area B took place during several phases (Arcadis Inc., written commun., December 2019). From 1974 to 1992, 51 groundwater monitoring wells were installed across Area B and sampled to characterize the presence of hazardous constituents. Limited lithologic and well construction data are available for some of these wells. Additional wells were installed during remedial investigation activities from 1994 to 1997; 16 wells were installed in 1994–95, and 6 wells were installed in 1997. In 1998, two additional wells were installed to expand the groundwater monitoring network. In 2011–12, 29 new monitoring wells were installed at Area B, including 8 nested well pairs. In 2013–14, eight monitoring wells were installed in six borehole locations in Area B and on private property to the southwest of Area B and in between Areas A and B. In addition to bedrock wells, 12 shallow offsite groundwater wells were installed in 2017 at or above the soil-bedrock interface. In 2017–18, 5 replacement wells and 16 new monitoring wells were installed across Area B.

Geophysical data has been collected at Fort Detrick Area B to characterize the hydrogeologic system. Several studies have attempted to characterize a geologic contact that occurs in the center of Area B between units of Cambrian and Triassic age, due to its potential influence on groundwater flow. The first was a seismic refraction study (Llopis and Simms, 1994), which also characterized the depth to bedrock and identified the contact between geologic units within the center of the site. The second was a surface geophysical survey undertaken in 2000 using magnetic, frequency-domain electromagnetics, electrical resistivity, and spontaneous potential methods (Arcadis Inc., written commun., December 2019). A third geophysical survey using an azimuthal square-array direct current resistivity method in the southeast corner of Area B was used to identify the presence of a fault or unconformity using these methods (Tepper and Phelan, 2001). An airborne geophysical survey using magnetic and electromagnetic sensors was completed in 2001 (Oak Ridge National Laboratory, 2001). Borehole geophysical data were collected in 2011–12; these data are discussed in a following section.

Studies to characterize the regional controls on groundwater and surface-water flow and identify preferential flow have also been undertaken. In 2000, the U.S. Army Corps of Engineers Research and Development Center performed lineaments analysis on a predevelopment aerial photograph of the Fort Detrick area from 1937 (U.S. Army Corps of Engineers, 2000). In this analysis, the orientation of linear features was recorded, as these can represent the orientations of fractures, faults, joints, and other linear structural features that can serve as preferential groundwater flowpaths. Approximately 17 percent of the features were field checked to ensure they were not anthropogenic. In 2000–01, the USGS undertook a seepage gain/loss study to identify stretches of Carroll Creek and other streams that considerably gained or lost flow. At streams surrounding Fort Detrick, the stream discharge, specific conductance, water temperature, and pH were collected. The results of the streamflow data collection were published in tabular format in James and others (2002).

In karst terranes, where the flow of groundwater can be rapid and complex because of dissolution features, dye tracing is often used as a tool to identify hydrologic boundaries, the rate of groundwater movement, and features relevant to contaminant transport including those determining transverse spread. Three qualitative dye tracer studies have been conducted at Fort Detrick, two at Area B and one at Area A. In these studies, activated charcoal passive samplers were used as the primary means for dye detection, which are useful for detecting the presence of dye at low levels but are not as effective as sample collection or fluorometric monitoring for determining the actual concentrations of tracer in the water (Taylor and Greene, 2008).

During the first tracer study in 1995, eosine, fluorescein, and rhodamine WT were injected at five locations across the eastern and northern portions of Area B and were monitored at 139 groundwater wells, springs, and surface water locations for 13–17 weeks (Ozark Underground Laboratory, Inc., 1997; White and others, 2015). The dye injection points included four 6-inch (in.) diameter shallow injection wells designed to inject water to the base of the epikarst, a zone of solution weathering that was estimated to include the uppermost 30 feet (ft) of the bedrock. The fifth injection point (called EDIP5) was a sinkhole that serves as a resurgence point for a stream (named “Stream 1”) and received rhodamine WT dye. Each injection point accepted 500 gallons of water at a minimum rate of 1 gallon per minute (gpm). Subsequent monitoring for residual dyes in 2008 at seven locations was conducted using activated charcoal passive samplers, with the goal of evaluating the presence of dyes injected in the 1995 test.

In 2001, a second groundwater tracing study was undertaken at Area A, to the east of Carroll Creek. Pyranine dye was injected at well ARW190-1, a 25-ft-deep well with a screened interval of 19 ft. The dye injection consisted of 10 pounds (lb) of 77-percent pyranine flushed with nearly 1,500 gallons of water at a rate of 3.5 gpm. Monitoring occurred for 139 days at 47 locations, including the injection well, Area A monitoring wells, offsite wells, Carroll Creek, and springs (Ozark Underground Laboratory, Inc., 2002). Activated charcoal passive samplers were used as the primary means to detect the dye, with a secondary reliance on grab samples.

In 2013, the most recent dye tracer study was implemented with the intent of characterizing the deeper flow system and identifying if deeper groundwater discharges to Carroll Creek (Arcadis Inc., 2018). On May 21, 2013, 15 lb of fluorescein were injected into well BMW68E with a screened interval of 313–328 ft below land surface (bls), at a previously identified void. The following day, 23 lb of eosine dye was introduced into well BMW67C, which is screened at a depth of 140–155 ft bls. Samples were collected at 127 locations, including monitoring wells, springs, streams, ponds, and private wells. The monitoring period for the study was 9 months.

In addition to the work described above, in the course of remedial investigation activities from 2010 to 2017, horizontal flow meter surveys at 81 well locations, synoptic groundwater sampling events for potential groundwater contaminants, spring and stream surveys, and potentiometric mapping were performed (Arcadis Inc., written commun., December 2019).

Purpose and Scope

This report supports refinement of the conceptual site model of Area B at Fort Detrick using hydrologic data collected at the site from 2017 to 2020 and the results of prior investigations at the site. The report includes a review of the regional and site hydrogeology, the results of a 2019 synoptic groundwater chemistry sampling event, and a water budget analysis using a groundwater and stream-monitoring network established at the site as part of this study. The results of this study are interpreted to refine the following aspects of the hydrogeologic conceptual site model:

  • Evaluation of structural controls on regional groundwater flow.

  • Characterization of groundwater system responsiveness to hydrologic events ranging from individual storm events to interannual variability.

  • Evaluation of potential for regional interbasin groundwater flow that bypasses Carroll Creek.

  • Characterization of the sources and timing of groundwater recharge.

  • Determination of the degree of vertical connectedness and geochemical similarity between different regions of the aquifer.

  • Evaluation of in situ monitoring methods to quantify residual dye concentrations from prior dye tracer tests.

The aspects of the conceptual site model listed above advance the understanding of the complex karst aquifer underlying the site at a range of temporal (daily to interannual) and spatial (sitewide to regional) scales. Informed decision-making about remedial strategies, the interpretation of contaminant concentration trends, and the assessment of contaminant fate and transport benefit from a conceptual model that characterizes the physical hydrology using multiple independent approaches and that is refined as new data are collected (Interstate Technology and Regulatory Council, 2017). In this study, regional groundwater flow assessments including the evaluation of interbasin groundwater flow, structural controls on groundwater flow, and the sources of groundwater recharge provide information on the ultimate destination of water recharged at Fort Detrick Area B on long timescales. Site-scale groundwater flow assessments including responsiveness to hydrologic events, variability in groundwater gradients, and the degree of vertical connectivity within the aquifer provide information to support the interpretation of contaminant patterns and trends collected in separate monitoring programs at the site.


The following sections describe the geologic and hydrologic setting of Fort Detrick, a description of karst hydrogeologic processes, and the findings of regional karst hydrogeologic studies. The previously described conceptual site model for Fort Detrick Area B is introduced and opportunities for refining this conceptual model are identified.

Geologic Setting

Fort Detrick Area B is located at the west edge of the Fredrick Valley synclinorium and in the Piedmont physiographic province in Maryland, a region characterized by a humid subtropical climate, rolling hills, and well-developed fertile soils. The geologic setting is summarized in figure 2. The geologic units underlying Fort Detrick Area B are shown in figure 2A. The dashed line in A indicates the approximate boundary between the conglomerate of the New Oxford Formation and the Rocky Springs Station Member of the Frederick Formation, which is refined in the study area based on boring logs and field outcrop observations. The physiographic regions in the vicinity of Fort Detrick are shown in figure 2B. The extent of Cambrian and Triassic units within the western Frederick Valley are shown in figure 2C, and figure 2D presents a stereonet with the bedding orientation of select units in the western Frederick Valley shown as poles to the plane. The locations of these bedding observations are shown in figure 2C.

In the north-south trending Frederick Valley, the Cambrian units dip to the east and
                        the Triassic units dip to the west
Figure 2.

Geologic units and structural features in the vicinity of Fort Detrick Area B adapted from Southworth and others (2007), including A, site map of geologic units; B, regional physiographic provinces; C, extent of Cambrian and Triassic units within the western Frederick Valley; D, the orientations of western Frederick Valley bedding; and E, a geologic cross-section adapted from Southworth and others (2007).

To the west of Area B lies the Bull Run Mountain Fault, a regional normal fault that marks the boundary between the Piedmont and Blue Ridge physiographic provinces. Near Area B, the strike of the Bull Run Mountain Fault is 72°. Immediately to the west of Area B are Early Cambrian rocks of the Harpers Formation and Weverton Formation, which were metamorphosed during the late Paleozoic Alleghenian orogeny (Southworth and others, 2007). The carbonate and sedimentary rocks that compose the Frederick Valley synclinorium were deposited in a shallow marine environment in the Late Cambrian (Southworth and others, 2007). Much later, extensional forces created the Bull Run Mountain Fault, which was active throughout the Triassic and Jurassic. This fault created two half-graben basins near Area B: the Gettysburg basin to the north and the Culpeper basin to the south. Arkosic sandstones, claystones, and breccias were deposited in both basins and the units currently dip 20–30° W. (Southworth and others, 2007).

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.

Blue Ridge Units

In the western portion of the study watersheds, the Blue Ridge is underlain by the Lower Cambrian Weverton Formations of the Chilhowee Group. The Weverton Formation is stratigraphically lower and crops out further to the west within the Blue Ridge anticlinorium. Its members (Owens Creek, Maryland Heights, and Buzzard Knob) are all primarily quartzite with minor interbedded siltstone. The mapped Chilhowee Group unit within the study watersheds (fig. 2A) is the Harpers Formation, a foliated metasiltstone and biotite-chlorite-muscovite-quartz phyllite (Southworth and others, 2007).

Frederick Valley Synclinorium Units

The Frederick Formation is composed of thick intervals of limestone and dolostone with minor amounts of shale and sandstone. The unit was divided into three separate members by Reinhardt (1974). These are, in ascending order, the Rocky Springs Station Member, the Adamstown Member, and the Lime Kiln Member. The basal Rocky Springs Station Member is the surficial unit within the southern portion of the Area B boundary. It is a clast-containing, dark-gray, lumpy-bedded dolomitic limestone that was deposited at the base of the continental slope. Southworth and others (2007) describe the upsection intervals of the Rocky Springs member as a flaggy-bedded limestone with dusky-yellow to light-olive-gray, silty, dolomitic partings and laminations. Intervals of medium-dark-gray polymictic breccia (cemented angular clasts from multiple rock types) as thick as 30 ft are present as a consequence of repeated continental slope submarine sides; these intervals grade upsection into planar-bedded arenaceous limestone. Clasts range from sand size to 20 in. in diameter on the west side of the Frederick Valley. Southworth and others (2007) mapped a contiguous interval of grayish-black shale as interbedded within the Rocky Springs Station Member that is as much as 60-ft thick. This unnamed shale unit is shown to crop out on the more steeply dipping east side of the Fredrick Valley synclinorium. In cross-section they show the shale unit to be present at depth of less than (<) 300 ft on the west side of the syncline and the unit does not crop out at Fort Detrick Area B. The bedding of the Rocky Springs Station Member of the Frederick Formation west of the city of Frederick is shown to consistently dip east-southeast at ~30°. The Adamstown and Lime Kiln Members, also limestone units, crop out further to the east of the study area. The bedding information for the Cambrian units between the Potomac River and the north edge of the mapped quadrangle, between the Bull Run Mountain Fault and the central axis of the Frederick Valley syncline are shown in figure 2D; these units dip to the southeast.

Triassic Basin Deposits

The sandstone and conglomerate units of the Triassic New Oxford Formation were deposited within the half-grabens formed by the Bull Run Mountain Fault, and the southern part of the Gettysburg basin is identified as the New Oxford Formation. There have been some differences in how authors have described the geologic units in the Gettysburg and Culpeper basins. Southworth and others (2007), following previous workers, designated the sandstone unit in the Culpeper basin as the Manassas Sandstone. Both the New Oxford Formation and Manassas Sandstone are described as containing an upper reddish-brown, medium-coarse grained sandstone and claystones underlain by a basal carbonate conglomerate to breccia. Southworth and others (2007) group the depositional units in the Culpeper and Gettysburg basins using the nomenclature established for the Culpeper basin. In this report we use the nomenclature for the Gettysburg basin, as that is the Triassic basin that underlies Fort Detrick Area B. A description of the geologic nomenclature after that of Weems and Olsen (1997) and Southworth and others (2007) is shown in table 1. Bedding information for the Triassic units in both basins between the Potomac River and the north edge of the mapped quadrangle, between the Bull Run Mountain Fault and the central axis of the Frederick Valley syncline, is shown in figure 2D; these units dip gently to the west. The difference in orientation between the Triassic units and Cambrian units is also illustrated in cross section in figure 2E.

Table 1.    

Relation of geologic units underlying Fort Detrick, Area B described by Weems and Olsen (1997) and Southworth and others (2007).

Gettysburg basin Culpeper basin Description
New Oxford Formation Manassas Sandstone, Poolesville Member Arkosic sandstone interbedded with siltstone that contains muscovite mica. At the top of the unit there is a rooted mudstone with small, light gray carbonate nodules.
New Oxford Formation, conglomerate Manassas Sandstone, Tuscarora Creek Member Basal conglomerate that overlies the Frederick Formation, containing limestone and dolostone clasts.
Table 1.    Relation of geologic units underlying Fort Detrick, Area B described by Weems and Olsen (1997) and Southworth and others (2007).

Site Hydrologic Setting

Area B is located within the Monocacy River Basin and is bordered by multiple streams and springs (fig. 3). On its east edge, Area B is bordered by Carroll Creek. A small perennial stream is located along the south edge that is referred to as “Stream 2” by prior investigations (Arcadis Inc., 2014). The central portion of Area B contains a local depression with an ephemeral wetland area; during wet intervals, water flows through channels in the tall grass at the site and during dry intervals the channels become dry or stagnant. Previous investigations have named these channels “Stream 3” and “Stream 4.” Standing surface water is also present at times along the west edge of Area B, where water from culverts underneath Kemp Lane enters the site. In previous investigations, water entering the western portion of the site is called Stream 1. Stream 1 terminates in a swallet, where the stream enters the groundwater system directly.

The springs studied are along the southeast edge of Area B. Stream 6 is to the north.
                        Stream 1 enters from the west, then transitions to Stream 3
Figure 3.

Hydrologic features near Fort Detrick, Frederick, Maryland.

Karst Hydrogeology

Karst features such as sinkholes, subsurface conduits, and springs are the result of geologic, geomorphic, anthropogenic, and climatic factors. Within Maryland’s Frederick Valley (in which Area B is located), the presence of karst features is determined by the principal factors of lithology, fracturing, and ancestral drainage patterns (Brezinski, 2004, 200710). The early Paleozoic carbonates (including the Rocky Springs Station Member of the Frederick Formation) and the Triassic units that comprise the Culpeper and Gettysburg basins (including the New Oxford Formation) have similar physiographic and solubility characteristics. However, there is a considerable difference in the predominant orientation of fractures and joints between the Ordovician Grove Formation (within the Frederick Valley synclinorium) and the Triassic Leesburg Member of the Bull Run Formation of Southworth and others (2007), within the Culpepper basin (Brezinski, 2007). The difference reflects the tectonic stresses acting on these units during the formation of these units. The Grove Formation has a predominant joint strike of 288°, nearly perpendicular to the Bull Run Mountain Fault, whereas the Leesburg Member has a predominant joint strike of 87°, subparallel to the fault. Karst features, including sinkholes, often align with the dominant joint directions (Brezinski, 2007). These structural controls are inferred to enhance the formation of solutionally enlarged features by way of preferential routing of water along the predominant joint directions.

Studies of groundwater movement within karst in the region surrounding Fort Detrick provide a framework for understanding groundwater flow. Although outside of Area B, intensive field studies, dye tracing, geospatial analysis, and modeling work has taken place in the Valley and Ridge Province of Virginia, a region underlain by carbonate rocks and containing karstic features that continues into Maryland in the Hagerstown area. The Valley and Ridge Province is composed of different geologic units with a differing geologic history; however, the similarity of climate and geomorphic influences means that these studies are useful analogues for the regions of the Frederick valley underlain by carbonates.

A series of studies have informed groundwater flow in the vicinity of Leetown, West Virginia and within the Hopewell Run watershed area (Kozar and others, 2008), which has been extended to the larger Opequon Creek watershed (Kozar and Weary, 2009; Yager and others, 2013) located in Clarke County, Virginia (Nelms and Moberg, 2010), and draws upon previous assessments of groundwater flow developed in Frederick County, Virginia (Harlow and others, 2005) and in Jefferson and Berkely Counties of West Virginia (Kozar and others, 1991; Shultz and others, 1995). This body of work has contributed to the conceptualization of groundwater flow in these karst environments.

The major findings of these analogous studies indicate that groundwater flow in the Valley and Ridge Province is complex. Groundwater primarily resides within fractures and solutionally enlarged bedding planes. The fractures and conduits that transmit much of the water are a small proportion of the overall aquifer volume. The orientation and size of these features, and therefore the rate and direction of groundwater flow, is strongly controlled by complex geologic structures such as cross-strike faults and overturned or tightly folded bedding. Where bedding planes dip steeply, groundwater may circulate to deeper depths than is typical in karst groundwater systems. Low permeability geologic units can act as barriers to flow, increasing flow through fractures, faults, and other sources of secondary and tertiary porosity. The locations of groundwater discharge to the surface, including springs, are typically associated with structural features or lithologic contacts with low-permeability units. Groundwater recharge occurred primarily as a distributed process across a broad area, except during intense localized rainfall. Under those conditions, focused recharge to sinkholes can occur.

The surficial weathered zone of the aquifer underlying the soil, called epikarst, is an important control on shallow groundwater movement. The epikarst is characterized by a high degree of carbonate dissolution, high-angle joints that allow rapid infiltration into the aquifer, and relatively fast groundwater flow (on the order of tens to hundreds of feet per day) as indicated by tracer tests, and is discussed in a following section. The reported depths of the epikarstic zone are between 30 and 60 ft. Below the epikarst is a zone of fractured carbonate bedrock (referred to as the “intermediate zone” by Kozar and others [2008]) with decreasing fracture size and density with depth that has a decreased hydraulic conductivity relative to the epikarst. Estimates of the groundwater age in this zone are approximately 15 to 50 years (McCoy and Kozar, 2008). Below a depth of approximately 250 ft, the fracture apertures decrease and aquifer-test data show that the hydraulic conductivity decreases by half from the intermediate zone. Limited information is available on the groundwater age within the deeper portions of karst aquifers in the Valley and Ridge Province of Virginia. In the fractured and folded karst of the Opequon Creek watershed, Yager and others (2013) reported that 20 sampled springs had a median apparent age of 4.4 years and a median of 14 percent of the water in the samples was recharged prior to 1953.

Dye tracer tests conducted within the Valley and Ridge Province of Virginia and Maryland provide a basis of comparison for the properties of the epikarst and for the results of previous tracer tests at Fort Detrick. Dye tracer tests were conducted at the USGS Leetown Science Center in 1979–80 (Kozar and others, 2008); at the Leetown Pesticide Site and Jefferson County Landfill Site (NUS, 1986); in Jefferson County, West Virginia (Kozar and others, 1991); in the Leetown area in 2004–05 (Kozar and others, 2008); and at the Central Chemical site in Hagerstown, Maryland in 2014–15 (Field, 2017). These studies contain some common observations. Dye recovery at springs typically takes place on the order of weeks and detections can take place several miles from the injection point. Estimated groundwater interpreted linear velocities range from (1) 171 ft/day (NUS, 1986); (2) 50–840 ft/day (along strike) to 30–235 ft/day (perpendicular to strike) (Kozar and others, 1991); (3) 12.5–610 ft/day with a median of 50 ft/day (Kozar and others, 2008); and (4) 16–66 ft/day (Field, 2017). Although these linear velocities are much more rapid than other types of aquifers, the rates are much slower than many reported in conduit-dominated karst aquifers where velocities of 1,000–10,000 ft/day are reported (Jones, 1997; Kozar and others, 2008). This has been interpreted to reflect the presence of a slow-flow component connected with the more rapid transport pathways. More rapid interpreted linear velocities and stronger recoveries were generally observed when dye was released in losing stream reaches (Kozar and others, 2008). Heavy precipitation events were observed to mobilize dye that was sequestered in the subsurface in relatively immobile zones. Dye resurgence was typically observed at multiple points. An interpreted linear velocity of 1,640 ft/day was estimated by Doctor and others (2011) for a karst spring in the northern Shenandoah valley, though only 12 percent dye mass recovery after 5 weeks of monitoring was interpreted to reflect a considerable interaction with slow-moving groundwater.

The groundwater budgets developed by previous studies in the region show that the rock type, and therefore the prevalence of solution conduits, has a strong influence on groundwater recharge rates. Groundwater recharge in the Opequon Creek watershed, as estimated from stream base flow records, were found to range from 6.6 to 9.8 inches per year (in/yr) for carbonate rocks and 3.6 to 5.5 in/yr for siliciclastic rocks (Yager and others, 2013). The water budget of the karstic Hopewell Run watershed reported per unit stream base flow as 17.73 in/yr (46 percent of precipitation) and evapotranspiration (ET) as 24.23 in/yr (63 percent of precipitation) for the period of October 2003–September 2005 (Kozar and others, 2008).

Previously Described Conceptual Site Model at Fort Detrick

Based on the previous work to characterize the site, a conceptual site model (CSM) for groundwater at Area B has been described (Arcadis Inc., 2014; Arcadis Inc., written commun., December 2019). The descriptions of the groundwater system contained within these two reports are referred to in the rest of the report as “the CSM” or “the previously described CSM.” This section details the primary features of this CSM concerning Area B hydrology and the results of studies used to form the hydrogeologic descriptions within the CSM.

In the previously described CSM, Area B hydrology is strongly influenced by karst features within the carbonate lithology, which is visible as sinkholes, springs, and in subsurface voids that are interpreted to be solutionally widened. The surface of Area B is described as mantled by a heterogeneous soil overburden layer that is locally as much as 30-ft deep and thins to a few feet deep near Carrol Creek. Seismic refraction surveys show that the depth to the bedrock-soil interface is locally highly variable (Llopis and Simms, 1994). Several attempts to characterize the contact between the limestone of the Frederick Formation and the basal conglomerate of the New Oxford Formation generally confirm the presence of this geologic transition in the center of the site, though the contact is interpreted as nonlinear and complex. The previously described CSM describes the contact as a steeply dipping strike-slip transfer fault, though a lack of wells intersecting this feature leave the physical nature of this contact unconfirmed.

In most of Area B, the water table is described as occurring within the soil overburden layer. Below the soil layer, weathered bedrock containing subsurface voids is present, with a decreasing frequency and size of solution features to a depth of approximately 150 ft. The greatest proportion of solution features are observed at depths of 32–65 ft, with decreasing proportions to 131–161 ft deep, after which very few features are observed (White and others, 2015). However, solution voids more than 300-ft deep are also observed in some locations. The two primary geologic units underlying Area B are the limestone of the Frederick Formation and the basal conglomerate of the New Oxford Formation. The units are described as equally influenced by karst dissolution processes and have similar bulk hydraulic properties, but multiple interpretations are presented as to whether the units are hydraulically connected or the geologic contact acts as a flow barrier. For example, the presence of springs in Carroll Creek at the location where it crosses the contact between the limestone of the Frederick Formation and the basal conglomerate of the New Oxford Formation has been suggested to result from differing hydrologic properties. A steepening of hydrologic gradients close to Carroll Creek and within the conglomerate of the New Oxford Formation are also attributed to a decrease in transmissivity in the conglomerate (Arcadis Inc., written commun., December 2019). Conversely, the results from dye tracing studies, further discussed in the following section, show transport of dye upward and downward across the geologic contact. In Arcadis Inc. (2014), the contact was described as “transparent” from a hydrostratigraphic perspective, with similar bulk hydrographic properties between the two geologic units resulting in no influence on groundwater flow. A third unit (the sandstone or claystone member of the New Oxford Formation) is observed in the northern part of Area B and is not karst-forming or as transmissive as the two other units.

The flow of groundwater within the bedrock, although complex and locally driven by the presence and orientation of fractures and voids, is generally from west to east towards Carroll Creek. The exception is a local hydrologic mound or high point in the northern portion of the site adjacent to an active landfill that is attributed to the local presence of the lower transmissivity of the claystone of the New Oxford Formation (Arcadis Inc., written commun., December 2019). A potentiometric surface map developed using water elevations reported in the draft final Remedial Investigation Report is shown in figure 4 (Arcadis Inc., written commun., December 2019). Groundwater gradients are generally steeper along the west edge of Area B than within central Area B. Groundwater gradients indicate convergence of flow along the Carroll Creek corridor, though interpreted flow paths do not converge precisely at the stream channel. Both upwards and downwards vertical gradients are observed in periodic potentiometric levels, though most reported gradients are upwards. Generally, the western region of Area B is described as an area of focused groundwater recharge as streams flow from the crystalline bedrock of Catoctin Mountain onto valley karstic units. Streams in this area are observed to lose flow to groundwater and at least one swallet is observed.

The gradient gradually descends, meeting at a low area slightly east of Carroll Creek.
                        The exception is a mound in north-central Area B
Figure 4.

Potentiometric surfaces of water elevations reported by Arcadis Inc., written commun., December 2019.

The previously described CSM identifies geologic bedding, fracture and joint sets, and geologic contacts as potential controls on groundwater flow and the development of karst features that act as preferential flow paths. The bedding orientation of the New Oxford Formation at Area B is described as gently westward dipping, whereas the bedding orientation of the Rocky Springs Station Member of the Frederick Formation is more complex and varied, with local dips exceeding 30 degrees, including under the waste disposal areas in the southwest of Area B. The previously described CSM defines a joint orientation perpendicular to the Bull Run Mountain Fault based on a regional geologic analysis conducted by Brezinski (2007). In the vicinity of Area B, 122 linear features were identified through photogeologic analysis, with a predominant northwest-southeast strike direction roughly perpendicular to the Bull Run Mountain Fault (U.S. Army Corps of Engineers, 2000). The presence of numerous springs within a concentrated area along Carroll Creek suggests an association with the contact between the New Oxford Formation and Frederick Formation, possibly as a result of preferential flow along the contact and a difference of hydrologic properties between the formations. The Bull Run Mountain Fault is not considered a preferential flow zone.

A key aspect of this CSM is all groundwater flow from Area B eventually discharges into Carroll Creek; deep recirculation may occur, but all flow converges into the so-called primary discharge zone. This proposed zone, located along Carroll Creek from Area B to the confluence of “Stream 2” shown in figure 3, contains numerous springs that discharge groundwater to Carroll Creek; the springs have been identified in previous investigations and are located both on and off the Area B footprint. There are several springs that reach the surface at a higher elevation than the banks of Carroll Creek, including South Armory Spring, RISP-3, Hospital Spring, and CC003. The documented springs are primarily on the west bank of Carroll Creek, though springs (including Hospital Spring) are also located on the east side. Springs are not reported within the stream along the south boundary of Area B. The location of gaining and losing reaches of Carroll Creek were assessed in the 2000–01 seepage run data reported by James and others (2002). The mapped results of this data are presented in figures 1.11.10 of appendix 1. In general, stream discharge, specific conductance, and pH increase along the stream network leading to the Monocacy River. The increasing specific conductance and pH is consistent with the addition of groundwater having greater contact time with the carbonate rocks in the drainage basin. Strongly gaining conditions are observed along the east edge of Area B, reflecting the input of numerous springs. These conditions are slightly more pronounced in April 2001 than September 2000. Losing conditions are observed at some locations along the south boundary of Area B during both September 2000 and April 2001. Streamflow was lower in September 2000 than April 2001.

The previously described CSM was developed, in part, on the interpretation of the three dye trace studies conducted at Area A and Area B from 1995 to 2013. During the 1995 tracer study, rhodamine, eosine, and fluorescein injected into the shallow borehole injection points were detected during the study period at wells and springs to the east of the dye injection points, indicating the west-to-east movement of groundwater toward Carroll Creek (Ozark Underground Laboratory, Inc., 1997; White and others, 2015). During the test, dye was not detected at monitoring points south of Area B. The travel rates calculated for fluorescein and eosine dye using a straight-line distance and the first arrival after introduction range from 79 to 246 ft/day, with a mean of 151 ft/day. No rhodamine dye was detected during the test, but a fluorescent dye was detected during a subsequent dye test 18 years later in 2013 at eight springs and one well. The fluorescence peak properties of this material were close to, but out of range of, typical rhodamine peaks. These peaks were interpreted to be rhodamine affected by a deaminoalkylation process that is perhaps related to the alteration of the dye in the environment over the 18-year span (White and others, 2015). The fluorescent dye could also be an unknown material with similar fluorescent properties. The detection of rhodamine-like dye was observed east of the rhodamine dye injection point and along Carroll Creek. A summary of the interpreted flowpaths for the 1995 test, with the assumption that the rhodamine-like observations originated from the rhodamine injection point are shown in figure 5. Dye injected in the limestone of the Frederick Formation was detected in the conglomerate of the New Oxford Formation, indicating hydrologic connectivity between those two units.

Injected dye in western Area B and detected along Carroll Creek indicate west to east
                        transport. Overlap is observed between locations
Figure 5.

Summary of injection locations and detections of dye during the 1995 dye tracer test at Area B, Fort Detrick, Maryland (modified from White and others, 2015).

During the 2001 dye tracer study on Area A, no pyranine dye was detected within any of the springs or surface water locations, which was interpreted as a slow transport rate within the aquifer, though no linear velocities were calculated in the study (Ozark Underground Laboratory, Inc., 2002). Dye was detected in a radial pattern from the injection point, with the greatest number of detections south of the injection point. A summary of detections during the 2001 dye tracer test is shown in figure 6.

Injected dye was detected at 11 sites between 4 and 111 days. Dye was not detected
                        at 14 groundwater sites, nor 16 sites surface water sites
Figure 6.

Summary of injection locations and detections of dye during the 2001 dye tracer test at Area A, Fort Detrick, Maryland.

During the 2013 dye tracer study, eosine was detected at more locations (34) than fluorescein (6). All dye was detected to the east of the injection locations, as shown in figure 7. The interpreted linear velocities of first arrival between the injection locations and springs along Carroll Creek were 91–168 ft/day (eosine from BMW67C) and 30 ft/day (fluorescein from BMW68E) (White and others, 2015). The mean linear velocity from injection locations to monitoring wells was 51 ft/day (eosine from BMW67C) and 10.5 ft/day (fluorsecein from BMW68E) (White and others, 2015). The detection of both dyes at springs indicated the presence of upward transport paths and a hydrologic connection between the deep groundwater and Carroll Creek springs. The dye tracer injected at intermediate well BMW67C (eosine) was detected at more monitoring locations and arrived more quickly than dye tracer injected at deep well BMW68E (fluorescein). Dye injected in the limestone of the Frederick Formation was detected in the conglomerate of the New Oxford Formation, indicating hydrologic connectivity.

Fluorescein and eosine dye were injected at in western Area B and detected at groundwater
                        and spring locations in central and eastern Area B
Figure 7.

Summary of injection locations and detections of dye during the 2013 dye tracer test at Area B, Fort Detrick, Maryland (modified from White and others, 2015).

Taken together, these three tracer tests are interpreted in the previously described CSM to indicate that groundwater from Area B discharges exclusively to Carroll Creek, though the karst aquifer has a notable storage component arising from poorly interconnected karst conduits that retains dye mass and results in dye detections over the course of many years following dye injection. The CSM interprets these studies and the body of hydrogeological work at Area B as supportive of “a standard hydrologic basin model in which groundwater and surface water flow are contained within a single coincident drainage basin” (Arcadis Inc., written commun., December 2019).

Opportunities To Improve the CSM

The CSM previously described by Arcadis Inc. (2014) and the draft final remedial investigation report (Arcadis Inc., written commun., December 2019) is based on a synthesis of a considerable body of site investigations spanning from the 1970s through 2017. As discussed in these documents, a robust conceptual site model is developed when possible from multiple independent lines of evidence. When multiple potential interpretations of the hydrogeologic data are present, a weight-of-evidence approach can synthesize independent results into a holistic conceptual model. Conceptual models for complex hydrogeologic systems are typically refined as additional data are collected at a site (Interstate Technology and Regulatory Council, 2017). A review of the previously described CSM identified several opportunities for further refinement through the collection of new data and independent evaluation of previously collected data.

The site geologic setting is complex and has a direct influence on the formation of karst features, the preferential flow of water within the fractured rock, and the location of springs. Although the previously described CSM discusses the local site geology in detail, regional structural influences on the groundwater flow system receive less consideration. Additionally, comparison of the Area B groundwater system to other well-studied karst groundwater systems within the mid-Atlantic will either further support the CSM or identify gaps in understanding. In this study we include an independent review of regional geologic structures, geophysical data, the state of knowledge on their influence on the Area B groundwater system, and provide a comparison with other karst groundwater systems to refine this aspect of the conceptual site model.

Karst systems are highly dynamic and rapidly respond to hydrologic changes. Depending on the conduit density, conduit connectivity, and structural controls within the aquifer, the magnitude and temporal scale of hydrologic responses to hydrologic events can vary. The responses can include changes in flow directions, changes to the locations of surface watershed and groundwatershed boundaries, and activation of different portions of the karstic and epikarstic conduit network (White, 1999). Understanding this variability is the key to characterizing any karst groundwater system. To date, groundwater gradients, water levels, and streamflow at Area B have primarily been characterized using synoptic (point-in-time) observations. These point-in-time observations may be suitable for characterizing seasonal or interannual responses but may not adequately characterize short-lived but important hydrologic responses. In this study, improved high-frequency hydrologic monitoring characterizes the event-driven changes to the aquifer.

One of the primary challenges of characterizing a karst groundwater system is to define the groundwatershed. In some systems this may be quite different from the topographically defined surficial watershed resulting in what is termed interbasin groundwater flow. Understanding the groundwatershed and its relation to the topographic watershed is important to characterize the ultimate destination of water recharged into the aquifer. Due to the orientation of Catoctin Mountain and the Monocacy River, regional groundwater flow from the west to east beneath Area B that bypasses the surficial watershed should be considered. The two prior dye tracer studies implemented at Area B were intended to help define the groundwatershed and determine the destination of water recharged within Area B contaminant source areas. As discussed in the draft final remedial investigation report (Arcadis Inc., written commun., December 2019), the design and interpretations of the study did not rule out a deeper regional groundwater flow path that bypasses Carroll Creek and flows to Monocacy. An independent assessment on the potential for interbasin groundwater flow will provide an additional line of evidence that can support the interpretation of the dye tracer results. In this study, a water budget approach is used to evaluate the potential for interbasin groundwater flow.

The depth of groundwater circulation and the influence of multiple porosities on connectivity within the aquifer underlying Area B has important implications for understanding contaminant transport and developing potential remedial strategies. One of the motivations for the 2013 dye tracer study was to investigate the depth of surface connection within the groundwater system and elucidate the connectedness of the deep and shallow flow within the aquifer by injecting dye into the deep portion of the aquifer. Dye injected deep within the aquifer was subsequently detected at surface springs, indicating a connectivity with the surface. However, observed dye retention at the injection locations also indicates that at depth within the aquifer there is the potential for storage in flowpaths relatively disconnected from the surface at the timescale of interest. Improved understanding of the degree of vertical mixing within the aquifer and identification of the connectivity between different locations within the aquifer using independent methods will support the interpretation of the past dye tracer studies. In this study, a synoptic age dating and geochemical analysis provides independent lines of evidence on the vertical structure, lateral connectivity, and sources of groundwater recharge across the aquifer.

The presence of dyes reported long after their injection informs the understanding of the hydrologic system. Continued evaluation for the presence of dyes from the 1995 and 2013 injection events on a periodic basis is beneficial for interpreting these studies and understanding the storage and release of the dyes. To date, passive sampling and grab sampling techniques have been employed. In situ monitoring methods for dye concentrations are also employed during dye tracer tests (for example, Field [2017]). These methods have an improved ability to characterize the timing of dye resurgence, though they can be limited by interference due to site-specific factors. In this study, we analyzed water samples collected during a routine sampling event for residual dyes, applied in situ monitoring to three spring sites, and performed benchtop laboratory experiments to quantify the detection limit of these methods at Fort Detrick.

Methods of Data Collection

This study included the collection of hydrologic and geochemical data to characterize the hydrogeologic system. A network of sensors was placed in groundwater wells, springs, and Carroll Creek to record the raw data used to evaluate groundwater gradients, calculate a water budget, estimate aquifer storage changes, and evaluate the presence of fluorescent dyes. A single geochemical sampling event in September-October 2019 collected samples from which the age and geochemical characteristics of the water were evaluated.

Hydrologic Monitoring Network

A telemetered hydrologic monitoring network was installed at Area B to monitor the responsiveness of the aquifer to changes in hydrologic conditions and to conduct the water budget analysis. Two USGS streamgages were installed on Carroll Creek and monitored for stream stage, discharge, specific conductance, and water temperature from October 2017 to October 2020. The upstream Carroll Creek streamgage has an area of 4.21 square miles (mi2) (10.90 square kilometers [km2]), and the lower Carroll Creek streamgage has an area of 7.30 mi2 (18.90 km2). A precipitation gage was installed at the upstream gage on the north edge of Area B. Water-level monitoring was conducted at eight monitoring wells from September 2018 to October 2020. An additional four wells were added in December 2019 to increase the extent of the network along the south and west edges of Area B. All water monitoring data were collected using standard USGS techniques and methods and published using the National Water Information System (NWIS; Cunningham and Schalk, 2011). A summary of the hydrologic monitoring conducted for this study and the locations of these monitoring points are shown in figure 8 and table 2.

Table 2.    

Summary of the hydrologic monitoring network.

[Depths are rounded to the nearest foot. ID, identification; USGS, U.S. Geological Survey; NO, New Oxford Formation; NOC, New Oxford Formation conglomerate; FFRS, Frederick Formation Rocky Springs Station Member; HF, Harpers Formation; GC, geochemistry synoptic sample; WL, water level monitoring; DF, discrete fluorometric monitoring; CF, continuous fluorometric monitoring; Q, stream discharge; P, precipitation; SC, specific conductance; T, temperature; —, not applicable]

Fort Detrick ID USGS site name USGS site number Site type Geologic formation Depth, feet below land surface Screened interval, feet below land surface Monitoring conducted
AMW 568-13A FR Dd 238 392548077255101 Well FFRS 40 30–40 GC
AMW 568-13B FR Dd 239 392548077255102 Well FFRS 197 177–197 GC
BMW35 FR Dd 229 392603077263701 Well NOC 97 87–97 GC
BMW35D FR Dd 230 392603077263702 Well NOC 167 155–165 GC
BMW47D FR Dd 246 392604077274201 Well FFRS 102 92–102 GC, WL
BMW5 FR Dd 241 392625077270301 Well NO 35 25–35 GC
BMW53B FR Dd 219 392600077271001 Well FFRS 87 77–87 GC, WL
BMW53F FR Dd 220 392559077271002 Well FFRS 311 295–305 GC, WL, DF
BMW5D FR Dd 242 392625077270302 Well NO 119 109–119 GC
BMW65B FR Dd 222 392608077273602 Well FFRS 313 298–313 GC, WL
BMW70A FR Dd 231 392553077264501 Well FFRS 36 31–36 GC
BMW70C FR Dd 232 392553077264503 Well FFRS 231 221–231 GC
BMW71C FR Dd 240 392605077263001 Well NOC 232 222–232 GC
BMW78A FR Dd 233 392610077270801 Well FFRS 78 68–78 GC
BMW78B FR Dd 226 392611077270801 Well FFRS 150 140–150 GC, WL
CC003 FR Dd 244 392553077262601 Spring NOC 0 GC
County Shallow FR Dd 234 392552077261601 Well FFRS 109 99–109 GC
FR Dd 243 (PW) FR Dd 243 392617077275801 Domestic well HF 200 Open hole GC
RISP-3 FR Dd 213 392556077263301 Spring NOC 0 GC, CF, SC, T
South Armory Spring FR Dd 245 392605077262801 Spring NOC 0 GC, CF, DF, SC, T
BMW20D Well FFRS 138 118–138 DF
BMW24D Well FFRS 175 164–175 DF
BMW32 FR Dd 248 392549077271401 Well FFRS 38 23–38 WL
BMW32D FR Dd 249 392549077271402 Well FFRS 85 75–85 WL
BMW47 FR Dd 247 392604077274202 Well FFRS 72 72–72 WL
BMW56D Well FFRS 91 71–91 DF
BMW58D Well FFRS 130 110–130 DF
BMW59D Well FFRS 183 163–183 DF
BMW65A FR Dd 221 392608077273601 Well FFRS 127 112–127 WL
BMW66B FR Dd 223 392556077263801 Well NOC 150 135–150 WL
BMW66D FR Dd 224 392556077263802 Well NOC 275 270–275 WL
BMW67C Well FFRS 155 140–155 DF
BMW69A Well FFRS 32 22–32 DF
BMW77 FR Dd 225 392625077270901 Well NO 71 56–71 WL
Carroll Creek above Rock Creek at Fredrick, MD 01642198 Stream gage 0 Q, P, SC, T
Carroll Creek near Frederick, MD 01642199 Stream gage 0 Q, SC, T
Hospital Spring FR Dd 178 392552077262201 Spring FFRS 0 CF, DF, SC, T
Table 2.    Summary of the hydrologic monitoring network.
Groundwater monitoring locations span the entirety of Area B, with spring sites to
                        the southeast, wells and streams are around Carroll Creek
Figure 8.

Hydrologic, geochemical, and fluorometric monitoring locations within this study. USGS, U.S. Geological Survey.

Fluorometric Monitoring

Continuous water-quality monitoring for fluorescent dyes at three spring sites using in situ fluorometric monitoring was undertaken to determine if detectable amounts of dye were present from previous dye tracing efforts undertaken in 1995 and 2013. These locations are named FR Dd 213 (RISP-3), FR Dd 178 (Hospital Spring), and FR Dd 245 (South Armory Spring) and their locations are shown in figure 8. These three sites were selected due to the suitability (in other words, sufficient water depth) for in situ monitoring and the prior detection of dyes following the 1995 and 2013 dye tracer tests. At each spring site, a calibrated Turner Scientific C6 submersible fluorimeter was deployed along with a YSI Multiparameter Sonde for an approximately 3-month period to measure concentrations of rhodamine, fluorescein, dissolved oxygen, pH, specific conductivity, temperature, and turbidity. The unit was calibrated prior to each 3-month deployment using an eight-point calibration. The calibration was checked at the end of each 3-month deployment to identify sensor drift. When sensor drift was identified, the values were corrected using standard USGS procedures described by Wagner and others (2006). FR Dd 213 (RISP-3) was monitored from November 6, 2019, to January 30, 2020; FR Dd 178 (Hospital Spring) was monitored from February 14 to May 21, 2020; and FR Dd 245 (South Armory Spring) was monitored from June 18 to October 14, 2020.

Discrete dye observations were also collected as 1-liter split samples from 12 sites within the Fort Detrick quarterly landfill monitoring network December 1–9, 2020. These samples were collected following well purging performed during the periodic sampling for site contaminants. Hand-dipped grab samples from FR Dd 245 (South Armory Spring) and FR Dd 178 (Hospital Spring) were also collected at that time. The concentration of rhodamine, fluorescein, turbidity, and the temperature were recorded using a calibrated Turner Scientific C6 submersible fluorimeter in a laboratory setting. Eosine dye was also injected during previous dye tracer tests. To evaluate the potential influence of the presence of eosine dye on apparent detections of rhodamine and fluorescein, a 5,000 micrograms per liter (µg/L) standard eosine solution was evaluated on the fluorimeter used in this study. Apparent dye concentrations of 480.0 µg/L fluorescein and 400.0 µg/L rhodamine were reported by the instrument, indicating eosine has an interfering effect at high concentrations.

In situ fluorometric monitoring can be affected by the interference of turbidity within the water column, as suspended sediment either blocks or reflects light into the optical sensor on the fluorimeter (for example, Saraceno and others [2017]). To understand the influence of suspended particulates on the observed continuous monitoring results, a series of benchtop experiments were conducted by systematically varying the turbidity of water samples between 0 and 100 Nephelometric Turbidity Units (NTU) using standard solutions, while recording the apparent rhodamine and fluorescein concentrations. This evaluation was repeated in March and October 2020.

2019 Synoptic Groundwater Well and Spring Sampling

In September and October 2019, 20 samples were collected from 16 monitoring wells, 3 springs, and a domestic supply well for a suite of constituents including tritium (3H), dissolved noble gasses, sulfur hexafluoride (SF6), trace metals, major ions, nutrients, stable isotopes of water, and carbon isotopes. The sampling locations are shown in figure 8. Prior to sampling, wells were evaluated using slug tests to ensure wells had a high degree of connectivity with the aquifer. A submersible pump was used to purge at least three well volumes from the well, then monitoring for stability in field parameters including temperature, specific conductance, turbidity, pH, dissolved oxygen, and water level using a YSI Multiparameter Sonde, flow-through cell, and electric water-level tape in accordance with USGS groundwater sampling procedures (U.S. Geological Survey, 2006). In two wells, excessive drawdown in the wells indicated a poor connectivity between the well and the aquifer and alternate well locations were selected. Copper tube pinch clamp samples (Weiss, 1968) were collected by running sample water in-line through refrigeration grade copper tubing. Clear tubing was placed upstream to monitor for bubbles and the copper tubes were clamped under back pressure to limit gas exsolution during sample collection. Water was pulled through the sampling tubing using a peristaltic pump operating at an extremely low setting to avoid bubble creation within the tubing. Springs were sampled by inserting the copper tube directly into the flowing stream orifice as far as practical to ensure water was sampled before reaching the atmosphere. Samples at the domestic well were collected from a point upstream of water treatment systems within the home.

Water analyses were performed at the USGS National Water Quality Laboratory, the USGS Groundwater Dating Laboratory, the USGS Stable Isotope Laboratory, Woods Hole Oceanographic Institute, the USGS Denver Noble Gas Laboratory, and the USGS Menlo Park Tritium Laboratory. The data are available in a companion data release (Goodling, 2023).

Methods of Analysis

The following sections describe the methods used to interpret previously collected geologic and geophysical data, perform a water budget analysis, and interpret the geochemical data collected in this study. These methods were applied to (1) evaluate structural controls on groundwater flow, (2) characterize groundwater system responsiveness to hydrologic events ranging from individual storm events to interannual variability, (3) evaluate the potential for regional interbasin groundwater flow that bypasses Carroll Creek, (4) characterize the sources and timing of groundwater recharge, and (5) determine the degree of vertical connectedness and geochemical similarity between different regions of the aquifer.

Geological Data Interpretation

Previously collected geologic and geophysical data from Area B were evaluated primarily to provide an overview of the presence and orientation of features that control the magnitude and direction of groundwater flow. These features include solutionally enlarged fractures, cavities, and voids as well as the orientation of bedding and unweathered fractures.

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.

Seven geophysical logs were collected in the north or northeast of Area B and six
                           logs were collected in the southwest
Figure 9.

The location of geophysical log data interpreted in this study. Geologic units shown are adapted from Southworth and others (2007) and are shown in figure 2A.

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.

Hydrological Data Analysis

The following sections describe the methods applied to interpret the streamflow, water level, and precipitation data collected in this study at Fort Detrick Area B. The interpretations support a water budget analysis to evaluate the potential for interbasin transfer. The methods applied also characterize the hydrologic system during the observation period.

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

  1. Assumes that stream flow is a mixture of two components: base flow and quickflow;

  2. Assumes that the groundwater system acts as a linear reservoir;

  3. Assumes the SC of quickflow (derived primarily from precipitation) is constant, whereas groundwater SC varies in a prescribed manner through time; and

  4. Provides a measure of model fit by comparing observed and simulated values of SC.

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:

Q B j = 1 β α Q B j 1 + 1 α β Q j 1 α β

α and β

are dimensionless adjustable parameters,


is streamflow in units of length cubed per unit time [L3/T],


is base flow in units of length cubed per unit time [L3/T], and


is an index representing the time step (typically a day).

Base flow at time step j is restricted to Q B j Q j . The parameter β is identical to the parameter BFImax 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 (QB) and quick flow (QS), while values of SC are proposed for the base flow (CB) and quick flow (CS) components. These variables, as well as streamflow (Q), are used to estimate stream SC (Csep) by

C s e p j = C B Q B j β + C s Q s j β Q j  
The root-mean-square error, E, between the estimated stream SC and observed stream SC (Cobs) is calculated as
E β = j = 1 n C o b s j C s e p j β 2 1 / 2
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 (CB) (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

C B j = C B ¯ + C B * s sin 2 π t j t 0 365.25 + C B * c
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 (tj) with an annual period, described by a mean value ( C ¯ B ), amplitudes ( C B * s , C B * c ), and a beginning time/day (t0). In the optimization, the following values are estimated that minimize the root-mean-square error (RMSE; eqs. 13): β, C ¯ B , C B * s , t0, C B * c , and CS.

The second OHS model estimates CB 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 CB 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 mi2 in total area (mean watershed size was 96 mi2 and the median was 22 mi2). The model domain considered for the Fort Detrick water balance calculations was 14,960 mi2. 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.

Data Requirements for the SWB Model

In this model, commonly available gridded datasets that describe the climate, soil hydrologic characteristics, land use, and topography are used to estimate soil water budget components using a modified Thornthwaite-Mather soil-moisture accounting method. Input data for the SWB included tabular and gridded datasets of (1) daily meteorological parameters (precipitation and temperature), (2) land-use and land cover classification, (3) hydrologic soil groups, and (4) estimated available soil-water capacity, all of which are required soil water budget components on a cell-by-cell basis as described in detail in the SWB model documentation (Westenbroek and others, 2010).

Precipitation and Temperature

The SWB model code calculates soil water budget components from meteorological data on a daily time step and compiles this information to a user-defined period (daily, monthly and [or] annually). Meteorological data can be gridded or stacked grid files such as NetCDF, which allow for time series analysis. Two different meteorological data sources, Parameter-elevation Regressions on Independent Slopes (PRISM) (PRISM Climate Group, 2021) and Daymet (Thornton and others, 2020), were initially used and compared. An analysis between these and the USGS precipitation gage indicated that the PRISM model agreed better with observed precipitation data. PRISM outputs were therefore selected for the model simulations. Data were downloaded and run in the SWB model using NetCDF grids. Data from PRISM were supplied as model inputs using its native 4 x 4 km grid cell resolution.

Land-Use Classification

Land cover information is required for the SWB model to calculate recharge, as differing land uses will have different impervious cover and percolation rates for various soil types. Land-use cover within the PAMD SWB model area is primarily forest, agriculture, and urban. Land-use classification inputs to the model are assigned using the 30-m pixel resolution 2016 National Land Cover Dataset (Dewitz, 2019).

Hydrologic Soil Groups

The SWB model requires information on hydrologic properties of soils to determine if there is runoff or infiltration on any given piece of land. The U.S. Department of Agriculture’s Natural Resources Conservation Service (NRCS) has classified soils across the United States based principally on their infiltration capacity and runoff potential (Natural Resources Conservation Service Staff, 2009). The PAMD SWB model used soil hydrologic grouping data provided in the NRCS Soil Survey Geographic Database (SSURGO) dataset. Soil hydrologic data compiled to a 90-m raster resolution by Wieczorek (2014) was used in the model. The soils data was numbered to correspond to various lookup tables that the SWB model used to calculate recharge.

Soils are generally assigned one of seven combinations of letter designations (using A–D) with A having the greatest infiltration rate and D having the lowest infiltration rate and highest runoff potential (Natural Resources Conservation Service Staff, 2009). The data provided for the whole model domain of the soil-water balance had seven different soil types: A, B, C, D, A/D, B/D, and C/D. The dual soil classes exist in areas where a near-surface water table or other feature can delay infiltration. Summary statistics on the model domain indicated that a considerable portion of the cells held dual hydrologic classes. Previous reports using the SWB dealt with the hydrologic soil groups by classifying the dual groups as either the higher or lower of the two. During the comparisons with regional streamflow for this SWB model, the authors determined that the higher soil class provided a better fit. Thus, a soil identified as A/D would be considered an A for modeling purposes. The relation between the NRCS soil classes and the properties used in the SWB model are shown in table 3.

Table 3.    

Relation of soil classes and soil groups used in soil-water balance modeling in this study.

[SWB, soil-water balance; NRCS, Natural Resources Conservation Service. Letters A–D refer to soil classes defined by the Natural Resources Conservation Service (Natural Resources Conservation Service Staff, 2009), with A having the greatest infiltration rate and D having the lowest infiltration rate. Dual soil classes (for example, C/D) exist where a near-surface water table or other feature can delay infiltration]

Designation Category
SWB soil group 1 2 3 4 5 6 7
NRCS soil class B C/D D C B/D A A/D
SWB soil class assigned B C D C B A A
Table 3.    Relation of soil classes and soil groups used in soil-water balance modeling in this study.

Soil data grids typically have missing values. The missing value areas are generally under large bodies of water but can also occur on land. To avoid issues with running the model, no-data grid cells were filled by averaging the soil characteristics around the cells missing data. This was done by passing a rectangular shape over the raster. Any cell without data would then receive the mean value of the surrounding five raster cells. This process was done five times until all no-data values were removed.

Available Water Capacity

The SWB model calculates the available water capacity of individual grid cells by comparing land cover information with soil physical properties. Available water is defined as the water held by soil in between field capacity (water held in soil despite gravity) and the permanent wilting point (the point at which plants cannot remove additional water). The SWB calculates available water for plant consumption/recharge by multiplying the rooting depth by the available soil water, which is dependent on soil texture. Information on the moisture holding capacity of various soil types has been measured in soil surveys and published by the U.S. Department of Agriculture for the conterminous United States. The PAMD SWB model utilized soil hydrologic data characteristics from the U.S. General Soils Map available from SSURGO and compiled by Wieczorek (2014). Missing data in the available water capacity grids were filled using the same methods as the hydrologic soil lookup tables.

Lookup Tables

The SWB model uses user-modified tables of land use, soil, and crop-specific information to assign model cell values to calculate recharge. The land-use lookup tables include information such as rooting depth for different types of vegetation, runoff curve numbers, and maximum possible recharge based on the land-use designation.

The PAMD SWB model’s land-use lookup tables included literature values for the rooting depth of vegetation, which is used to determine the amount of water available for ET. Values of rooting depth were taken from literature and applied across various soil types. Runoff curve values used in SWB lookup tables are empirically derived parameters originally developed by the U.S. Department of Agriculture Soil Conservation Service, which estimated the percentage of rainfall running off of various land uses (Hjelmfelt, 1991).

Spatial and Temporal Resolution

Gridded inputs to the SWB model were available at differing scales, requiring a regridding to produce the SWB outputs. Fine resolution inputs included the 30-m resolution land use data and 90-m resolution soil data; coarse resolution inputs included the 4-km resolution precipitation and temperature datasets. The output resolution of the PAMD model was an intermediate 250-m resolution. SWB version 2.0 performs this data assimilation process by way of nearest-neighbor resampling (Westenbroek and others, 2018). In the case of the PAMD model, the 4-km grid resolution climate data that corresponds to the centroid of the 250-m cell is used. This type of interpolation can produce the illusion of greater output resolution than is supported by the underlying data. In this analysis, we only interpret the outputs at the watershed spatial mean scale. The 10.9 km2 upper Carroll Creek watershed has an area equivalent to a 3.3 × 3.3 km2 and the 18.9 km2 lower Carrol Creek watershed has an area equivalent to a 4.3 × 4.3 km2. To illustrate the output resolution, we present the gridded SWB outputs at the 250-m resolution, though we do not interpret the patterns at the subwatershed scale. SWB model outputs are computed at a daily time step, though the results for individual days may differ appreciably from observed conditions due to misrepresentation in timing from the climate datasets. In this study we limit interpretation for the SWB model outputs to periods aggregated to 1 water year or greater.

Parameter Adjustment

The PAMD model parameters were adjusted by comparing estimated recharge with base flow estimates from reference watersheds in the surrounding region. As the Carroll Creek gaging station was only in operation from 2017 onward, a larger group of gages with longer records was required to compare SWB recharge estimates. Of 55 total available gages, the authors selected 34 that (1) were less than 1,000 mi2 in area, (2) weren’t dominated by urban areas, and (3) lacked flow-regulation features such as hydroelectric generation. The 34 reference watersheds, in addition to the 2 Carroll Creek streamgages, provided information on base flow conditions. After the initial SWB setup, parameters were adjusted to match base flow as best as possible. These decisions included the hydrologic soil groups, noted above, where dual soil classes were determined to follow the higher soil class. In addition, parameters within the lookup tables were adjusted, such as rooting depth, readily evaporable water, and total evaporable water. As with previous SWB reports, the authors found that altering the rooting depth in areas with soil group B had the greatest impact on recharge. The adjustment process involved running the entire SWB model domain, which included the 36 watersheds (inclusive of the two Carroll Creek watersheds), from January 1, 2000, through December 31, 2020. Recharge and ET estimates produced by SWB model were then compared to base flow estimates and parameters were adjusted manually for the best fit.

Sensitivity Analysis

To determine model sensitivity to model parameters, a series of 140 sensitivity runs were completed for the entire model domain. Parameters altered to evaluate sensitivity were contained in the land use and irrigation lookup tables. The parameters in the sensitivity analysis included the runoff curve number, rooting depth, readily evaporable water, and total evaporable water; these parameters were highly sensitive in a recent SWB modeling study (Nielsen and Westenbroek, 2019). The values of the four sensitivity parameters were adjusted sequentially for each soil type to 80, 90, 100, 110, and 120 percent of the original values.

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.

+ ∆
= 0.

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.

P E T I Q c c t o t A     Δ S t o t   = ε t o t  


is the spatial average precipitation within the watershed in units of length per unit time [L/T],


is the spatial average evapotranspiration from the watershed in units of length per unit time [L/T],


is the spatial average interception by vegetation in units of length per unit time [L/T],


is the spatial average groundwater recharge in units of length per unit time [L/T],


is the topographically defined area of the watershed in units of length squared [L2],

Q c c t o t

is the total stream flow from Carroll Creek at the watershed outlet in units of length cubed per unit time [L3/T],


is the spatial average change in storage of water within the watershed in units of length per unit time [L/T], and


is the spatial average residual from the total water budget calculation in units of length per unit time [L/T].

A schematic shows how precipitation infiltrates the groundwater system, is intercepted
                           by tree canopy, or is evapotranspirated by plants
Figure 10.

The water budget components considered within this study.

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.

E 1 2 +   E 2 2 +   E n 2 =   E t o t  


is error associated with the ith water budget component; and


is the error associated with the water budget residual term, which will provide uncertainty on the residual (Kampf and others, 2020).

If Etot 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.


Precipitation represents the input of water to the water budget and was evaluated through a combination of direct observations and published gridded datasets. Precipitation was recorded at the upstream USGS streamgage location using an Ott Pluvio 2 weight volume gage along with a Water Log tipping bucket rain gage as a backup in case of malfunction. Gridded estimates of precipitation in this study are from the PRISM model (PRISM Climate Group, 2021). This model estimates daily precipitation at 4-km grid cells for the conterminous United States and is the official spatial climate dataset of the U.S. Department of Agriculture. In order to perform the water budget calculation, the spatial mean of the precipitation within each watershed was extracted from the gridded precipitation dataset using the R package “exactextractr” version 0.3.0 (Baston, 2021). Using this method, grid cells with partial intersection with the watershed boundary were weighted by the fraction of intersecting cell area. To estimate the uncertainty on the precipitation estimates, the observed precipitation is compared to modeled precipitation for the location of the precipitation gage, adjacent to the upstream Carroll Creek streamgage.


Evapotranspiration (ET) within the water budget is estimated using the SWB model. In SWB, gridded estimates of actual ET are created for given air temperatures, vegetation types, soil classes, and antecedent precipitation events. In the SWB definition, actual ET is the water extracted from storage in the soil reservoir and is distinct from water evaporated from vegetation surfaces, which is classified as interception and described in the following section (Westenbroek and others, 2018).

Independent gridded estimates of ET exist for the study area, though these are derived from remote sensing products and may not capture smaller scale variability or may require bias correction for application in small-scale water budget studies. The operational Simplified Surface Energy Balance (SSEBOP; Senay and others, 2013) is a 1-km resolution estimate of actual ET derived from six hourly weather datasets, the standardized Penman-Monteith equation to estimate the daily short grass reference ET, and a parameterized fraction of ET. The MOD16A2 version 6 ET product is a 500-m resolution estimate of ET using the Penman-Monteith equation and a blend of inputs from daily meteorological data and remotely sensed inputs including vegetation property dynamics, albedo, and land cover (Mu and others, 2011). TerraClimate is a 4-km resolution monthly water balance dataset (Abatzoglou and others, 2018). The TerraClimate dataset first calculates a reference ET using the Penman-Monteith equation, then a modified one-dimensional Thornthwaite-Mather water balance model is constructed to estimate actual ET. Although useful as a reference or first-order approximation, the large spatial scale ET estimates from individual remotely sensed products may not be suitable as a direct input into a small water budget without additional evaluation (Allen and others, 2011).

In many water budget studies, ET is estimated by measuring all other components of the water budget and assuming that interbasin transfer is negligible. Although we do not take that approach in this study, the values of ET from the SWB, SSEBOP, MOD16A2, and TerraClimate datasets are compared to the value of total precipitation from PRISM minus the basin-normalized total streamflow measured at the USGS streamgages and the change in storage. This value is representative of a simple estimate of ET assuming no interbasin transfer and does not account for the uncertainty in streamflow, precipitation, or storage; it is provided as a reference only.

Actual ET is difficult to estimate at an appropriate scale for a water budget analysis and is subject to a high degree of uncertainty (Winter, 1981; Healy and others, 2007). A review by Allen and others (2011) identifies typical errors for field instruments including lysimeters (5–15 percent error), soil water measurements (10–30 percent error), and eddy covariance measurements (15–30 percent error) as well as for remote sensing methods (10–40 percent error).


Interception is the process by which precipitation is captured on the leaves of vegetation and does not reach the ground surface before evaporating. Although some definitions of ET may include this process (Savenije, 2004), the definitions used in the SWB model for ET and interception are separate and therefore both must be accounted for to ensure all budget components are included.

Total Streamflow

The total streamflow in Carroll Creek is reported by the two USGS streamgages, one located at the upstream boundary of Area B and one located within Waterford Park, south of Area B (fig. 8). At each location, a rating curve was developed and a continuous record of streamflow was produced for the period from October 1, 2017, to December 31, 2020.

The uncertainty in total streamflow measurements at USGS streamgages is estimated using a quality code system, where discharge measurements that are rated as excellent, good, or fair are considered to be within 5, 10, and 15 percent, respectively, of the actual daily streamflow 95 percent of the time. Poor records have less accuracy than fair, though the percentage is not estimated. The accuracy of the discharge record is determined by both the stability of the stage-discharge relation established at the site and the accuracy of the observations of stage and measurements of discharge in the field. The downstream Carroll Creek gage (station ID 01642199) discharge for the period of record was rated good, except for discharges above 350 cubic feet per second (ft3/s) and estimated discharges, which were rated as poor. The upstream Carroll Creek gage (station ID 01642198) discharge for the period of record was rated as fair, except for discharges above 200 ft3/s, which were rated as poor.

Changes in Storage

In this study the changes in total water storage are approximated to be equal to the changes in groundwater storage. This approximation is justified by the lack of large water bodies within the study watersheds. The changes to the amount of water stored in the shallow unconfined aquifer system is quantified by examining the changes in groundwater levels within the hydrologic monitoring network. In karstic and fractured rock settings, well hydrograph variations may not represent changes in water level (and groundwater storage). For example, this can occur when fractures and voids are rapidly filled with recharging water while the rock matrix remains unsaturated (Healy and Cook, 2002). Although considerable uncertainty will be present in estimates of groundwater storage using water level data, several techniques in this study are applied to minimize these to the extent possible.

Due to the variability of groundwater responsiveness from the heterogeneous karstic aquifer, this analysis uses water level data from across the monitoring network to reduce the influence of individual wells and characterize the changes in storage across Area B. The change in groundwater storage is calculated by multiplying the change in water level by the specific yield (Healy and Cook, 2002). Stream base flow can be another measure of the changes in groundwater storage. For a closed system with no interbasin transfers, stream base flow discharge will be a function of the groundwater storage; in some cases, this is approximated with a power law relation on the basis of theory and observations (Wittenberg, 1999). However, as a result of the complexity of the karst storage processes at the scale of interest and the potential for interbasin transfers, we do not apply any restriction to the functional form of the discharge-storage relation. Instead, we use a regression approach to identify the relation between storage derived from groundwater levels and stream base flow. A suitable regression equation between these two quantities will enable the groundwater storage to be estimated. The magnitude of the storage change is determined by the specific yield value estimated for the wells in this study (described in a previous section). The following steps are used to estimate changes in groundwater storage from the hydrologic monitoring network:

  1. A mean monthly groundwater level deviation was calculated using the groundwater monitoring network for the longest possible interval. To construct this mean, wells with shorter records (BMW47D, BMW47, BMW32, and BMW32D) and wells affected by pumping (BMW65A) were removed. The long-term mean groundwater level was subtracted from each groundwater-level time series to create a groundwater level deviation. The daily mean deviation was created by averaging across all seven wells. Then the monthly mean deviation was created by averaging the value across all days in a month. The resulting time series characterizes monthly change in the groundwater level across the hydrologic network.

  2. The monthly mean groundwater level deviation was converted to groundwater storage deviation by multiplying by the specific yield determined for the well network. The methodology to determine this specific yield is described in a previous section. A unit conversion from feet to inches was applied.

  3. The monthly mean groundwater storage deviation was related to stream base flow, as determined using the OHS method using regression. Two regression models were created and the base flow time series from both Carroll Creek streamgages were used. The lm() and summary() functions from R version 4.1.0 (R Core Team, 2018) were used to create the regression equations using ordinary least squares, and the package hydroGOF() version 0.4-0 (Zambrano-Bigiarini, 2020) was used to calculate goodness-of-fit statistics.

Using the two regression models, the estimated storage at the beginning and end of the water budget period was estimated. The change in these two values is used as the ∆Stot term in equation 7. The ∆Stot for the downstream Carroll Creek gage is used, as the wells are within this watershed. The change in groundwater storage estimated using the methodology outlined above has several sources of uncertainty. The approximation that the mean groundwater level changes within the observed monitoring network at Fort Detrick represents changes in groundwater throughout the water budget watersheds introduces uncertainty. The monitored wells are primarily in karstic geologic units that do not span the entirety of the study watersheds. However, the log-linear regression with stream base flow is a way to evaluate the representativeness of this assumption. Uncertainty is also introduced with the use of the regression model to extend the mean water level time series. However, this uncertainty is small when compared to the uncertainty introduced by assuming a single specific yield value for the entire watershed. In this analysis, an uncertainty of 50 percent is applied to the specific yield value. Assuming that all other uncertainties are negligible when compared to this uncertainty, the resulting change in storage value also has an uncertainty of 50 percent.

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.

Soil-Water Balance Net Infiltration

The application of the USGS SWB 2.0 model (Westenbroek and others, 2018) produces estimates of the net infiltration of water on a daily time step for the model domain. In this deterministic model, commonly available gridded datasets that describe the climate soil hydrologic characteristics, land use, and topography are used to estimate soil water budget components using a modified Thornthwaite-Mather soil-moisture accounting method. Daily gridded estimates for net infiltration, which is the water that infiltrates past the plant rooting zone and represents the estimated groundwater recharge, are an output of the SWB model created for this study. The output grid cells have a spatial resolution of 250 × 250 m. These gridded groundwater recharge estimates were extracted for the Carroll Creek watershed areas using a spatial average and included in water budget calculations. The Thornthwaite-Mather soil moisture accounting approach applied by SWB for this study neglects the effects of focused groundwater recharge that bypasses the soil zone. In the karstic environment within the Carroll Creek watersheds, where focused groundwater recharge is an important process, the model may overestimate the amount of surface runoff and underestimate the amount of groundwater recharge.

Episodic Master Recession Analysis

Groundwater recharge occurs as both a constant-rate and episodic process. Episodic recharge from precipitation occurs during discrete storm events that cause increases in the static water table elevation and is particularly evident in regions with well-drained soils, small unsaturated zones, and preferential flowpaths including in karst areas. The water-table fluctuation method (for example, Healy and Cook [2002], Nimmo and others [2015]) estimates the amount of water recharged to the groundwater system during discrete events using time series of water table elevations. Water table rises are quantified by calculating the difference between water level peaks and the water table elevation that would have been present without the episodic recharge. In between recharge events, the water table in an unconfined aquifer decreases along a recession curve as water drains to stream base flow and ET occurs. The rate of recession is typically highest during wet periods with high water-table altitudes and lower during drought periods with low water-table altitudes. A master recession curve is the relation between the rate of recession and the water table elevation and may be linear or a piecewise function depending on the underlying structure and function of the aquifer.

In this study, the Episodic Master Recession (EMR) algorithms introduced by Nimmo and Perkins (2018) were applied to time series of water-table elevation recorded at site wells. Two wells (FR Dd 221/BMW65A and FR Dd 222/BMW65B) were excluded from the analysis because the time series lacked evidence of episodic rises and falls related to precipitation events. The wells may be disconnected from the secondary and tertiary porosity within the aquifer that rapidly transmits changes in hydraulic head within the aquifer.

The EMR method requires the use of parameters to accurately select discrete recharge events and build the master recession curve. A manual adjustment process was implemented to identify combinations of parameters that reasonably captured all major recharge events without any false detections of events. The number and elevations of the piecewise splits in the master recession curve were manually determined for each well based on visual inspection of the slope (dH/dt) versus water table elevation (H) graphs.

Specific Yield Estimation

The specific yield is a unitless groundwater aquifer property that describes the volume of water drained by gravity for a given volume of aquifer material; it is typically expressed as a fraction between zero and one or as a percentage. The specific yield is an important aquifer property in this study because it is used to (1) estimate the episodic recharge using the EMR analysis described in the preceding section and (2) estimate changes in groundwater storage as described in a following section. The specific yield in fractured rock and karst settings is a combination of the drainable matrix, fracture, and void specific yields. Fractures and voids generally comprise a low proportion of the total aquifer storage but drain much more rapidly than the aquifer matrix. In settings such as these, the fluctuation of water levels in wells may be more representative of storage changes in fractures and voids than of the matrix (Healy and Cook, 2002). The application of a single time-invariant specific yield in these settings is subject to a high amount of uncertainty and should be evaluated with care (Healy and Cook, 2002).

The water budget method described by Olmsted and Hely (1962) and Risser and others (2005) was used to determine the apparent specific yield for each well. As noted by Healy and Cook (2002), this method is the most widely used technique for estimating the specific yield in fractured rock settings. Risser and others (2005) applied the method to a 2.8 mi2 research watershed in fractured bedrock in the Valley and Ridge Province of Pennsylvania and Olmsted and Hely (1962) applied it to a 287 mi2 watershed in the Piedmont Province of Pennsylvania. Using this method, the apparent specific yield is shown in equation 9:

/ ∆


is the apparent specific yield [unitless],


is stream base flow volume per unit area during a recession period consisting of only groundwater discharge and limited evapotranspiration [inches], and


is the average decline in water-table altitude during the recession period [inches].

The method assumes that the water yielded under gravity drainage (gravity yield) during short periods is equivalent to the specific yield. Though the gravity yield is likely less than the specific yield due to incomplete drainage of the geologic material, the specific yield estimated using this method may be suitable because this complete drainage is unlikely under natural (unstressed) conditions. In this study the specific yield derived from this method is the apparent specific yield. The method also assumes a groundwater system with insignificant interbasin transfer during the recession periods and a linear groundwater system that can be characterized with a single specific yield. Delayed storage effects are not included, so the response time of the aquifer to changes in storage are estimated to occur within 1 day.

Risser and others (2005) report that by using a weighted average of 10 upland wells to calculate ∆h, they derived similar apparent specific yield estimates to what was determined by a separate study that included lysimeter percolate (Gburek and Folmar, 1999). When the calculated apparent specific yield was multiplied by changes in water-table altitude to estimate recharge, the method was found to estimate mean annual recharge comparable to estimates from seven other methods examined by the authors (Risser and others, 2005). The water table fluctuation method tended to underestimate recharge, though the mean annual recharge was statistically significantly different from only one of the seven methods (RORA computer program; Rorabaugh, 1964). However, the authors do caution that the variability of water-level responses in fractured bedrock environments and the errors in estimating specific yield mean that, when possible, estimates of recharge from multiple wells should be used.

To apply this method with the Fort Detrick hydrologic monitoring network, recession periods with minimal ET and precipitation were selected for which streamflow and groundwater level data were available. Evapotranspiration was assumed to be negligible from October–May, as in Risser and others (2005). In this analysis, only the downstream Carroll Creek gage (01642199) was used since the monitoring wells are not within the upstream watershed drainage basin. Wells FR Dd 222 (BMW65B) and FR Dd 223 (BMW65D) were excluded, as muted and delayed water level responses suggest disconnection from surficial recharge pathways and are not suitable for this method. The summed daily mean stream base flow discharge was used to calculate S in equation 9, while the difference in daily mean water level in each well was used as ∆h. Recession periods were identified as periods (1) with assumed negligible ET (between October 1 and May 1), (2) at least 5 days without precipitation at the USGS precipitation gage, (3) where stream base flow (determined by OHS) was approximately equal to total streamflow, and (4) with no irregularities in the groundwater level time series.

A linear model was constructed between ∆h and S for each period of declining water level and was forced through 0, since conceptually in a closed system with no streamflow there would be no change in water level. The R version 4.1.0 functions lm(), confint(), and summary() were used in this analysis. The slope of the line from the linear model represents the specific yield. A minimum threshold of 10 qualifying recession events was considered suitable for creating a regression model. Only regression models where the p-value calculated using the t-statistic from the t-distribution was less than 0.05 were retained (indicating 95-percent confidence that the slope is not zero) and where the adjusted R2 value exceeded 0.66 (indicating greater than 66 percent of the variance is explained by the linear model). The estimated 95-percent confidence interval on the specific yield is also reported for the accepted regression models.

RORA Computer Program

An additional method of estimating groundwater recharge is using the recession-curve displacement method incorporated in the RORA computer program (Rorabaugh, 1964). The method estimates the groundwater recharged for each peak in a streamflow hydrograph and is based on a one-dimensional mathematical equation describing the instantaneous and uniform rise in groundwater hydraulic head throughout the aquifer contributing to a stream. The method assumes a closed system with no interbasin transfers of water. The method represents the net groundwater recharge, which is the total recharge minus any evaporation from the groundwater system through plant roots. Due to the model assumptions regarding the hydrologic system and timing, Rutledge (2000) recommends comparing the results of the recession-curve displacement method to other methods of estimating groundwater recharge and at a timescale greater than 3 months.

To apply the method, a basin recession index, K, describing the time required for stream discharge to decline one log cycle is required. The recession index, α, determined for the Carroll Creek streamgages in the application of OHS using the program RECESS, was also applied in the application of RORA, with a unit conversion as described in Raffensperger and others (2017). Sensitivity tests in Rutledge (2000) indicate a relative insensitivity to variations in the recession constant. In this analysis, the function rora() from the DVStats R package (Lorentz and DeCicco, 2021) was used to perform the analysis. Aggregation of recharge events to the water year were used for comparison with other groundwater recharge estimation methods.

Methods of Geochemical Interpretation

The geochemical data collected at 20 locations during the 2019 synoptic sampling event were interpreted to inform the overall conceptual understanding of the physical hydrologic system. The focus of the geochemical interpretation is identifying spatial patterns and depth relations that may indicate groundwater recharge sources and similarity of flow-path geochemical processes. Groundwater samples were analyzed for major ions and trace metals, stable isotopes of oxygen and hydrogen in water (δ18O, δ2H) for recharge source identification, carbon isotopes of dissolved inorganic carbon (δ13C, 14C) for age dating, dissolved gases (Ar/N2, Ne, Ar, Kr, Xe, 3He, 4He, SF6) for characterizing recharge conditions and age dating, and tritium (3H) for age dating. The geochemical analysis consisted of a combination of statistical and geochemical modeling analyses.


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 NO3, manganese, iron, and then SO4, 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 N2 gas resulting from denitrification; anoxic or mixed conditions were also used as indicators during 14C 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 (14C) is a useful tracer for old (hundreds to thousands of years old) groundwater. Use of 14C 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 14C and δ13C 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 δ13C values for the soil zone gas and the aquifer carbonates were estimated from the plot and further refined by checking corrected final 14C 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 O2, CO2, N2, and CH4 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 (Ra; 1.38 × 10−6). By definition, air equilibrated water (AEW) R/Ra 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 (χ2) between measured and modeled gas concentrations. Noble gas samples were modeled both including and excluding N2 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 3H (3Hetrit), and terrigenic sources from radioactive decay (4Heterr). 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, N2, O2, CO2, and CH4 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 N2/Ar bottles were chilled to slow such activity. CE-modeled concentrations of DO, CO2, and CH4 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 3H (tritium), 3He (helium-3), 4He (helium-4), SF6 (sulfur hexafluoride), and 14C (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 3H in the atmosphere at the study location determines if 3H 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 3H 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 3H 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 3H 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 (3H) to daughter product (3Hetrit) and is not influenced by the variable concentration of 3H in the atmosphere. If the components of the helium budget are reasonably constrained, 3H/3He ages are insensitive to initial 3H concentration and to mixing with groundwater containing no tritium. Comparing the estimated initial 3H (3H0; sum of 3H and 3Hetrit) to the atmospheric 3H concentration for a given 3H/3Hetrit estimated recharge validates the helium budget calculations and identifies nonatmospheric 3H. Nonatmospheric 3H affects the parent-daughter ratio; therefore, the 3H/3Hetrit ages provide a minimum age estimate when nonatmospheric 3H is present. The 3H/3Hetrit 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 3H/3Hetrit 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 3H-dead (negligible 3H) premodern groundwater being specifically excluded.

As a quantitative estimate of premodern groundwater contributions, 3H 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 3H modern fraction is calculated using a cut-off year of 1950 and the 3H atmospheric history corrected to the respective sample date. The estimated values of 3H modern fraction for samples with likely 3H contamination from nonatmospheric sources are greater than 1 and are not meaningful.

Terrigenic 4He (Heterr) 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 4He 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 4He is used to qualitatively identify the presence of old groundwater.

Sulphur hexafluoride (SF6) 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.

Groundwater age is conceptualized as a mixture of older and younger groundwater
Figure 11.

Conceptual diagrams showing the representation of groundwater age used in this study including A, the conceptualization of groundwater age; B, a hypothetical age distribution; C, the 3H fraction modern metric; D, the lumped parameter model mean age metric; E, apparent 3H/3Htrit age metric; and F, susceptibility index metric.

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.

Susceptibility Index

The susceptibility index (SI) provides a quantitative estimate of susceptibility of a well to land surface contamination. More specifically, it is a relative measure of how soon a conservative contaminant in the groundwater recharge area would arrive at the well. The SI is calculated as the normalized difference between the LPM-derived cumulative distribution function (CDF) of age and a reference CDF (Solder and others, 2020) as described by the following equations.

S I = 1 1 + D H
where DH is the Hellinger distance calculated as
D H = i = 1 k P i Q i 2 2


is the cumulative probability of the age distribution of interest for a given time-step,


is the cumulative probability of the reference age distribution for a given time-step,


is the time-step of the cumulative distribution function, and


is the age corresponding to the last time-step.

The reference CDF was defined in this study as a piston-flow model with a mean age of 1 year. The SI is unitless and ranges between 1 (indicating young ages and a narrow age distribution) and approaching 0 (indicating older ages and a broad age distribution). Implicit in the calculated SI are the limitations of the LPMs in this study.

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).


The following sections describe the observed geologic, hydrologic, and geochemical data collected at Fort Detrick Area B. These data were interpreted to characterize the potential controls on groundwater flow, characterize the responsiveness of the coupled groundwater and surface water system, evaluate the potential for interbasin transfer by conducting a water budget analysis, and inform the characteristics of the aquifer underlying Area B through a synoptic groundwater age dating and geochemical analysis.

Geological Data Interpretation

The geological and geophysical data interpreted in this study were collected in previous investigations at Fort Detrick Area B (Arcadis Inc., 2014; Arcadis Inc., written commun., December 2019). The data were interpreted in this report to identify the potential geological controls acting on the groundwater system and to characterize the potential controls that are not well described by the available data.

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.

Reported soil thicknesses generally decrease from west to east across Area B, from
                           values of 30–50 feet to values of 5 to 15 feet
Figure 12.

Reported soil thickness in drilling logs, Fort Detrick, Frederick, Maryland. Boring log data reported by Arcadis Inc., written commun., December 2019. Ft, feet.

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.

The number of reported voids is most commonly zero across Area B, though high numbers
                           (three or greater) are reported in southwestern Area B
Figure 13.

Voids reported during drilling, Fort Detrick, Frederick, Maryland. Boring log data reported by Arcadis Inc., written commun., December 2019.

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.

The depths collected across Area B are generally greatest in southwestern Area B,
                           with values exceeding 100 feet reported at six locations
Figure 14.

Interpreted depth to bottom of epikarstic zone, Fort Detrick, Frederick, Maryland. Boring log data reported by Arcadis Inc., written commun., December 2019. Ft, feet.

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.

Nine logs indicate a gentle northeastern dip. Thirteen have no clear pattern. Three
                           of the 13 have deep eastern bedding dips
Figure 15.

Stereonets indicating the orientation of planar features identified in the geophysical log data.

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).

Hydrologic Monitoring Results

The hydrologic data collected from the monitoring network of streamgages, wells, and precipitation characterizes the responsiveness of the system to individual storm events, seasonal variation, and interannual variability. The data collected in this study for the period from October 1, 2017, to December 31, 2020, are shown in figure 16, including the Palmer Drought Severity Index (Palmer, 1968), as reported by the GridMET dataset (fig. 16A; Abatzoglou, 2013). This metric is derived from estimates of precipitation and temperature data to estimate the cumulative deviation from normal hydrologic conditions. Negative values indicate drier than normal conditions and positive values indicate wetter conditions, with the magnitude indicating the severity of the deviation. The monthly total precipitation recorded at the USGS gage adjacent to the upstream streamgage on Fort Detrick Area B is shown in figure 16B. During the monitoring period at Fort Detrick, an abrupt transition from slightly drier than normal conditions to very wet conditions was driven by historically high precipitation beginning in May 2018. Through the rest of the monitoring period, conditions remained wetter than normal, but the severity of the departure began to decrease beginning in mid-2019. The daily mean total streamflow and stream base flow for the two Carroll Creek USGS streamgages are shown in figure 16C; note the logarithmic y axis. Base flow is calculated using the OHS method (Raffensperger and others, 2017). The daily mean stream temperature and specific conductance observed at the two Carroll Creek USGS streamgages are shown in figure 16D. The daily mean water table elevations observed at 12 groundwater monitoring wells on Fort Detrick Area B are shown in figure 16E. The monitoring period for the wells began in October 2018.

From October 1, 2017, to May 1, 2018 saw monthly precipitation totals increase, then
                        decline from January 1, 2019, to December 31, 2020
Figure 16.

Hydrologic data collected at Fort Detrick Area B, Frederick, Maryland including the A, Palmer Drought Severity Index; B, monthly observed total precipitation; C, daily mean discharge; D, daily mean stream temperature and specific conductance; and E, daily mean water table elevations. A logarithmic axis is used in panel C and a dual y-axis is used in panel D. USGS, U.S. Geological Survey.

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; R2=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.

Stream baseflow per unit area and stream quickflow per unit area sharply increased
                           in May 2018 but declined by the end of 2018
Figure 17.

Monthly mean stream hydrograph separation results for the upstream and downstream Carroll Creek streamgages. USGS, U.S. Geological Survey.

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.

Table 4.    

Results of optimal hydrograph separation (OHS) and Collischonn and Fan (CaF) for accepted models.

[USGS, U.S. Geological Survey; MM, two-digit month; DD, two-digit day, YYYY, four-digit year; SC, specific conductance; SCfit, peak-fitting base flow OHS SC model; OHS Beta, optimized value of beta using optimal hydrograph separation; CaF Beta, beta determined using the method of Collischonn and Fan (2013); CS, quickflow SC; RMSE, root-mean-square error of modeled and observed SC in OHS model; μS/cm, microseimens per centimeter; ft3/s, cubic feet per second; NS, Nash-Sutcliffe efficiency of modeled and observed SC in OHS model; —, no value]

USGS station ID Number of days with discharge data Number of days with SC data Beginning date (MM/DD/YYYY) End date (MM/DD/YYYY) Model Computed alpha (α) OHS beta (β) CaF beta (β) CS RMSE (μS/cm) NS Average daily streamflow (ft3/s) Base flow index, CaF Base flow index, OHS
01613000 11,050 10/1/1990 12/31/2020 CaF 0.971 0.58 4379.3 0.521
01613030 3,368 10/13/2011 12/31/2020 CaF 0.949 0.611 8.6 0.542
01613050 10,685 10/1/1991 12/31/2020 CaF 0.924 0.55 14.4 0.5
01613095 5,206 10/01/2006 12/31/2020 CaF 0.935 0.504 122.9 0.449
01613525 5,206 10/01/2006 12/31/2020 CaF 0.958 0.564 227.2 0.509
01613900 7,032 10/01/2001 12/31/2020 CaF 0.934 0.474 16.1 0.41
01614000 6,028 07/01/2004 12/31/2020 CaF 0.945 0.488 217.8 0.428
01614500 11,050 10/01/1990 12/31/2020 CaF 0.972 0.626 665.4 0.575
01615000 9,243 10/01/1990 12/31/2020 CaF 0.953 0.502 59.6 0.435
01616100 6,707 08/22/2002 12/31/2020 CaF 0.984 0.761 10.3 0.711
01616400 3,366 10/15/2011 12/31/2020 CaF 0.975 0.753 19.3 0.713
01616425 2,158 02/04/2015 12/31/2020 CaF 0.984 0.82 9.9 0.786
01616500 9,587 10/01/1994 12/31/2020 CaF 0.98 0.638 279 0.587
01617000 5,375 04/13/2006 12/29/2020 CaF 0.981 0.77 12.1 0.732
01617800 14,661 11/11/1980 12/31/2020 CaF 0.976 0.843 11.9 0.818
01618000 4,483 10/01/1990 12/31/2020 CaF 0.968 0.606 7551 0.549
01618100 4,542 03/27/2008 12/31/2020 CaF 0.988 0.897 13.8 0.874
01619000 5,665 06/29/2005 12/31/2020 CaF 0.982 0.76 108.5 0.722
01619500 14,447 06/13/1981 12/31/2020 CaF 0.987 0.784 316.9 0.749
01636316 6,706 08/23/2002 12/31/2020 CaF 0.986 0.824 21.6 0.782
01636464 4,657 04/02/2008 12/31/2020 CaF 0.99 0.856 20.4 0.829
01636500 9,589 10/01/1994 12/31/2020 CaF 0.973 0.666 2981.2 0.611
01636690 7,032 10/01/2001 12/31/2020 CaF 0.956 0.669 15.3 0.622
01637500 14,038 07/27/1982 12/31/2020 CaF 0.956 0.613 82.4 0.567
01638350 7,124 07/01/2001 12/31/2020 CaF 0.944 0.54 38.3 0.486
01638420 7,105 07/20/2001 12/31/2020 CaF 0.944 0.63 26 0.583
01638480 11,050 10/01/1990 12/31/2020 CaF 0.954 0.556 104.9 0.502
01638500 17,754 05/24/1972 12/31/2020 CaF 0.974 0.607 10332.6 0.55
01639000 11,050 10/01/1990 12/31/2020 CaF 0.932 0.398 239.6 0.351
01639500 14,581 01/30/1981 12/31/2020 CaF 0.969 0.612 122.7 0.563
01642190 5,206 10/01/2006 12/31/2020 CaF 0.965 0.521 970.2 0.471
01642198 1,219 1,192 08/31/2017 12/31/2020 Scfit and CaF 0.965 0.737 0.714 328.368 1740.572 0.717 6.6 0.668 0.685
01642199 1,178 1,136 10/11/2017 12/31/2020 Scfit and CaF 0.979 0.763 0.664 393.74 2030.648 0.691 15.3 0.605 0.672
01643000 17,738 06/09/1972 12/31/2020 CaF 0.969 0.519 1075.8 0.467
01643395 6,184 01/27/2004 12/31/2020 CaF 0.957 0.656 1.6 0.592
01643500 17,717 06/30/1972 12/31/2020 CaF 0.968 0.646 73.9 0.598
01643590 7,071 08/23/2001 12/31/2020 CaF 0.945 0.52 9.7 0.45
01643700 9,439 10/01/1990 12/31/2020 CaF 0.952 0.584 142.9 0.545
01643805 7,103 07/19/2001 12/31/2020 CaF 0.953 0.572 48.4 0.511
01643880 7,107 07/18/2001 12/31/2020 CaF 0.946 0.484 52.8 0.45
01644000 11,050 10/01/1990 12/31/2020 CaF 0.96 0.577 364.3 0.529
01644280 7,032 10/01/2001 12/31/2020 CaF 0.937 0.344 130.8 0.286
01644371 5,935 10/01/2004 12/31/2020 CaF 0.96 0.693 0.7 0.623
01644372 5,905 11/01/2004 12/31/2020 CaF 0.926 0.603 0.5 0.512
01644375 5,936 10/01/2004 12/31/2020 CaF 0.951 0.526 2.2 0.451
01644380 5936 10/01/2004 12/31/2020 CaF 0.955 0.676 1 0.602
01644388 2,760 06/12/2013 12/31/2020 CaF 0.945 0.501 4.3 0.442
01644390 3,716 10/28/2010 12/31/2020 CaF 0.951 0.497 6.2 0.431
01645000 17,747 05/31/1972 12/31/2020 CaF 0.967 0.616 141.8 0.556
01645704 4,830 4,663 10/01/2007 12/31/2020 Scfit and CaF 0.918 0.532 0.442 376.082 18156.766 0.587 8.9 0.369 0.423
01645762 4,839 4,560 10/01/2007 12/31/2020 Scfit and CaF 0.96 0.506 0.557 147.911 2225.673 0.569 3.2 0.477 0.445
01646000 4,840 3,268 10/02/2007 12/31/2020 Scfit and CaF 0.963 0.436 0.444 262.522 9251.594 0.518 71.6 0.378 0.373
Table 4.    Results of optimal hydrograph separation (OHS) and Collischonn and Fan (CaF) for accepted models.

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.

The observed and simulated daily specific conductance generally plots near the 1:1
                           line. The timeseries has an annual winter spike pattern
Figure 18.

Optimal hydrograph separation (OHS) results using the SCfit model for base-flow specific conductance (SC) for U.S. Geological Survey (USGS) gage 01642198, the upstream Carroll Creek gage, showing, A, simulated compared to observed SC, B, the streamflow hydrograph, and C, time series of simulated and observed SC values and simulated SC base flow. CaF, method developed by Collischonn and Fan (2013); CS, quickflow SC; Model p, model significance level of regression; R2, coefficient of determination.

The observed and simulated daily specific conductance generally plots near the 1:1
                           line. The timeseries has an annual winter spike pattern
Figure 19.

Optimal hydrograph separation (OHS) results using the SCfit model for base-flow specific conductance (SC) for U.S. Geological Survey (USGS) gage 01642199, the downstream Carroll Creek gage, showing, A, simulated compared to observed SC, B, the streamflow hydrograph, and C, time series of simulated and observed SC values and simulated SC base flow. CaF, method developed by Collischonn and Fan (2013); CS, quickflow SC; Model p, model significance level of regression; R2, coefficient of determination.

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. 20AC) and (2) the base flow index (fig. 20DF). 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 PART method has slightly higher estimates for both baseflow per unit area and
                           baseflow index when compared to OHS and CaF
Figure 20.

Comparison of computed base flow per unit area per year (AC) and base flow index (DF) for 52 streamgages in the region surrounding Fort Detrick, Area B for water years 2017–20. OHS, optimal hydrograph separation; CaF, method developed by Collischonn and Fan (2013).

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.

Water Level Results

Groundwater levels observed within the groundwater monitoring network reveal a highly responsive groundwater system at Fort Detrick Area B. Daily mean groundwater levels across the site (fig. 16E) show rapid rises to individual storm events that can be on the order of several feet. The water levels time series show a similar pattern to the stream base flow (fig. 16C). The least responsive water level time series were collected at FR Dd 221 (BMW65A) and FR Dd 222 (BMW65B), located on the east edge of Area B. Large drawdowns in water level that take several days to recover are observed when these locations are sampled, leading to gaps in the observational record. These drawdowns, along with limited responsiveness to individual storm events may indicate these wells are not as connected to larger fractures and voids. However, the seasonal to interannual variability within these wells does agree with the rest of the monitoring network.

Specific Yield

Thirteen periods of recession were identified from which the apparent specific yield could be calculated (table 5). The recession period streamflow from the lower Carroll Creek gage (01642199) is described in table 5.

Table 5.    

Recession periods from the lower Carroll Creek streamgage (01642199) used to calculate specific yield.

[MM, two-digit month; DD, two-digit day, YYYY, four-digit year]

Recession period index Recession period start (MM/DD/YYYY) Recession period end (MM/DD/YYYY) Number of days Change in streamflow (cubic feet per second) Cumulative streamflow during recession period (inches per unit area)
1 10/02/2018 10/10/2018 9 −11 0.95
2 12/04/2018 12/12/2018 9 −3.9 0.76
3 03/27/2018 04/04/2018 9 −0.05 0.4
4 03/12/2019 03/20/2019 9 −4.6 0.76
5 10/10/2019 10/15/2019 6 −0.49 0.12
6 11/13/2019 11/23/2019 11 −0.55 0.29
7 03/26/2019 04/04/2019 10 −8.8 0.94
8 02/14/2020 02/24/2020 11 −3.4 0.64
9 10/03/2020 10/10/2020 8 −0.15 0.17
10 10/16/2020 10/25/2020 10 −0.44 0.27
11 11/02/2020 11/08/2020 7 −0.58 0.2
12 11/17/2020 11/22/2020 6 −0.56 0.2
13 12/06/2020 12/13/2020 8 −2.9 0.33
Table 5.    Recession periods from the lower Carroll Creek streamgage (01642199) used to calculate specific yield.

Using the criteria defined for this project, 6 of the 12 available groundwater-level time series had suitable data completeness and regression model fit statistics so that the apparent specific yield could be calculated. The fit statistics for the linear regression models are shown in table 6; the average of the calculated specific yields is 0.067 and the median value is 0.06. Only one well, FR Dd 219 (BMW53B), had an apparent specific yield confidence interval that included the mean and median values, perhaps indicating a bimodal distribution of apparent specific yields, although the number of available values was too few to make this determination. Olmsted and Hely (1962) calculated a value of 0.063–0.074 using wells within the Brandywine Creek basin in Pennsylvania, a 287 mi2 watershed underlain by fractured rock and carbonate units. Risser and others (2005) reported an average specific yield of 0.013 in fractured and folded shales, sandstones, and siltstones, though the apparent specific yields in upland wells ranged from 0.0035 to 0.035. Delle Rose and others (2018) apply this method to a deep, highly developed karst aquifer in the Salento peninsula, Italy, a region with extensive cave development, and report a much higher apparent specific yield of 0.45–0.49. In Berkely County, West Virginia, Shultz and others (1995) calculated specific yields using stream base flow records and groundwater level fluctuations in a 60 mi2 area underlain by carbonate rocks; the calculated specific yields ranged from 0.044 to 0.049. Based on the results from these studies, the observed median value of 0.06 and a range of 0.04 to 0.11 is consistent with the storage and release of water from secondary porosity (fractured rock) and relatively poorly developed karstic tertiary porosity (with drainage from solutionally enlarged fractures but not large caves). Due to the variability in observed specific yields, an uncertainty of ±50 percent to the median apparent specific yield value is applied in the water budget analysis, a range that includes the lowermost 95-percent confidence interval value (0.03) and the uppermost 95-percent confidence interval value (0.12).

Table 6.    

Calculated specific yield values for six monitoring wells, Fort Detrick, Area B.

[R2, coefficient of determination; p, significance level of regression; <, less than]

Name Number of recessions Apparent specific yield 95-percent confidence interval lower bound 95-percent confidence interval upper bound Adjusted R2 p-value
FR Dd 220 (BMW53F) 11 0.04 0.04 0.05 0.95 <0.05
FR Dd 219 (BMW53B) 12 0.07 0.05 0.11 0.73 <0.05
FR Dd 223 (BMW66B) 12 0.11 0.1 0.12 0.99 <0.05
FR Dd 224 (BMW66D) 12 0.09 0.08 0.11 0.96 <0.05
FR Dd 225 (BMW77) 12 0.04 0.03 0.04 0.91 <0.05
FR Dd 226 (BMW78B) 12 0.05 0.04 0.05 0.97 <0.05
Table 6.    Calculated specific yield values for six monitoring wells, Fort Detrick, Area B.
Groundwater Gradients

Hydraulic gradients observed between monitoring locations at Fort Detrick Area B were consistently from west to east and generally upwards, although the vertical well nest containing FR Dd 220 (BMW53F) and FR Dd 219 (BMW53B) showed a consistent downward gradient. Daily mean horizontal hydraulic gradients, daily mean vertical hydraulic gradients, and a schematic showing the well locations are shown in figure 21A–C, respectively. The greatest west-to-east hydraulic gradients were observed during the period of high-water levels in late 2018 and early 2019. These gradients decreased in late 2019 during a drier period, then steepened again slightly in the spring of 2020 during a period of higher water levels and stream base flow (fig. 16). During the fall of 2020, gradients again declined. These results indicate an influence of hydrologic regime on the observed west-east hydraulic gradient.

From mid-2018 to the end of 2020, west to east gradients decline in magnitude in eight
                              select well pairs. One of four well nests decreases
Figure 21.

A, horizontal hydraulic gradients and, B, vertical hydraulic gradients between, C, selected well pairs. In A, negative is from west to east, and in B, positive is upwards. Ft, feet.

Vertical hydraulic gradients at five nested well pairs are shown in figure 21B. The nested well pair of FR Dd 222 and FR Dd 221 (BMW65B and BMW65A) are influenced by pumping and show the most delayed and muted of the water-level responses within the hydrologic network (fig. 16E). The highly variable vertical hydraulic gradient in response to storm events is explained by a lag in the response time for deeper well FR Dd 222 (BMW65B) relative to shallower well FR Dd 221 (BMW65A). During most of the monitoring period, the shallow well responded first, briefly resulting in a decreased vertical upward gradient until the deeper well responded. Small downward gradients were observed at this well network in late 2018, though for most of the period upward gradients were observed. The upward gradient was greatest during relatively dry periods (fall 2019 and 2020) and least during relatively wet periods (early 2019 and 2020), which stands in contrast to the pattern observed in three well pairs: FR Dd 224/223 (BMW66B/D), FR Dd 219/220 (BMW53B/F), and FR Dd 248/249 (BMW32/BMW32D). In those well nests, during relatively wet periods the upward gradient steepens and during relatively dry periods the upward gradient decreases or reverses entirely. Wells FR Dd 219/220 (BMW53B/F) had a slight upward vertical gradient during the extremely wet conditions at the beginning of the monitoring period. Since that time, a slight downward gradient was observed. This well nest is located within a low-lying ephemerally wet marshy area in the center of Area B. In late 2019, artesian conditions were observed in well FR Dd 220 (BMW53F). Consistent upward gradients were observed in well pair FR Dd 246/247 (BMW47/BMW47D). The magnitude of these vertical gradients generally agrees with the well nest FR Dd 221/FR Dd 222 (BMW65A/B) that is also located on the west edge of Area B.

EMR Results

EMR analysis was successfully applied to 10 of 12 groundwater-level time series, allowing estimates of groundwater recharge to be calculated from the well data. In addition to quantifying the episodic recharge entering the aquifer during storm events and overall aquifer responsiveness, this analysis also characterizes the physical properties of the aquifer through recession slope analysis. Wells FR Dd 221 (BMW65A) and FR Dd 222 (BMW65B) were excluded from this analysis, as the muted water-level responses and influence of groundwater pumping make these water-level records unsuitable for evaluation. The master recession curves for all wells are shown in figure 22.

Seven of the master recession curves are linear relations. Three indicate that at
                              high elevations a much steeper recession slope is observed
Figure 22.

Master recession curves developed for 10 water-level time series collected at wells at Fort Detrick Area B.

Seven of the ten wells have simple linear relation between recession slope and water level elevation, which can be fit with linear regression (fig. 22). However, three wells show a break in slope in the master recession curve, with higher water table elevations having a much faster rate of recession than when the water table is lower: FR Dd 223 (BMW66B), FR Dd 224 (BMW66D), and FR Dd 226 (BMW78B). A break in slope in the master recession curve can be related to different physical processes, including the drainage of lithologies with different specific yields and the draining of fractures of voids within shallow karst aquifers (Nimmo and Perkins, 2018). A water table starting within the epikarst or a zone of high fracture transmissivity and lowering to a zone with less transmissivity can produce these effects. Due to the highly localized nature of preferential flow features, which may or may not intersect the monitoring well screen, predicting the recession behavior of unmonitored wells or unmonitored portions of the aquifer is problematic. However, this analysis shows the presence of this process.

Water-level fluctuations in response to storm events provide information about the water recharged to the groundwater system episodically during individual storm events. Nimmo and Perkins (2018) introduced a set of codes in the R language environment that can rapidly and reproducibly convert records of water level and precipitation observations to estimates of groundwater recharge using a set of adjustable parameters. For this study, the R scripts were applied to 10 groundwater-level time series recorded at Fort Detrick Area B. Wells FR Dd 221 (BMW65A) and FR Dd 222 (BMW65B) were excluded from this analysis, as the muted water level responses and influence of groundwater pumping make these water level records unsuitable for this analysis. The precipitation data used in this analysis was recorded by the USGS at the upstream Carroll Creek streamgage. Prior to the analysis, which requires continuous groundwater level records, the instantaneous records of groundwater level and precipitation were resampled to 15-minute intervals using linear interpolation and converted to fractional days since October 1, 2018, the beginning of many of the groundwater observations. Following the procedure outlined in Nimmo and Perkins (2018), a manual adjustment process for each site took place to ensure that (1) recession periods were adequately characterized to create a master recession curve relating the water table elevation and recession slope, and (2) recharge events were adequately characterized to capture the water table rise amounts, recorded as a linear distance. The parameters used to perform this analysis are shown in table 7. The master recession curves for all wells in this analysis are shown in figure 22. The specific yield is needed to convert the water table rises to estimates of recharge. The specific yields reported in table 6 were used for the six wells in the table; the remaining wells were assigned the median specific yield value of 0.06.

Table 7.    

Parameters used to apply the methods of Nimmo and Perkins (2018) to 10 water-level time series collected at Fort Detrick Area B.

[MRC, master recession curve; EMR, episodic master recession, NAVD88, North American Vertical Datum of 1988; dH/dt, change in elevation (dH) per unit of time (dt); —, not applicable]

Parameter Units FR Dd 223
FR Dd 224
FR Dd 225
FR Dd 220
FR Dd 219
FR Dd 226
FR Dd 246
FR Dd 247
FR Dd 248
FR Dd 249
MRC parameters
Piecewise master recession curve? Yes Yes No No No Yes No No No No
Breakpoint of master recession curve Feet above NAVD88 322.5 326.5 347.5
Force fit through origin? No No No No No No No No No No
Minimum duration of interval between significant precipitation and recession Days 0 0 0 0 0 0 0 0 0 0
Negligibility-of-precipitation criterion Inches 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
Duration of slope elements for linearization Days 1 1 1 1 1 1 1 1 1 1
Maximum allowable total uptick within a linearized slope element Feet 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
Desired bin size (units of response) to bin the calculated dH/dt for fitting Feet 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25
Maximum allowable response-increase for linearized slope elements to be used in fitting MRC Feet/day 0 –0.02 0 0 0 0 0 0 –0.02 –0.02
EMR parameters
Lag time between input and response Days 0 0 0 0 0 0 0 0 0 0
Fluctuation tolerance (the maximum dH/dt value allowable as system noise) Feet/day 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
Minimum amount of precipitation within precipitation window to allow inclusion as an episode Inches 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25
Parameter used in determining the end of an episode, time allowed for stabilization Days 0.1 0.1 0.3 0.3 0.3 0.1 0.3 0.3 0.3 0.3
Number of points used to smooth hydrograph-slope values 4 4 4 4 4 4 4 4 4 4
Table 7.    Parameters used to apply the methods of Nimmo and Perkins (2018) to 10 water-level time series collected at Fort Detrick Area B.
Changes in Groundwater Storage

The relation between the monthly groundwater storage deviation and the monthly mean base flow computed using OHS is well-fit with a log-linear regression (fig. 23), indicating that the representation of groundwater storage in the groundwater monitoring network and stream base flow are highly related. The suitability of the log-linear regression ensures that changes in storage from the well network reasonably approximate the storage changes observed throughout the portion of the watershed draining to Carroll Creek. The regression equations shown in figure 23 allow the monthly groundwater storage to be estimated during the entire water budget period (fig. 24). The estimated change in groundwater storage for the water budget periods using the regression equations is shown in table 8, along with the estimated uncertainty. The estimated uncertainty is derived from the 50-percent uncertainty on the specific yield.

Two linear regression relations indicate that increasing storage is associated with
                              increasing baseflow discharge
Figure 23.

Regression relations between monthly groundwater level storage deviation and monthly stream base flow. RMSE, root-mean-square error; OHS, optimal hydrograph separation; R2, coefficient of determination.

Deviation increases from April to May 2018, then declines from May to October 2019.
                              A smaller increase is observed in January–April 2020
Figure 24.

Monthly estimated groundwater storage at Fort Detrick Area B, shown as deviation from the mean value.

Table 8.    

Change in storage for key water budget periods.

Name Change in storage from downstream regression, in inches per year Upper uncertainty bound, in inches per year Lower uncertainty bound, in inches per year
Water year 2018 5.79 8.69 2.9
Water year 2019 −5.86 −8.79 −2.93
Water year 2020 0.13 0.2 0.07
Water year 2018–20 0.02 0.03 0.01
June 1, 2017–December 31, 2020 −1.32 −1.98 −0.66
October 1, 2017–December 31, 2020 0.94 1.41 0.47
Table 8.    Change in storage for key water budget periods.

The values reported in figure 24 and table 8 show that historically wet weather in early to middle 2018 led to a large increase in groundwater storage from the beginning of the water budget period, which by the end of 2020 had mostly returned to the original level of groundwater storage. High levels of groundwater storage are observed throughout the second half of 2018 and the first half of 2019. By the second half of 2019 through the end of the monitoring period, an annual pattern typical of the mid-Atlantic is shown; low amounts of groundwater storage in the fall and winter and increased groundwater storage in the spring. The net change in groundwater storage from the beginning (October 2017) to the end of the observation period (December 2020) was 2.97 in., a net increase of approximately 1 in/yr.

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.

The modeled and observed values scatter about the 1:1 line and within an error bound
                           of 15 percent, except for values below 5
Figure 25.

Monthly total observed precipitation at the U.S. Geological Survey (USGS) precipitation gage and modeled by Parameter-elevation Regressions on Independent Slopes (PRISM) from May 1, 2018, through December 31, 2020, the period for which both measures are available.

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.

Discrete Monitoring

The 1-liter samples collected as split samples from quarterly monitoring and spring grab samples were analyzed for rhodamine and fluorescein on the Turner Scientific C6 Fluorometer in a laboratory setting. The apparent dye concentrations associated with varying levels of turbidity were also evaluated in a benchtop experiment, which indicated that low-level interference from turbidity was observed with this instrumentation (fig. 26A and B). The apparent rhodamine observed from turbidity interference is greater than the apparent fluorescein, though the relation between rhodamine and turbidity is better fit with a linear function than fluorescein. Although there were slight differences in the slope of the regression lines, overall, the magnitude of the interference for both the March 2019 and October 2020 tests were similar, showing good reproducibility of this evaluation.

Apparent concentrations in the discrete data are similar to expectations, except for
                              well BMW67C, which is much higher
Figure 26.

Results of benchtop turbidity experiment and discrete fluorometric sample analysis. The line shown in C and D is from the combined regression equations in A and B. The sample from BMW67C is not plotted in C and D; in that sample, 487.47 micrograms per liter (μg/L) fluorescein and 406.28 μg/L rhodamine were observed, and the sample was visibly pink. R2, coefficient of determination.

The observed dye concentrations in discrete samples collected for this study as a function of the sample turbidity are shown in figure 26C and D; the results are available in tabular format in table 9. The combined regression lines shown in figure 26A and B are also shown in figure 26C and D. All samples except BMW67C have low-level observed apparent dye concentrations that are explained by the natural particulates within the water; dye is not detected at a level above the instrument’s ability to resolve the dye. Concentrations of 487.47 μg/L fluorescein and 406.28 μg/L rhodamine were recorded in BMW67C. This well was the site of Eosine dye injections in the 2013 dye tracer study. It is likely that the recorded concentrations of fluorescein and rhodamine reflect the interference of Eosine dye in this well. To evaluate this, a 5,000 µg/L standard solution of eosine was evaluated on the fluorometer used in this study and the concentrations of apparent fluorescein and rhodamine were recorded. Apparent concentrations of 480.0 µg/L fluorescein and 400.0 µg/L rhodamine were observed, which indicates that Eosine interferes with the detection of these dyes in this well.

Table 9.    

Discrete sample results from split well samples analyzed for fluorescein and rhodamine dyes.

[MM, two-digit month; DD, two-digit day; YYYY, four-digit year; µg/L, micrograms per liter; FNU, Formazin Nephelometric Unit; °C, degrees Celsius; —, no note recorded]

Sample name Sample date and time (MM/DD/YYYY and 24-hour time) Sample type Notes Measured fluorescein (µg/L) Measured rhodamine (µg/L) Turbidity (FNU) Temperature (°C) Fluorescein explainable by turbidity using regression
y = 0.0054*turbidity + 0.0985
Rhodamine explainable by turbidity using regression
y = 0.0254*turbidity + 0.0665
BMW67C 12/01/2020 1525 Well split Visibly green and pink 478.47 406.28 0.01 18.49 0.1 0.07
BMW32D 12/04/2020 1035 Well split Fine orange sediment 0.49 1.44 66.71 18.49 0.46 1.76
BMW69A 12/03/2020 1305 Well split Fine orange sediment 0.45 0.84 35.14 18.54 0.29 0.96
Hospital Spring #1 12/09/2020 0902 Spring 0.25 0.85 36.6 18.78 0.3 1
BMW56D 12/01/2020 1250 Well split 0.16 0.11 0.05 18.61 0.1 0.07
BMW58D 12/03/2020 1310 Well split 0.16 0.14 0.2 18.51 0.1 0.07
South Armory Spring #2 12/09/2020 0850 Spring Duplicate sample 0.14 0.38 17.52 18.72 0.19 0.51
Site 5 Spring #1 12/09/2020 0912 Spring 0.12 0.24 8.13 18.76 0.14 0.27
BMW59D 12/02/2020 0940 Well split 0.1 0.11 0.05 18.37 0.1 0.07
BMW77 12/02/2020 1105 Well split 0.08 0.08 0.09 19.55 0.1 0.07
Hospital Spring #2 12/09/2020 0902 Spring Duplicate sample 0.08 0.12 1.45 18.77 0.11 0.1
BMW70A 12/02/2020 1450 Well split 0.06 0.05 0.08 18.41 0.1 0.07
Site 5 Spring #2 12/09/2020 0912 Spring Duplicate sample 0.06 0.17 6.31 18.97 0.13 0.23
BMW20D 12/02/2020 1405 Well split 0.05 0.06 0.52 18.35 0.1 0.08
South Armory Spring #1 12/09/2020 0850 Spring 0.05 0.16 6.83 18.43 0.14 0.24
BMW53F 12/03/2020 1055 Well split 0.03 0.03 0.05 18.57 0.1 0.07
BMW71C 12/03/2020 1545 Well split 0.03 0.03 0 18.5 0.1 0.07
BMW24D 12/03/2020 1005 Well split 0.03 0.05 0 18.52 0.1 0.07
Table 9.    Discrete sample results from split well samples analyzed for fluorescein and rhodamine dyes.
Continuous Monitoring

Continuous monitoring at the three springs was conducted sequentially between November 2019 and October 2020. The complete time series of data collected at each site is shown in figures 2729. All three springs had stable temperatures between 12.5 and 14.5 °C and pH values between 7.1 and 7.4, which is consistent with groundwater discharging from carbonate aquifers (Hem, 1985). Two of the springs (RISP-3 and Hospital Spring) are located within completed spring boxes, whereas South Armory Spring is in a channel originating approximately 20 ft from Carroll Creek. To maintain a sufficient water level for in situ monitoring, a temporary weir from sandbags was constructed. Intermittent leaks or rearrangement of this weir from Carroll Creek storm events led to some gaps in the monitoring period. The difference between the springs in completed box structures and the South Armory Spring is evident in the dissolved oxygen time series. RISP-3 and the Hospital Spring had diurnal fluctuations in dissolved oxygen consistent with plant and algal activity, whereas no such pattern exists in the South Armory Spring. The difference is also observable in the turbidity; the completed spring boxes had turbidities below 5 formazin nephelometric units (FNU) for the duration of the monitoring period and South Armory Spring had reported turbidities exceeding 60 FNU.

Overall, there is general consistency in all observed elements during the monitoring
                              period. The dye concentration remained at a low level
Figure 27.

A, observed temperature, B, turbidity, C, specific conductance, D, dissolved oxygen, E, pH, F, apparent fluorescein, and, G, apparent rhodamine concentrations at spring FR Dd 213 (RISP-3). Two sensors were deployed at this location; the results of both sensors are provided.

Overall, there is general consistency in observed elements outside some short-lived
                              storm responses. The dye concentration remained low
Figure 28.

A, observed temperature, B, turbidity, C, specific conductance, D, dissolved oxygen, E, pH, F, apparent fluorescein, and, G, apparent rhodamine concentrations at spring FR Dd 178 (Hospital Spring). Two sensors were deployed at this location; the results of both sensors are provided.

Temperature increased, turbidity had one multi-day increase, four short-lived decreases
                              in specific conductance, but consistency in all else
Figure 29.

A, observed temperature, B, turbidity, C, specific conductance, D, dissolved oxygen, E, pH, F, apparent fluorescein, and, G, apparent rhodamine concentrations at spring FR Dd 245 (South Armory Spring). Two sensors were deployed at this location; the results of both sensors are provided.

Hospital Spring had a much higher specific conductance (775–875 μS/cm) than the other two springs (400–500 μS/cm) for the duration of the monitoring period. Freshening events from rainfall were apparent in all three springs, though these were much more pronounced in the South Armory and Hospital Springs. These freshening events may represent a mixture between groundwater and recharge from precipitation that represents quickflow or lateral vadose zone flow within the soil and (or) epikarst.

The observed apparent fluorescein and rhodamine concentrations are very low, with all concentrations below 1 μg/L for both dyes. These apparent detections must be considered in the context of the influence of suspended particulates (turbidity), which is explored further in the preceding section. The low-level and spiky nature of the apparent dye detections suggest that dyes are not observed in situ at a level greater than the ability of the sensor to detect the dyes in natural waters found at Fort Detrick.

Soil-Water Balance Model Results

The SWB model created for this study used the best available inputs of land cover, land use, climate, topography, and process-based information controlling the routing of water within the Fort Detrick watersheds. Additional detail on the implementation of the model including the complete outputs are available in an accompanying data release (Soroka, 2023). For the water budget calculation in this study (eq. 7), two water budget components (ET and interception) are assessed using the SWB. An additional SWB output, net infiltration, is presented though is not used in the water budget. This section is focused on the SWB model results at the watershed scale and timescales of 1 year or greater, as this is the scale relevant to the water budget calculation

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 mi2 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.

Precipitation and infiltration values agree with observed streamflow while simulated
                           net infiltration and baseflow do not
Figure 30.

Comparisons between soil-water balance model outputs, streamflow, and stream base flow computed at, A, 36 watersheds within the region surrounding Fort Detrick for, B, observed streamflow, C, stream base flow using the CaF method and, D, stream base flow using the PART method. USGS, U.S. Geological Survey; RMSE, root-mean-square error; NSE, Nash-Sutcliffe efficiency coefficient; CaF, method developed by Collischonn and Fan (2013).

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.

The budget values for precipitation show a northwest-southeast increasing gradient
                           while the actual values are more variable
Figure 31.

A, annualized gridded values for gross precipitation, B, actual evapotranspiration, C, net infiltration, and, D, interception from October 1, 2017, to December 31, 2020, for components used in the water budget calculation.

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 and net infiltration model outputs are sensitive to changes
                           in the soil curve numbers used for soil groups 1 and 4
Figure 32.

Results of a sensitivity analysis examining the sensitivity of changes of select model parameters to select model 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 = PQ / 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 PQ / 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 PQ / A – ΔS estimates.

The estimates derived from the SWB model are of a comparable magnitude to an estimate
                           derived from other measured water budget components
Figure 33.

Comparison between estimated actual evapotranspiration between the soil-water balance (SWB) model used in this study and four alternative estimates of evapotranspiration. Aet, actual evapotranspiration; intercep., interception; P, precipitation; Q, stream discharge; A, watershed area; ΔS, change in groundwater storage.


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.

Estimated recharge varies, with SWB consistently the lowest, then stream base flow
                           per unit area, RORA, and the episodic master recession
Figure 34.

Comparison between estimated net infiltration from the soil-water balance (SWB) model used in this study to other data sources. EMR, episodic master recession; OHS, optimal hydrograph separation; RORA, RORA computer program.

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.

Water Budget Results

The water budget described in equation 7, as well as the estimated uncertainties used to calculate the uncertainty in the water budget residual (eq. 8) were calculated in this study for six time intervals. Four of these are defined by the water year, which begins on October 1 and is often used in water budget studies as autumn is typically a period of low streamflow. These four intervals are (1) the 3 water-year period (October 1, 2017–September 30, 2020), (2) water year 2018, (3) water year 2019, (4) water year 2020. The budget was also computed for the full monitoring period (October 1, 2017–December 31, 2020), and the period following the anomalous May 2018 flooding (June 1, 2018–December 31, 2020).

Uncertainty in the water budget components was estimated using available information in the hydrologic monitoring, SWB model sensitivity evaluation, and when possible, through comparisons with independent datasets. Percentage uncertainties are assigned on the basis of this available information in a manner comparable to other studies (Maupin and Weakland, 2009; Winter, 1981). An uncertainty of 10 percent was assigned to streamflow totals, 5-percent uncertainty to precipitation totals, 30-percent uncertainty to actual ET calculated by the SWB model, and 10-percent uncertainty to interception calculated by the SWB model. Uncertainty in groundwater storage changes was estimated to be 50 percent.

The calculation results are shown in table 10 for the six time periods. For each time period, the calculation was completed for three spatial areas: the upstream Carroll Creek watershed, the downstream Carroll Creek watershed, and the watershed difference area, which represents the area within the downstream watershed that is not in the upstream watershed. The table presents the annualized budget components in inches per year. The water budget residual, as defined in equation 7 and the uncertainty on the residual as defined in equation 8 are presented for each time period and each spatial area. The interpretation of the water budget residual is only meaningful when the absolute value of the residual exceeds the absolute value of the uncertainty. Meaningful residuals are indicated with a footnote in table 10.

Table 10.    

Water budget summary table.

[in/yr, inches per year; P, precipitation; ET, evapotranspiration; I, interception; ΔS, change in groundwater storage; Q, stream discharge; R, residual; SWB, soil-water balance; PRISM, Parameter-elevation Regressions on Independent Slopes]

Components Upstream watershed Downstream watershed Watershed difference area
Equation term Data source Uncertainty applied (percent) Value
Water years 2018–2020
(October 1, 2017–September 30, 2020)
P PRISM 5 51.4 2.6 51.3 2.6 51.1 2.6
ET SWB 30 23.4 7 23.5 7.1 23.7 7.1
I SWB 10 2.6 0.3 2.7 0.3 2.8 0.3
Q Observed 10 22.6 2.3 29.6 3 39.3 3.9
ΔS Calculated 50 0 0 0 0 0 0
(P − ET − I − Q − ΔS = R)
(eq. 7)
Sum of squared errors
(eq. 8)
2.8 7.8 −4.5 8.1 −14.7a 8.5
Water year 2018
(October 1, 2017–September 30, 2018)
P PRISM 5 67.8 3.4 67.1 3.4 66.1 3.3
ET SWB 30 30.2 9.1 30.5 9.2 30.9 9.3
I SWB 10 3.7 0.4 3.8 0.4 3.9 0.4
Q Observed 10 24.2 2.4 33.7 3.4 46.9 4.7
ΔS Calculated 50 5.8 2.9 5.8 2.9 5.8 2.9
(P − ET − I − Q − ΔS = R)
(eq. 7)
Sum of squared errors
(eq. 8)
3.9 10.4 −6.7 10.7 −21.5a 11.3
Water year 2019
(October 1, 2018–September 30, 2019)
P PRISM 5 55.3 2.8 55.9 2.8 56.7 2.8
ET SWB 30 28.1 8.4 28.2 8.5 28.4 8.5
I SWB 10 2.9 0.3 3 0.3 3.2 0.3
Q Observed 10 32.1 3.2 37.7 3.8 45.6 4.6
ΔS Calculated 50 −5.9 −2.9 −5.9 −2.9 −5.9 −2.9
(P − ET − I − Q − ΔS = R)
(eq. 7)
Sum of squared errors
(eq. 8)
−2 9.9 −7.3 10.1 −14.7a 10.5
Water year 2020
(October 1, 2019–September 30, 2020)
P PRISM 5 31.1 1.6 30.9 1.5 30.6 1.5
ET SWB 30 11.8 3.5 11.8 3.5 11.8 3.6
I SWB 10 1.2 0.1 1.3 0.1 1.3 0.1
Q Observed 10 11.6 1.2 17.3 1.7 25.3 2.5
ΔS Calculated 50 0.1 0.1 0.1 0.1 0.1 0.1
(P − ET − I − Q − ΔS = R)
(eq. 7)
Sum of squared errors
(eq. 8)
6.4a 4 0.4 4.2 −8a 4.6
Total observation period
(October 1, 2017–December 31, 2020)
P PRISM 5 47.4 2.4 47.3 2.4 47.1 2.4
ET SWB 30 21.6 6.5 21.9 6.6 21.9 6.6
I SWB 10 2.4 0.2 2.5 0.2 2.6 0.3
Q Observed 10 21.6 2.2 28.6 2.9 38.3 3.8
ΔS Calculated 50 0.9 0.5 0.9 0.5 0.9 0.5
(P − ET − I − Q − ΔS = R)
(eq. 7)
Sum of squared errors
(eq. 8)
0.9 7.2 −6.6 7.6 −16.6a 8
Period following May 2018 floods
(June 1, 2018–December 31, 2020)
P PRISM 5 47.1 2.4 47.3 2.4 47.4 2.4
ET SWB 30 22.5 6.7 22.7 6.8 22.9 6.9
I SWB 10 2.6 0.3 2.6 0.3 2.8 0.3
Q Observed 10 23.3 2.3 29.6 3 38.3 3.8
ΔS Calculated 50 −1.3 −0.7 −1.3 −0.7 −1.3 −0.7
(P − ET − I − Q − ΔS = R)
(eq. 7)
Sum of squared errors
(eq. 8)
0.1 7.6 −6.3 7.8 −15.2a 8.2
Table 10.    Water budget summary table.

Meaningful residual, such that the absolute value exceeds the absolute value of the associated uncertainty.

Water budget residuals differ between the upstream, downstream, and watershed difference areas during the six water budget time intervals. These differences during a given time period are largely a result of differences in the streamflow per unit area. For example, during the water year 2018–water year 2020 period, the range in the streamflow per unit area varied from 22.6 in/yr for the upstream watershed to 39.3 inches per area for the watershed difference area. The precipitation, ET, and interception used in this budget vary little between these three areas and the change in storage is assumed to vary uniformly across the watersheds.

The watershed difference area consistently has a large magnitude negative water budget residual across all six water budget time periods, indicating a greater loss of water per unit area than gain per unit area. The significance of this residual indicates that a water budget term is missing or misrepresented. There are two primary possible explanations for the large magnitude residual. The first is that focused groundwater recharge within karst features located in western Area B allows water to quickly enter the groundwater system, leading precipitation to bypass vegetation and the soil zone, and ultimately discharging to Carroll Creek through springs. Rapid recharge in this process would result in the ET component used in this study to be an overestimate. An ET component with approximately half the magnitude would be required for the watershed difference area residual to be less than the uncertainty. The second possible explanation is that the interbasin transfer of water brings water from outside the watershed difference area defined by the surface topography. As no anthropogenic transfers are present, this would likely occur through interbasin groundwater flow where the groundwater is recharged from an area outside the topographic watershed boundary. Depending on the time interval, the interbasin transfer component would need to be 3.4 to 10.1 in/yr per unit area for the watershed difference area residual to be less than the uncertainty.

The water budget residuals for the upstream Carroll Creek watershed were not significant for the 3 water-year period, water year 2018, and water year 2019. However, during water year 2020, the driest interval, the upstream watershed has a significant positive water budget residual. Similar to the above discussion, the cause of these significant residuals could be a result of a misrepresentation of one of the water budget components or interbasin transfer. In the case of the upper watershed, a positive residual is observed, suggesting that either ET is underestimated or water leaves the basin through groundwater flow. During the water budget periods examined, the downstream watershed did not have any significant water budget residuals.

Geochemical Interpretation Results

The geochemical interpretation in this study seeks to evaluate the recharge source(s) and degree of connectivity between groundwater flow-paths underlying Fort Detrick Area B. The geochemical data used in this study includes major ion composition and groundwater tracers including dissolved noble gases, isotopes of hydrogen and oxygen in water, and age tracers. Major ion composition was used to characterize the primary water-rock interactions controlling groundwater chemistry and identify patterns of groundwater chemistry across the site. Tracers were used to identify recharge sources, estimate residence times, and further refine the understanding of groundwater flow-paths. Results for each of the approaches provide specific information about the system and were interpreted holistically to inform the conceptual site model.

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.

Table 11.    

Select groundwater field parameters, major ion concentrations, trace metal concentrations, and redox classification, Fort Detrick Area B, Maryland.[—Left]

[USGS, U.S. Geological Survey; ID, identification; °C, degrees Celsius; µS/cm, microseimens per centimeter at 25 °C; DO, dissolved oxygen; mg/L, milligrams per liter; N, nitrogen; P, phosphorus; SiO2, silica dioxide; µg/L, micrograms per liter; <, less than; Cl:Br, mass ratio of chloride to bromide; —, no note recorded]

USGS site number Sample name Site ID Temperature (°C) Specific conductance (µS/cm) pH Field DO (mg/L) Bicarbonate (mg/L) Nitrate as N (mg/L) Orthophosphate as P (mg/L) Calcium (mg/L) Magnesium (mg/L) Sodium (mg/L) Potassium (mg/L) Chloride (mg/L) Bromide (mg/L) Sulfate (mg/L) Fluoride (mg/L)
392556077263301 RISP-3 FR Dd 213 13.3 479 7.7 5.7 223 4.97 0.039 64.3 17 11.4 1.66 29.9 0.023 9.2 0.04
392600077271001 BMW53B FR Dd 219 14.2 493 7.9 4.7 246 3.93 0.012 62.4 18.2 14 1.58 34.4 0.024 6.94 0.03
392608077273602 BMW65B FR Dd 222 17.2 1,250 8.3 0 266 0.04 0.026 30.1 4.5 200 3.15 24.5 0.111 293 0.11
392611077270801 BMW78B FR Dd 226 15.1 453 8 4.9 235 3.28 0.017 58.5 19.9 8.36 2.14 18.4 0.021 11.7 0.04
392559077271002 BMW53F FR Dd 220 13.6 470 8 5.1 213 4.23 0.009 54.1 20.8 10.3 1.09 38.7 0.029 7.13 0.02
392552077261601 County Shallow FR Dd 234 14.8 1,030 7.6 5 287 4.36 0.014 116 14 76.7 4.16 164 0.072 26.1 0.07
392553077264503 BMW70C FR Dd 232 13.8 734 8 2.3 296 3.09 0.004 110 13.7 27.5 3.18 65.2 0.031 25.2 0.03
392610077270801 BMW78A FR Dd 233 14.4 449 8.3 5.6 239 3.87 0.022 67.1 16.1 8.22 2.68 18.3 0.018 10.5 0.04
392553077264501 BMW70A FR Dd 231 13.8 627 7.7 0.3 288 2.31 0.009 97.6 9.43 21.9 2.29 42 0.029 21.1 0.03
392603077263701 BMW35 FR Dd 229 13.9 482 7.7 5.5 219 5.12 0.041 65.6 17 11.3 1.81 30.2 0.022 9.4 0.04
392603077263702 BMW35D FR Dd 230 13.7 476 7.7 5.5 213 4.98 0.047 63 16.1 11.1 1.78 30.8 0.021 9.28 0.03
392548077255101 AMW 568-13A FR Dd 238 16.6 674 7.9 0.8 256 2.29 0.009 71.7 9.43 55.6 1.94 67.7 0.025 25.4 0.24
392548077255102 AMW 568-13B FR Dd 239 16 2,180 7.9 0.2 259 0.04 0.01 108 51.8 250 1.66 498 0.117 134 0.27
392605077263001 BMW71C FR Dd 240 14.9 437 8.2 5.5 190 4.24 0.018 45.3 26.8 4.6 1.37 26.2 0.022 8.39 0.04
392617077275801 FR Dd 243 (PW) FR Dd 243 22.3 175 7.1 2.9 67.5 1.69 0.007 24.6 2.56 6.05 0.82 7.99 0.03 9.52 0.05
392625077270301 BMW5 FR Dd 241 13.6 515 7.6 3.2 244 1.49 0.05 68.1 19.4 10 2.28 21.4 0.048 19 0.04
392625077270302 BMW5D FR Dd 242 13.9 531 7.8 3.4 305 1.53 0.026 73.5 22.5 8.57 1.85 16.7 0.044 18.7 0.03
392553077262601 Site 5/FR Dd 244 FR Dd 244 13.3 493 7.8 4.8 214 4.74 0.033 65.5 16.6 12.4 1.64 33.2 0.021 9.48 0.03
392605077262801 South Armory Spring FR Dd 245 13.9 513 7.7 5.7 237 5.39 0.043 68.5 17.4 12.8 1.6 35.2 0.02 10 0.04
392604077274201 BMW47D FR Dd 246 14.8 553 7.8 2.6 274 4.05 0.027 66.4 26.1 12.5 1.08 40.7 0.027 3.57 0.02
Table 11.    Select groundwater field parameters, major ion concentrations, trace metal concentrations, and redox classification, Fort Detrick Area B, Maryland.[—Left]

Table 11.    

Select groundwater field parameters, major ion concentrations, trace metal concentrations, and redox classification, Fort Detrick Area B, Maryland.[—Right]

[USGS, U.S. Geological Survey; ID, identification; °C, degrees Celsius; µS/cm, microseimens per centimeter at 25 °C; DO, dissolved oxygen; mg/L, milligrams per liter; N, nitrogen; P, phosphorus; SiO2, silica dioxide; µg/L, micrograms per liter; <, less than; Cl:Br, mass ratio of chloride to bromide; —, no note recorded]

USGS site number Sample name Site ID Silica as SiO2 (mg/L) Manganese (µg/L) Iron (µg/L) Barium (µg/L) Boron (µg/L) Strontium (µg/L) Vanadium (µg/L) Lithium (µg/L) Selenium (µg/L) Uranium (µg/L) Redox category Na:Cl molar ratio Cl:Br mass ratio Cl:Br note
392556077263301 RISP-3 FR Dd 213 11.9 <0.4 <10 61.5 11 69.6 0.3 0.78 0.17 0.119