Purpose and Scope

The primary purpose of this report is to present results of hydrogeologic data collection for 2019–22 and to document the development of a regional groundwater-flow model for 2000–19. The developed model was also used to simulate potential future impacts to the groundwater system in Wake County through 2070. Specifically, this report aims to enhance the understanding of local groundwater resources, help identify critical hydrogeologic parameters using a developed modeling tool, and assess the groundwater system response to future climate variability and urban development. This report presents the description of the development and calibration of a groundwater-flow model, a discussion of both historical and future simulation results, and the limitations of the model.

The scope of work included compilation of existing data, such as historical water levels and groundwater use, and collection of new data to characterize hydraulic properties and geochemical conditions of the fractured-rock groundwater system. Results from 17 monitoring well slug tests and 21 vertical borehole-flow measurements were used to estimate bulk aquifer hydraulic properties. Water samples for 19 sites were collected for laboratory analysis of naturally emplaced stable isotopes of water, dissolved gases, and compounds used as tracers for age dating of groundwater. Analytical results were used to describe the source composition of local precipitation, groundwater, and surface water, as well as identify potential groundwater mixing along flow paths. Building on the established regional hydrogeologic framework, this study integrated local data from Wake County that culminated in a groundwater-flow model for the greater Wake County area to assess current groundwater conditions. As part of this work, a data release with the groundwater model and associated compiled datasets was made publicly available. The model also served as a tool that evaluated current groundwater conditions and how the sensitive components of the groundwater system might respond to various future scenarios.

Description of the Study Area

Although the area of study for this report is primarily focused on Wake County (referred to hereinafter as the “study area”) (fig. 1), the model area used in this study extends into portions of neighboring counties (Durham, Orange, Alamance, Chatham, Harnett, Johnston, Nash, Franklin, Granville, and Person Counties) (fig. 2) to capture regional hydrologic influences and groundwater flow into and out of the Wake County study area. The Wake County study area has a total area of about 857 square miles (mi²), representing about 26 percent of the entire model area (3,325 mi²). Twelve municipalities, including the State capital City of Raleigh, are located within Wake County: Morrisville, Cary, Apex, Holly Springs, Fuquay-Varina, Garner, Knightdale, Wendell, Zebulon, Rolesville, and Wake Forest (fig. 1). As of 2020, the population of Wake County exceeded 1 million people (1,129,393; U.S. Census Bureau, 2020).

Wake County is located predominantly within the Piedmont physiographic province (Fenneman, 1938) (fig. 1). A small section of southeastern Wake County (about 35 mi²) extends across the fall line into the Coastal Plain physiographic province. The county is characterized by a hill and valley topography with rounded slopes, where land elevation ranges from 137 to 558 feet (ft) above the North American Vertical Datum of 1988 (NAVD 88) (fig. 2). Wake County has a humid subtropical climate with daily-mean air temperatures ranging from 47 to 71 degrees Fahrenheit (°F) and annual precipitation of about 46 inches (National Oceanic and Atmospheric Administration, 2020). Based on 2019 data, 42 percent of land cover in Wake County is classified as developed, while only 23.3 percent of the entire model area is classified as developed (fig. 3; table 1) (Dewitz and U.S. Geological Survey, 2021). Other predominant land cover classes in the county study area include forest (36 percent), pasture/hay (6.5 percent), cultivated crops (4.6 percent), and wetlands (4.3 percent).

Wake County lies within the Neuse River and Cape Fear River Basins (fig. 1). More than 80 percent of the county is within the Neuse River Basin with major streams that include Crabtree Creek, Walnut Creek, and Middle Creek. Several Wake County lakes are within the Neuse River Basin including Falls Lake, which serves as one of the primary sources of water supply for the City of Raleigh, as well as smaller lakes such as Lake Crabtree, Lake Wheeler, and Lake Benson. The Cape Fear River Basin in the southwestern corner of the county includes B. Everett Jordan Lake and Shearon Harris Reservoir. The USGS operates a streamgage network within both basins that includes four continuous streamgage sites used for hydrologic analysis in the model area (figs. 1 and 2).

Compilation of USGS 5-year national water-use estimates for the years 1995–2015 show that total annual surface water and groundwater use in Wake County ranged from about 56.3 to 88.9 Mgal/d (fig. 4A) (Solley and others, 1998; Hutson and others, 2004; Kenny and others, 2009; Maupin and others, 2014; Dieter and Maupin, 2017). Surface water sources provided about 40.4 Mgal/d in 1995 and about 72.8 Mgal/d in 2005, throughout the county. Total annual groundwater withdrawals in Wake County ranged from about 15 Mgal/d in 2000 to about 19.8 Mgal/d in 2010, composing about 20–30 percent of all water use in the county (fig. 4A). Groundwater withdrawals were highest for domestic self-supply followed by public supply, which together were increasing through 2010 and then remained nearly constant through 2015 (fig. 4B). The USGS data include an estimated withdrawal coefficient of per capita use from domestic wells in Wake County of 70 gallons per day (Dieter and Maupin, 2017). Previous studies have assumed that most groundwater withdrawn from domestic wells is returned to the aquifer if the residence has an onsite wastewater (septic) system, given some consumptive loss to evapotranspiration (Daniel and Harned, 1998; Horn and others, 2008; Goode and Senior, 2020). There does remain uncertainty regarding the fate of septic discharge through the shallow subsurface to recharge deeper bedrock.

Hydrogeologic Setting and Conceptual Model of Flow System

The Wake County study area is located at the eastern edge of the Piedmont physiographic province and is underlain by regional geologic formations that generally trend northeast to southwest (Horton and Zullo, 1991). The rock types include stratified sedimentary basin rock, crystalline rock with varying degrees of metamorphism, and younger unconsolidated coastal plain sediments that overlie fractured bedrock eastward toward the Piedmont and Coastal Plain boundary (figs. 1 and 5). Lithologic units defined across North Carolina were compiled into distinct hydrogeologic units for the Piedmont and Blue Ridge physiographic provinces by Daniel (1989) with an associated map by Daniel and Payne (1990) (fig. 5). Hydrogeologic unit classification is associated with the water-bearing potential of the rock types, where primary and secondary porosity is related to composition, texture, and weathering rates of rock. The study area is underlain by massive and foliated crystalline rocks and sedimentary basin rocks. A large part of Wake County primarily is composed of closely folded metasedimentary and metavolcanic rocks, including phyllite, graphite-bearing mica schists, biotite-hornblende gneiss, quartzite, and amphibolite (Parker, 1979). Coarse-grained granitic plutons are found within the study area (fig. 5, IFI unit), including one of the largest in the Eastern United States that stretches across eastern Wake County (Farrar, 1985). Unmetamorphosed igneous intrusive rocks are present as diabase dikes and granitic veins. The Durham subbasin in the Triassic Deep River Basin is located on the western edge of Wake County (fig. 5, TRI unit) and is composed of stratified fanglomerate and red silty sandstone to mudstone sourced from the erosion of uplifted crystalline rocks after Mesozoic rift faulting (Parker, 1979; Bain and Brown, 1981). Triassic sedimentary rocks are bounded to the east by the north-northeast trending Jonesboro Fault (Campbell and Kimball, 1923), separating the basin from the metaigneous and metavolcanic rock types to the east (May and Thomas, 1968; Parker, 1979) (fig. 5). Onlapping coastal plain sediments of gray to white sand with interbedded clay lenses lie nonconformably atop metamorphic and igneous rock in the southeastern part of the county (May and Thomas, 1968).

A conceptual model of groundwater flow within the Piedmont of North Carolina was described by LeGrand (1967), Heath (1980, 1984), and LeGrand and Nelson (2004) and was adapted for this study. Groundwater flow generally follows the hydraulic gradient, moving from topographic highs to lows and forming local flow systems within small drainage basins that exhibit some topographic relief. Water-table levels are typically highest near the highest points of land-surface altitude near drainage basin boundaries, such as hills and ridges, whereas the lowest water-levels occur near streams. Recharge to the hydrogeologic units originates from precipitation that infiltrates, through the soil and root zone, down to the water table, continuously moving downgradient toward streams until discharge. Some groundwater may be lost to evapotranspiration as it flows near the surface, before being discharged as seepage to gaining streams. In areas with minimal topographic relief, these local flow systems are less distinct, where slower moving intermediate and regional groundwater flow can be predominant (Tóth, 1963).

Groundwater occurs in two primary layers of the flow system: shallow regolith and deeper fractured bedrock. The regolith consists of soil, local alluvium, residuum (sandy clay from the weathering of feldspar and mica minerals), and in situ weathered bedrock, also known as saprolite. Saprolite may retain some relict structural features of the parent bedrock, such as foliation and filled fractures. Regolith thickness varies geographically and is influenced by factors such as parent rock, topographic setting, and geologic history (LeGrand and Nelson, 2004). In the Piedmont, regolith typically is thinnest near streams and thicker on slopes, with variable thickness along hilltops and ridges. The lower regolith may contain a transition zone from saprolite to bedrock as weathering decreases with depth, often less gradually and more abruptly in highly foliated parent rock (Harned and Daniel, 1992). Where present, the transition zone has been described as more permeable than the upper regolith and soil layers, a condition likely occurring because intermediate weathering has not progressed to the formation of clay (Stewart, 1962; Nutter and Otton, 1969; Daniel and Harned, 1998). Groundwater within the regolith layer flows both vertically and laterally through intergranular spaces before discharging to nearby streams. Lateral flow primarily is driven by the rolling topography of the landscape that creates short flow paths in local flow systems. This lateral flow can be enhanced by the presence of clay horizons and relict structural features in the regolith, which impede some vertical flow. At the base of the regolith, the partially weathered material of a transition zone can provide increased permeability and act as a conduit for lateral groundwater movement. Some groundwater in the transition zone flows downward into the underlying bedrock, moving through a network of connected permeable fractures before discharging upwards into downgradient streams (Heath, 1984; Daniel and Harned, 1998).

Groundwater flow in the fractured bedrock depends on the hydraulic connection with the base of the regolith in recharge and discharge areas, as well as the connectivity of fractures. Most groundwater storage occurs within the overlying regolith, as the fractured bedrock has little to no storage capacity because of its relatively impermeable rock matrix (Heath, 1980). In the Piedmont, fracture features primarily form by fault displacement, shear (typically occurring in zones), and stress relief that can result in orthogonal vertical fracture sets and nearly horizontal sets, depending on the rock fabric. The various mechanisms of fracture development and nonuniform distribution of fracture networks contribute to the heterogeneity in aquifer properties of the bedrock and provide preferential flow paths that can cause groundwater flow to deviate from the typical flow paths observed in more homogeneous, porous media. These flow pathways within the connected fracture networks are largely part of local flow systems, yet some pathways may be part of broader flow systems that cross small drainage basin boundaries. Most of the groundwater flow in the Piedmont is constrained by increasing lithostatic pressure at depth, limiting most groundwater circulation primarily to the upper 800 ft of bedrock, with the most productive fractures within 350 ft from the top of bedrock (Daniel, 1989; Daniel and others, 1997; LeGrand and Nelson, 2004).

Previous Investigations

Previous hydrogeologic work in the study area was summarized in Antolino and Gurley (2022), which is adapted for parts of this section. The summary by Antolino and Gurley included a detailed groundwater study in Wake County by May and Thomas (1968), who identified correlations between well yields, rock types, and topographical locations. Welby and Wilson (1982) examined the connection between geological factors and well yields, noting that an increase in the proportion of impervious surface from future development may decrease groundwater recharge and result in local water-supply problems. A 2003 groundwater study by CDM (2003) estimated water budgets across major drainage basins in the study area and highlighted the need for information that would help to quantify local effects of increased development on groundwater sources. Chapman and others (2011) investigated declining groundwater levels in northern Wake County following reports of dry domestic wells in 2005 and 2007. Findings from the study revealed substantial aquifer storage loss from increased withdrawals for domestic irrigation and CWS supply, as well as identified groundwater-level data trends that could be correlated with local dominant geological structural features.

Antolino and Gurley (2022) compiled and analyzed well depth, casing depth, and reported well yield for well sites in Wake County. They collected geophysical logs and determined the orientation of geologic structure in 15 bedrock wells in different hydrogeologic units in the county. They used the Soil-Water-Balance (SWB) model (Westenbroek and others, 2010, 2018) to calculate the spatial and temporal variation of water budget components in the greater Wake County area from 1981 to 2019, including net infiltration below the root zone (or potential groundwater recharge). Estimates of mean-annual recharge determined with the SWB model ranged from 5.6 to 9.2 inches across the gaged basins. These estimates of potential recharge were compared with stream base flows derived from hydrograph separation for drainage basins with USGS streamgages, within the model area. The comparison showed that the model produced recharge values generally within 2 inches of base flow estimates for most basins. The SWB model also provided estimates of daily evaporation rates that ranged from 0.04 to 0.13 inches per day and daily runoff rates that ranged from 0.02 to 0.05 inches across the gaged basins. The monthly mean-daily rates for the major inputs and outputs determined for the SWB model area by Antolino (2022) are shown in figure 6.

Chapman and others (2005, 2007) studied the fractured-rock hydrogeology at the Lake Wheeler Road Field Laboratory, situated at a North Carolina State University research farm, including groundwater-quality analysis. Analytical results indicated that groundwater nitrate concentrations were highest in the shallower regolith and transition wells, whereas concentrations of fluoride, calcium, sulfate, molybdenum, uranium, radon, and arsenic were highest in the deeper bedrock wells. Using various groundwater age-dating tracers, the study identified longer flow paths with greater residence times in groundwater discharge areas and estimated that modern groundwater in the regolith is younger than deeper groundwater in the bedrock by several decades. McSwain and others (2013) conducted a study in the vicinity of the Neuse River Resource Recovery Facility (located near WK–334, fig. 5) to examine potential effects of historical use of biosolids applications in nearby fields on water quality in the subsurface and in the Neuse River. High nitrate levels were found in both the shallow regolith and deeper fractured bedrock layers of the groundwater system, highlighting the direct influence of land-use activities on groundwater in the area. The study also documented substantial nitrate flux from groundwater at the site to the adjacent Neuse River, indicating important interactions between the local groundwater and surface water systems.

As part of the Appalachian Valleys and Piedmont Regional Aquifer-System Analysis study, Daniel and others (1997) developed an interpretive steady-state numerical model for the Indian Creek Basin in the southwestern North Carolina Piedmont to understand the relation between groundwater-flow components and streams within the fractured crystalline bedrock groundwater system. The multilayered model was developed assuming an equivalent porous media for bulk fractured-rock aquifer properties and provided insights into the regional groundwater system. Most of the simulated groundwater flow was along short flow paths within the shallowest layers to streams. The model simulation showed that groundwater flow decreased with depth by two orders of magnitude between the top and bottom model layers, about 700 ft below land surface. Particle tracking simulated movement of groundwater flow by using numerous particles distributed across the top model layer to indicate groundwater travel times through regolith to the stream that ranged from 10 to 20 years at horizontal distances less than 1,000 and 2,000 ft, respectively. Travel times within the bedrock increased with depth. For horizontal distances between 2,000 and 5,000 ft, flow paths in the upper 100 ft of bedrock had travel times that ranged from 20 to 90 years, whereas the few flow paths in the lower bedrock (675 ft below land surface) had travel times that exceeded 300 years.

Development-related future water demand across North Carolina, including Wake County, was projected from 2012 to 2065 by Sanchez and others (2018, 2020a). Water demand was quantified as the sum of public supply, industrial self-supply, and domestic self-supply uses from both surface water and groundwater sources obtained from the 2010 county-level USGS national water-use estimates (Maupin and others, 2014). To obtain locally relevant representation of the spatial distribution of water use, county-level water-use estimates were spatially disaggregated to census tract units by using a population weighted procedure. This modeling framework evaluated the joint effects of social and environmental conditions on urban water demand under different future climate scenarios and urban growth. To evaluate the effects of growing populations, expanding development, and temperature change, the modeling framework assumed that all other predictors remained at 2010 conditions. Probabilistic model projections from Sanchez and others (2020a, 2020b) indicated increasing trends for future development-related water demand across Wake County census tracts, with upper limit estimates showing an increase of 78 percent to more than 500 percent by 2065 (relative to 2010 estimates). Within Wake County, the forecast percentage change was highest in southeast and north Raleigh and urban peripheries of Morrisville, Cary, and Apex (fig. 1), which are all areas expected to experience rapid low-density development in coming decades (Terando and others, 2014; Sanchez and others, 2020a).

Well-Construction and Groundwater-Level Data

Datasets on well construction and groundwater levels were compiled for wells located throughout the model area. This study primarily relied on existing well-construction and groundwater-level data that were retrieved from the USGS National Water Information System (NWIS) database (U.S. Geological Survey, 2022) and the North Carolina Department of Environmental Quality Division of Water Resources (NCDEQ DWR) monitoring well database (North Carolina Department of Environmental Quality, Division of Water Resources, 2022). These well-construction and groundwater-level data were primarily used for groundwater-flow model development and calibration. Information on well-casing depths was retrieved from the USGS NWIS database (U.S. Geological Survey, 2022) for more than 2,800 wells located in the model area and used to estimate the top of the bedrock layer in the model simulation. Both historical and recent groundwater levels that were available for the period 1941–2021 from the USGS NWIS and NCDEQ DWR databases were used for calibrating the groundwater-flow model.

Discrete water-level measurements were made using a chalked graduated steel tape or electric water-level tape, consistent with techniques described in Cunningham and Schalk (2011), and yielded data having an accuracy within 0.01 ft. The water-level tapes are routinely checked for quality assurance by USGS field staff and the USGS Hydrologic Instrumentation Facility based on USGS policy and techniques (U.S. Geological Survey, 2015b; Fulford and Clayton, 2015). Water-level measurements were made from a measuring point established at both discrete and continuous sites. In instances where obstructions within the well (such as submersible pumps) prevented the use of established measurement methods, sonic water-level meters were utilized, providing a measurement accuracy within 0.1 ft. Prior to sonic water-level measurements, air temperature within the well was recorded using a thermistor to configure the meter to account for variations in sound velocity between wells (Ravensgate Corporation, 2025).

For this study, two synoptic water-level surveys were completed in and adjacent to Wake County to provide spatial snapshots of the groundwater system across different seasons. This approach provides information on the influence of seasonal variations on groundwater levels. The first water-level survey included discrete water-level measurements for 102 domestic and monitoring wells during November and December 2020, when evapotranspiration rates are typically lower (fig. 6). The second water-level survey included measurements for 100 of mostly the same domestic and monitoring wells during June and July 2021, when evapotranspiration rates are typically highest (fig. 6).

Continuous groundwater-level measurements also were collected to improve spatial and temporal records of groundwater-level data for Wake County and to serve as observations for the groundwater model. As described by Antolino and Gurley (2022), a network of 17 USGS wells was established in 2019 to monitor continuous groundwater levels in the fractured-rock groundwater system for Wake County. The USGS monitoring-network wells (4 regolith and 13 bedrock), as well as the additional two NCDEQ DWR wells used for project data collection, are shown in figure 5 and summarized in table 2. The NCDEQ DWR operated and maintained continuous data-collection activities for two well sites in the study area, 353320078483101 (WK–438) and 354703078424401 (WK–439), which correspond to NCDEQ DWR sites N41G3 and K40M1, respectively (table 2). Water-level data for sites WK–438 (N41G3) and WK–439 (K40M1) are available through the NCDEQ DWR monitoring well database (North Carolina Department of Environmental Quality, Division of Water Resources, 2022).

Continuous water levels were measured at 15-minute intervals in each well by using a vented submersible pressure transducer (Cunningham and Schalk, 2011). USGS data are recorded using a data-collection platform and are relayed hourly through satellite transmission for upload to the USGS NWIS database (U.S. Geological Survey, 2022). The water-level monitoring sites are inspected on a routine basis, and the recorded water-level data are field verified with discrete steel-tape measurements approximately every 6 to 8 weeks.

Aquifer Testing

Slug tests and borehole-flowmeter logging were conducted in selected monitoring wells to estimate bulk hydraulic properties of the fractured-rock groundwater system. Measurement and analysis of groundwater-level response to induced stress, along with the identification of vertical flow zones within the regolith and bedrock wells, provided data used to estimate values of transmissivity (T) and horizontal hydraulic conductivity (K_h) within the fractured-rock groundwater system.

Water-Quality Sampling

The primary focus for collecting water-quality samples for laboratory analysis was to better understand the geochemical and hydrologic conditions of the groundwater-flow system. Groundwater samples were collected during November and December 2019 from all 17 USGS wells in the monitoring network and WK–439 (NCDEQ DWR well K40M1), with the exception of WK–438 (NCDEQ DWR well N41G3), for a total of 14 bedrock wells and 4 regolith wells (fig. 5; table 2). Water temperature, specific conductance, pH, and dissolved oxygen (DO) were monitored in the field while the well was being purged of at least one well volume of stagnant water. Water-quality multiparameter sondes were calibrated each sampling day and used to monitor field parameters prior to sampling. Groundwater samples were collected after stabilization of these physical and chemical properties (U.S. Geological Survey, 2006). Wells were purged at pump rates ranging from 0.25 to 4 gal/min depending on well yield and water-level response. Submersible pumps were used to purge and collect samples at all wells except WK–436, which was sampled using a double-ball bailer after removing three well volumes of water. Field equipment was cleaned between sampling sites by following established methods (U.S. Geological Survey, 2006). Groundwater samples primarily were collected for analysis of stable isotopes of water, dissolved gases, and age-dating tracers.

Samples of surface water and precipitation also were collected for analysis of stable isotopes of water. These samples were collected at USGS monitoring station 0208735012 (Rocky Branch below Pullen Road at Raleigh, N.C.), which has a continuous streamgage site and colocated rain gage (fig. 5) (U.S. Geological Survey, 2022). Eight surface water grab samples from the midpoint of the stream were intermittently collected from Rocky Branch between April 2021 and January 2022. A total of 48 precipitation samples were collected about every 2–4 weeks from December 2019 to April 2022 by using a commercial version (Palmex Ltd., Zagreb, Croatia) of the rain collector described in Gröning and others (2012). The rain collector was installed next to the Rocky Branch rain gage to relate the accumulated rainfall amount in the composite precipitation sample. The rain collector is designed to thermally isolate the sample from sunlight and limit exposure to the atmosphere to minimize fractionation effects from evaporation. Precipitation samples that were collected and exposed to freezing conditions, that may have potentially contributed to potential isotopic fractionation effects, were excluded from analysis.

All water-quality samples were collected and processed using techniques described in the “National Field Manual for the Collection of Water-Quality Data” (U.S. Geological Survey, 2006). Analytical results for all water-quality samples collected during the study are available from the USGS NWIS database (U.S. Geological Survey, 2022).

Groundwater-Flow Model Development

Existing and new hydrologic data were integrated to develop a groundwater-flow model to simulate historical groundwater conditions and to forecast future groundwater availability in the Wake County study area. The model was constructed using the USGS modular three-dimensional finite-difference groundwater-flow model code MODFLOW–2005 (Harbaugh, 2005), specifically the Newton formulation version, MODFLOW–NWT (Niswonger and others, 2011). MODFLOW simulates an aquifer as porous media, meaning that water is assumed to move through pore space between granular materials. However, in the study area, groundwater within the fractured bedrock flows through fracture networks of discrete transmissive conduits within a relatively impermeable rock matrix. The application of MODFLOW for this study assumes that flow through the fractured-rock groundwater system is at the scale of miles and can be represented on a regional scale as a single-continuum porous-equivalent medium. This modeling approach using bulk effective aquifer properties has previously been used with discrete conduit groundwater flow in other groundwater systems (Tiedeman and others, 1997; McCoy and others, 2015; Kuniansky, 2016; Harte, 2021).

Physical Properties and Chemical Constituents

Results for physical properties and chemical constituents (water temperature, pH, specific conductance, and DO) measured during groundwater sample collection are summarized in table 7. Water temperatures ranged from 15.1 to 19.0 degrees Celsius (°C) in the regolith and 16.4 to 24.3 °C in the bedrock. Values for pH ranged from 4.8 to 5.7 in the regolith and 5.3 to 7.7 in the bedrock. Results for the four sets of paired wells tended to indicate higher values of pH (more basic) and specific conductance in the deeper bedrock wells as compared to the shallower regolith wells (table 7). Specific conductance ranged from 72 to 367 microsiemens per centimeter at 25 °C (µS/cm) in the regolith and 89 to 598 µS/cm in the bedrock. The highest specific conductance measurements were likely related to land use near the well site, namely agricultural practices at the WK–283 site (598 µS/cm) at the North Carolina State University Lake Wheeler Road Field Laboratory and land application of biosolids at the WK–334 site (549 µS/cm) at the Neuse River Resource Recovery Facility. Concentrations of DO ranged from 2.1 to 6.9 milligrams per liter (mg/L) in the regolith and 0.15 to 8.1 mg/L in the bedrock wells. The highest observed DO concentration of 8.1 mg/L was measured in bedrock well WK–368, caused by water from shallower fractures cascading down the well borehole because of depressed water levels from nearby pumping--cascading water along the borehole wall was confirmed by a well camera survey. Six of the bedrock wells had DO concentrations less than 1 mg/L, indicating more reduced geochemical conditions for these sites.

Stable Isotopes of Water

The results of the analysis of stable isotopes of water in precipitation, surface water, and groundwater samples are listed in table 8. The relation between the δ¹⁸O and δ²H isotopic composition in precipitation samples is plotted in figure 9A. The isotope values for the precipitation samples (n = 47) ranged from −10.05 to −1.92 and −67.20 to 4.25 per mil for δ¹⁸O and δ²H, respectively. The samples did not show strong seasonal variation, but those samples collected during spring (March–May) had the most variation in isotopic composition. This range in isotopic concentration reflects the dynamic seasonal transition from drier northern air masses to more humid air masses originating from the south bringing fluctuating temperatures with a mixture of precipitation events across the study area. For example, the March 17, 2021, precipitation sample slightly deviates from the general trend of the data that may reflect a cooler, drier continental air mass as the moisture source. An LMWL of

was determined for the Wake County study area by linear regression (R² = 0.93) by using precipitation data collected from the rain gage located at USGS monitoring station 0208735012 (Rocky Branch below Pullen Road at Raleigh, N.C.) (fig. 9A; U.S. Geological Survey, 2022). The LMWL has a slightly lower slope of 7.77 compared to the GMWL slope of 8.0 defined by Craig (1961). This difference in slope with the GWML likely reflects the specific climatic conditions and moisture sources of the Raleigh area in contrast to the broader range of global conditions captured in the GMWL. The LMWL has a steeper slope than the river water line developed for the State of North Carolina using isotopic data from river sites by Kendall and Coplen (2001), δ²H = 6.32 * δ¹⁸O + 2.9 ‰ (n = 87). The river water line, acting as a proxy for the isotopic composition of modern precipitation, is based on the assumption that most of the surface water represented base flow and recent precipitation, yet the shallower slope likely reflects some influence of evaporation effects in the water samples, possibly due to mixing with older, evaporated waters in the basins.

The δ²H and δ¹⁸O results for groundwater generally plotted along and slightly below the LMWL with a regression line with a slope equation of δ²H = 6.91 * δ¹⁸O + 7.49 ‰ (R² = 0.84) (fig. 9B). This pattern suggests that the groundwater is isotopically similar to modern local precipitation, indicating that the recharge from rainfall has undergone minimal evaporation before entering the aquifer. This isotopic signature highlights that local precipitation is the predominant source for all of the sampled wells across the county, implying that any changes in rainfall directly impact recharge to the groundwater system. Ranges for δ²H and δ¹⁸O in the groundwater samples were from −32.67 to −26.82 per mil and −5.88 to −5.04 per mil, respectively. Results for all groundwater samples fell within a relatively narrow isotopic range, reflecting an averaged isotopic composition where short-term and seasonal variations from similar precipitation events across the county have been smoothed out and integrated into a composite signature.

The isotopic results for the Rocky Branch stream samples plotted along a regression line with a slope equation of δ²H = 6.91 * δ¹⁸O + 6.04 ‰ (R² = 0.98) (fig. 9B). Except for the January 2022 sample, all stream samples were collected during low-flow conditions (less than 1 cubic foot per second) and generally plotted in the same isotopic range as the groundwater samples, indicating that groundwater discharge is the predominant source of streamflow during low-flow conditions. The larger spread in the isotopic concentrations for the stream samples likely reflects the influence of seasonal evaporation during warmer periods, as well as variability with water-source mixing of recent rainfall and groundwater discharge.

Groundwater Dissolved Gases

Based on field meter readings of DO in the bedrock wells, 10 wells were considered to have oxic groundwater conditions (greater than 0.5 mg/L), and 4 wells were considered to have anoxic groundwater conditions (less than 0.5 mg/L) (McMahon and Chapelle, 2008) (table 7). Most of the sampled wells having oxic conditions suggest some degree of connection with more oxygenated sources, possibly through connected fracture networks at the top of bedrock with shallow groundwater in the regolith. Anoxic conditions may be present in deeper fracture zones where oxygen is consumed along the flow path. The laboratory results of the concentrations of the dissolved gases methane, carbon dioxide, nitrogen, oxygen, and argon in the 14 bedrock groundwater samples are provided in table 9. Three of the well sites (WK–428, WK–429, and WK–431) had laboratory-analyzed dissolved-gas samples with reported DO concentrations greater than 0.5 mg/L, whereas the remaining 11 well sites had anoxic concentrations of reported DO around 0.1 mg/L in the samples (table 9). This discrepancy between the field and laboratory results suggests potential changes in DO within the sample during storage. Laboratory analyses of DO can often be biased toward lower concentrations because of rapid microbial oxygen consumption during sample storage, despite proper preservation techniques (U.S. Geological Survey, 2006). Therefore, field measurements of DO generally are considered more reliable when interpreting groundwater conditions. The other dissolved gases measured in this study are less reactive and, if preserved properly, maintain stability prior to analysis.

Trace levels of dissolved methane were measured in the samples from 6 of 14 bedrock wells, and relatively higher concentrations were measured in samples from wells WK–429 and WK–439 (table 9). The presence of methane is suggestive of reducing conditions within the fractured-rock groundwater system and aligns with anoxic levels seen in well WK–439. However, the groundwater sample from well WK–429 is oxic and has the highest levels of dissolved methane. As previously discussed in the section “Characterization of Aquifer Hydraulic Properties,” the largest source of inflow to this well occurs from a fracture at the top of bedrock below the well casing, possibly providing a source of more oxygenated groundwater. This inflow may be mixing with anoxic water where methane production may be occurring in more stagnant zones deeper in the well. This mixing may have been enhanced, or even induced in the case where the natural flow regime is reversed, during sample collection, by pumping in wells that have fractures within different redox zones. Additional geochemical information outside the scope of this study would be needed, such as hydrogen sulfide and acetate concentrations, to verify the occurrence of methanogenic conditions in the wells.

Dissolved carbon dioxide concentrations ranged across two orders of magnitude from 0.8219 mg/L in well WK–283 to 95.3804 mg/L in well WK–426 (table 9). The values near the lower end of this range are closer to atmospheric equilibrium, and the larger values suggest more direct connection with shallow groundwater that has infiltrated soil zones that are high in organic matter, where microbial activity is elevated. Additionally, in more urban and industrial areas, it is possible that these levels may be influenced by the aerobic degradation of contaminants, such as volatile organic compounds (Appelo and Postma, 2005; Clark, 2015).

Ranges in dissolved nitrogen and argon concentrations in groundwater samples were from 14.0851 to 33.8943 mg/L and 0.5170 to 0.8596 mg/L, respectively (table 9). In this study, average excess air values ranged from 0.3 to 11.5 cubic centimeters (at standard temperature and pressure per liter) with the highest levels in well WK–283. Because nitrogen can also be produced from denitrification processes in the subsurface, excess nitrogen within samples from the wells was assessed by comparing observed concentration to expected values for nitrogen based on the air-water equilibrium derived from argon (Böhlke, 2002). Estimates of excess nitrogen present in wells WK–430, WK–283, and WK–334 suggest active denitrification processes in the aquifer near these wells. Laboratory corrections were applied for the presence of excess nitrogen above expected solubility for the nitrogen values reported for these three wells in table 9.

Accounting for excess air in groundwater samples refines the temperature-dependent gas solubility, leading to more accurate estimates of recharge temperature and aiding in determining when groundwater was recharged. Recharge temperatures calculated by the USGS RGDL from the dissolved-gas concentrations of nitrogen and argon are listed in table 9. The average recharge temperature was 15.2 °C, consistent with the average annual temperature for the study area (National Oceanic and Atmospheric Administration, 2020). The calculated recharge temperature for well WK–368 was higher than those of other wells, likely influenced by significant drawdown from nearby supply well pumping, which may have altered local flow dynamics to draw more groundwater from shallower, recently recharged fractures. These shallow fractures likely have shorter flow paths, reflecting recent recharge from prior events before sampling, likely during the warmer summer months. Deeper fractures, typically with longer flow paths, may contain a mixture of groundwater from multiple seasons. Cascading water observed along the open borehole during a borehole camera survey suggests a shallow productive fracture zone below the bottom of casing contributing more recently recharged groundwater. Refined recharge temperature estimates improve the reliability of dissolved-gas and tracer concentrations used in groundwater age dating (Stute and Schlosser, 2000).

Groundwater Age Dates

Concentrations of CFC-11, CFC-12, and CFC-113 measured in samples from the 14 bedrock wells were used by the USGS RGDL to calculate apparent groundwater ages based on the assumption of piston-type flow within the aquifer. The apparent ages based in the CFC analyses ranged from the 1950s to the 1990s (table 10). The groundwater age for a sample is a function of recharge and groundwater flow rates, as well as source and flow path mixing. As with all dissolved constituents in groundwater, the CFC concentrations can be influenced by flow transport processes, including the mixing of older and younger waters from different flow paths. Because groundwater samples were pumped from deep, open well boreholes, each sample likely contains some mixture of groundwater having different recharge ages.

Assuming there is a binary mixing of older CFC-free water and younger water containing CFCs, the tracer age date will be based on the younger water fraction within the sample. The younger and older fractions within the binary mixture can be assessed from ratios of the CFCs. Based on laboratory evaluations of information derived from the CFC tracer ratios, the mixing of younger, more modern recharge water with older groundwater could be derived for 5 of the 14 bedrock wells (table 10). Both samples for wells WK–435 and CH–252 contained less than 25 percent of young water. Mixing of the sample for well WK–430 reflected nearly equal parts of older and younger water. Younger water made up about 75 percent of the sample from well WK–434 and nearly all the sample from well WK–431.

Plummer and Busenberg (2000) suggested that CFC-11 may be, at least partially, degraded in groundwater having DO concentrations less than 0.5 mg/L. Four of the 14 bedrock groundwater samples had field-measured DO concentrations less than or equal to 0.3 mg/L (table 7). Other studies have shown that CFC-11 is more sensitive to anerobic microbial degradation than CFC-113 and CFC-12, especially once redox conditions reach sulfate reduction (Semprini and others, 1990; Katz and others, 1995; Shapiro and others, 1997). Five of the wells that had trace concentrations of dissolved methane also had oxic levels of DO, which suggests the mobilization of methane from groundwater in deeper fracture zones mixing with more oxic water in shallower zones (tables, 7, 9, 10). All CFCs would have likely undergone some degradation in low sulfate, methanogenic environments (Oster and others, 1996; Deipser and Stegmann, 1997; Plummer and Busenberg, 2000). This degradation would result in a bias toward older recharge ages for samples with a large fraction of highly reduced, methanogenic groundwater. Higher levels of methane were detected in two groundwater samples, from wells WK–429 and WK–439 (table 9), but oxic levels were such that the effect on the groundwater age is minimal (table 7).

Concentrations of all three CFCs measured in groundwater samples from well WK–428 were more than the air-water equilibrium and were likely influenced by an unknown anthropogenic source. In addition to microbial degradation, contamination will also increase the uncertainty in CFC concentrations. The CFC-derived groundwater ages were examined in combination with results from the additional environmental tracers of ³H, ³He, ⁴He, and neon to reduce uncertainty in groundwater age-dating results.

The results for ³H, ³He, ⁴He, and neon are listed in table 11. Concentrations of ³H measured in the groundwater samples ranged from 0.05 to 3.92 TU. The samples contained tritiogenic ³He concentrations ranging from less than 1 to 19.2 TU. Of the 14 bedrock wells sampled, 5 wells had samples with very low concentrations of tritiogenic ³He. All samples contained measurable amounts of terrigenic helium, with most samples showing some mixing with tritiogenic-sourced helium. Apparent ages based on the ³H and ³He concentrations ranged from the mid-1930s to the late 2000s, reflecting a broader range of apparent age dates for groundwater within the fractured-rock groundwater system compared to the CFC tracer data (tables 10 and 11).

Conceptually, groundwater moves from porous flow in the regolith to complex fracture flow in the bedrock; therefore, a binary-mixing LPM was initially used to analyze both CFC and ³H/³He tracer data. The binary-mixing model was unsuccessful at producing optimized computed ages within the expected ranges based on historical atmospheric concentrations and the previously discussed results for both tracer datasets. A simplified model for piston flow was used to compute ages, despite assumptions in the model that cannot adequately account for flow-path mixing that exists in the bedrock. The apparent ages computed by TracerLPM along with the environmental tracers used in the model optimization are shown in table 12. The groundwater ages range from the early 1940s for well WK–439 to the late 1990s for wells WK–426 and WK-431.

Groundwater age provides insight into recharge dynamics across the study area, with implications for groundwater sustainability. Higher amounts of younger water indicate areas with more frequent recharge events and shorter flow paths to supply wells, whereas older water reflects areas with slower recharge rates and longer flow paths. Understanding this age distribution within the aquifer helps identify areas more prone to depletion or at higher risk of rapid transport of contamination from the surface. This information supports the development of management strategies aimed at protecting long-term water availability and quality.

Parameter Sensitivity

As part of the PEST calibration process, sensitivity values were determined for all the individual calibrated parameter values (n = 308) within the six examined hydraulic parameter groups. The parameter groups included in the sensitivity analysis were horizontal hydraulic conductivity, horizontal anisotropy, vertical anisotropy, specific yield, specific storage, and drain conductance. Sensitivity values for the 308 model parameters were calculated by PEST to identify those that had the greatest influence on model results and model fit and those that were insensitive or had little to no effect on the model. Parameters with higher sensitivity values have more impact on how the model matches water-level and flow observations, whereas parameters with lower sensitivity values may reflect that there are insufficient relevant data to constrain parameter estimates or simply that the corresponding property plays a limited role in the flow system.

Both model layers had sensitive pilot points for K_h in the crystalline bedrock and overlying regolith (figs. 10 and 11). The most sensitive parameters were within the bedrock (layer 2), indicating variability in fractured-bedrock hydraulic conductivity is an important parameter to simulate the flow system dynamics. Horizontal and vertical anisotropy parameters for the bedrock layer were also among the most sensitive calibration parameters in the model (table 14). The sensitivity of the anisotropy parameters reflects the directional flow of groundwater in the fractured bedrock, where flow is highly dependent on the orientation and connectivity of fractures, as well as the vertical connection with the regolith. The high sensitivity of horizontal anisotropy in the regolith layer highlights the importance of lateral flow toward streams in the model. This flow is influenced by the orientation of relict structural geologic features and clay horizons within the regolith, as well as topography-induced hydraulic gradients. For simulating streamflow observations, the storage parameters for the regolith/transition-zone layer (layer 1) had the highest composite sensitivities, a reflection of the regolith’s capacity to store and release water as the primary reservoir in the groundwater system--critical in maintaining groundwater balance and flow. The very low storage capacity in the bedrock layer is reflected in the relative insensitivity of the bedrock storage parameters in the model such that iterative changes to the parameter value had little to no effect on the simulation results. Drain conductance for the four streamgage drainage basins were among the most insensitive parameters, likely because of hydraulic head gradients primarily driving the flow of groundwater discharge.

The sensitivity analysis identifies the most influential parameters in this groundwater-flow model simulation. This analysis also helps to identify areas where new data collection could be beneficial for improving model accuracy and reducing uncertainty. For example, including additional regolith monitoring wells with associated aquifer testing data will provide more data points to enhance the resolution of the aquifer storage properties, which were found to be highly sensitive in layer 1. Similarly, focusing on local areas of interest and those with more known variation in hydraulic conductivity data could improve the model’s simulation of groundwater flow. Given that horizontal anisotropy and vertical anisotropy were the most sensitive parameters, new data collection could be prioritized where existing information for these properties is lacking, especially in the bedrock layer.

Residuals of Water Levels and Flow Rates

Calibration of the groundwater-flow model was executed across 82 stress periods for initial steady-state conditions and the transient period 2000–19 by using field observations that included 35,749 historical groundwater levels from 1,100 wells (28 regolith or transition-zone wells and 1,072 bedrock wells; fig. 12) and 560 quarterly stream base flow measurements from four streamgage sites in the model area (fig. 2). Most of these groundwater levels consisted of single measurements at domestic wells, with few at monitoring-network wells with higher frequency measurements. The steady-state stress period contained 6,853 water levels. The subsequent transient stress periods had a range of 18–1,833 water levels with a mean of 357 per stress period.

The 14,802 water-level observations from the 28 regolith wells ranged about 420 ft, from a low of 135.32 ft to a high of 555.56 ft above NAVD 88 (table 15; fig. 13). The absolute value of the mean residual (difference between observed and simulated water levels) for the regolith wells was 8.1 ft with a root mean square error (RMSE) of 9.4 ft. The 20,947 observed water levels from the 1,072 bedrock wells ranged about 713 ft, from a low of 31.38 ft to a high of 744.48 ft above NAVD 88 (table 15; fig. 14). The absolute value of the mean residual for the bedrock wells was 19.9 ft with an RMSE of 26 ft. The normalized standard deviation (the weighted residuals divided by the range of water levels) was 0.01 and 0.02 for the regolith and bedrock layers, respectively, indicating a good calibration fit that accounted for the variation in water-level data. A normalized standard deviation less than or equal to 0.1 reflects the ideal that most of the residuals fall within 10 percent of the observational range (Kuniansky and others, 2004; Fine and Kuniansky, 2014).

A target calibration criterion of ±20 ft was used for the observed water levels considering hydrogeologic context and data uncertainty (Hill and Tiedeman, 2007). This criterion balances the need to reasonably represent regional trends in water levels while acknowledging the complexities within a fractured crystalline bedrock groundwater system at this scale. The ±20 ft threshold reflects about 5 percent of the observed water-level range within the regolith layer and about 3 percent of the range within the bedrock layer. For the regolith groundwater levels (14,802 total), 3 percent were outside of the ±20 ft calibration target (table 15), with 288 measurements exceeding the −20 ft tolerance and 75 exceeding the +20 ft tolerance. For the bedrock groundwater levels (20,947 total), 38 percent were outside of the ±20 ft calibration target, with 4,945 measurements that exceeded the −20 ft tolerance and 2,957 measurements that exceeded the +20 ft tolerance. The higher residuals in the bedrock layer reflect the inherent challenges of modeling fractured-bedrock groundwater systems, where spatially variable hydraulic conductivity and fracture density complicate calibration efforts.

Several additional factors may also contribute to the observed residuals involving the data uncertainty of the well locations and altitudes and the hydraulic conditions at the time of measurement. Hydraulic condition data remarks were not available for most of the domestic well water-level measurements, which may have been obtained during active pumping, after being recently pumped, or influenced by nearby pumping conditions and not representative of ambient conditions. The horizontal coordinates for many of the historical well locations were established using topographic maps having 10-ft contour intervals; thus, there is some inherent uncertainty in the reported water-levels to land surface at those mapped locations. The water level below land surface measurements used for model calibration were updated to water-level altitudes above NAVD 88 by using a 10-meter (32.81 ft) DEM obtained from The National Map (U.S. Geological Survey, 2021) that was upscaled to the resolution of the model grid to minimize uncertainty.

The simulated water-level surfaces of model layers 1 and 2 for the end of the 2019 calibration period are shown in figures 15 and 16, respectively. Groundwater flows perpendicular to equal altitude contour lines in the direction of the steepest downward slope. The water-level surface reflects the gaining stream conditions that prevail within the study area. Contours that form V-shaped patterns around streams convey how groundwater interacts with streamflow, where contour patterns that point upstream suggest groundwater discharge into streams (gaining stream conditions) and those pointing downstream would indicate potential groundwater recharge (losing stream conditions). Groundwater flow paths can be somewhat inferred from contour line spacing, with denser contours indicating local flow systems driven by the topographic highs within small basins.

The water-level surfaces simulated for the regolith and bedrock layers have similar contour patterns (figs. 15 and 16), reflecting the hydraulic connection between the two layers. Some slight differences exist in upland areas where water levels in the bedrock are below water levels in the regolith, reflecting the downward vertical hydraulic gradient in recharge areas. Large surface water features such as B. Everett Jordan Lake, Falls Lake, the Neuse River, and the Cape Fear River show some influence on the fractured-rock groundwater system (figs. 1, 15, and 16). The western boundary of the Triassic sedimentary basin shows a steep drop in water-level altitude from the higher topography underlain by crystalline bedrock in the northwest part of the model area into the sedimentary basin as groundwater flows toward B. Everett Jordan Lake and Falls Lake (figs. 2, 15, and 16). Areas near public supply wells show drawdown--resulting in cones of depression--ranging from 20 to more than 100 ft below regolith water levels and extending out to 1,500 ft from the supply well within the model simulation. The largest drawdowns were simulated in southeastern Wake County.

Two wells within the regolith model layer, WK–332 and OR-691 (USGS station 354315078300101 and 355944079013401; U.S. Geological Survey, 2022; table 2 and fig. 15) were selected for comparison between simulated and observed hydrographs (fig. 17). Well WK–332, screened within the regolith, is located in eastern Wake County near the Neuse River. Observed groundwater levels for well WK–332 were available from 2005 to 2009 and beginning in 2019 when continuous monitoring resumed, with a water level range of 169–174 ft above NAVD 88 (fig. 17A). Simulated groundwater levels at this location generally match the observed levels and ranged from 177 to 183 ft above NAVD 88. This agreement indicates that the model adequately represents hydraulic conditions in the regolith near this major discharge feature.

Well OR–691, screened within the transition zone between the regolith and bedrock and grouped with the regolith model layer, is in eastern Orange County (fig. 15). Observed groundwater levels for well OR–691 from 2013 to 2019 ranged from 443 to 451 ft above NAVD 88 (fig. 17B). Simulated groundwater levels at this location also show reasonable agreement, with a range of 429–450 ft above NAVD 88, which indicates that simulated hydraulic gradients were adequate where the transition zone contributes to the groundwater-flow system.

For the bedrock model layer, two wells (WK–438 and WK–368; table 2; fig. 16) were also selected to compare simulated and observed hydrographs (fig. 18). Well WK–438, located near the southern Wake County border, is open within the fractured bedrock with observed groundwater levels that ranged from 276 to 293 ft above NAVD 88 over the period 1982–2019. Simulated groundwater levels for well WK–438 are within 20 ft of the observed variability, ranging from 301 to 312 ft above NAVD 88, which would indicate somewhat reasonable model performance for long-term water level trends in this area of the model (fig. 18A).

Well WK–368, also open within the fractured bedrock, is located in northern Wake County in an area affected by supply well withdrawals (Chapman and others, 2011). Observed groundwater levels for well WK–368 for 2008, 2011–12, and 2019 show large variability, ranging from 187 to 315 ft above NAVD 88 (fig. 18B). However, the simulated groundwater levels at this location, with a range of 356 to 366 ft above NAVD 88, are much higher than the observed levels and do not reflect the observed drawdown effects from nearby pumping wells. This discrepancy highlights the model’s limitation in simulating local drawdown impacts away from the immediate location of the supply well. Refined grid size with higher resolution in the hydraulic conductivity field of the bedrock may improve the simulation of these localized pumping impacts.

The spatial distribution of groundwater-level residuals across the model domain for model layer 1 (regolith and coastal plain sediments) and layer 2 (bedrock) are shown in figures 19 and 20, respectively. Positive residuals indicate that simulated values are lower than observed values, and negative residuals indicate that simulated values are higher than observed values. Positive and negative residuals were generally well distributed across the model domain with no clear spatial pattern. Smaller residuals within the acceptable calibration criterion (−20 to 20 ft) tended to occur in areas near surface water features, such as the Neuse River, B. Everett Jordan Lake, and Falls Lake (figs. 1, 19, and 20). This pattern suggests better model performance in discharge areas where flow dynamics are more stable compared to upland areas with more variable recharge. Areas with limited observation data from single well measurements tend to have larger residuals--underscoring the importance of having a more robust temporal and spatial water-level dataset to further improve model calibration.

Simulated discharge to streams in the model was compared with stream base flow estimates. The monthly base flow estimates determined from the hydrograph separation analysis of four USGS streamgage sites with multiple years of daily streamflow records performed by Antolino and Gurley (2022) form the basis of the groundwater-flow calibration. A calibration criterion of ±50 percent of the observed flow in million cubic feet per day (Mft³/d) was applied to assess the fit between the observed and simulated flows. Three of the four drainage basins have, on average, acceptable fits for the observed and simulated data. Site 02087359, with the smallest of the four drainage basins (fig. 2), had observed base flows during the simulation period ranging from 0.19 to 2.2 Mft³/d (fig. 21A). Simulated base flows were within calibration criteria with a slightly smaller range from 0.8 to 1.9 Mft³/d. Streamgage site 02097314 had observed base flows that ranged from 0.9 to 13.7 Mft³/d, with simulated base flows that ranged from 1.5 to 5.2 Mft³/d and were within the calibration criteria (fig. 21B). Streamgage site 02087324 had observed base flows ranging from 1.4 to 12.1 Mft³/d and simulated flows that ranged from 3.2 to 8.1 Mft³/d (fig. 21C). Streamgage site 02088500, with the largest of the four drainage basins, had simulated flows outside the criterion that ranged from 6.9 to 22.9 Mft³/d yet were within the range of observed flows for the simulation period from 0.08 to 32.7 Mft³/d (fig. 21D). Hydrographs of simulated and observed flows for all four drainage basins show that simulated flows largely only capture the central tendency of the variability in observed base flow trends and low flows within each basin.

Groundwater Budget

The calibrated groundwater-flow model was used to simulate a regional groundwater budget for Wake County for 2000–19. The simulated budget highlights the major components of inflow and outflow and the resulting changes in aquifer storage for both the regolith and bedrock model layers. These components are recharge, groundwater discharge to surface water bodies (base flow), supply well withdrawals, and transfers between the regolith and bedrock layers, as well as with the regional groundwater system outside the county. The dynamic balance of these flow exchanges varies in response and magnitude to changing seasonal patterns across the simulation period, providing insights into the groundwater system under variable hydrologic conditions.

Flow rates, in million gallons per day, for the individual groundwater budget components simulated for all transient stress periods are summarized in table 16. The average budget component flows into and out of the regolith and bedrock layers are schematically illustrated in figure 22. For the overall budget, the mean total inflows equal the mean total outflows for the regolith (536.8 Mgal/d) and bedrock (181.9 Mgal/d) layers. On average, recharge (303.7 Mgal/d) composes the dominant inflow, whereas base flow (296.7 Mgal/d) represents the dominant outflow within the groundwater system. This balance reflects near-steady-state conditions for the groundwater system under average climate conditions.

Aquifer storage is substantially higher in the regolith layer than in the bedrock layer. The mean flow within the groundwater system is 64.6 Mgal/d into regolith storage and 1.3 Mgal/d into bedrock storage (table 16). The influence of seasonal extremes on the groundwater system is illustrated in figure 22 for a representative wet period (winter 2018–19) and a representative dry period (summer 2007). The above-average recharge of 753.2 Mgal/d under wet conditions for winter 2018–19 led to higher base flow to streams (430.1 Mgal/d) and an increase in regolith storage (315.9 Mgal/d). In contrast, during the dry conditions of summer 2007, recharge rates declined to 20.0 Mgal/d, which resulted in decreased base flow and large losses from aquifer storage. During wet periods, excess recharge replenishes regolith storage and, to a much lesser extent, bedrock storage. During dry periods, regolith storage acts as a buffer, supplying groundwater to sustain streamflow and other outflows when recharge is low.

The temporal variability of the simulated water budget components for each seasonal stress period over the 20-year transient simulation period is illustrated in figure 23. Periods of seasonal extremes, shown by peaks and lows in recharge inflow, further illustrate the dynamics of the groundwater budget. One of the wettest periods simulated occurred from summer 2018 to winter 2018–19, characterized by the highest seasonal recharge rates in the simulation period. This period followed high antecedent recharge conditions that increased aquifer storage and supported higher base flow. The driest period simulated, from spring 2007 to fall 2007, had low rates of recharge, causing lower baseflow and proportionally larger contributions from aquifer storage. These periods highlight the sensitivity of the groundwater system to changes in recharge and the critical role of regolith storage in buffering against climatic variability.

Recharge exhibited substantial variability, ranging from 5.4 Mgal/d in spring 2002 (driest period) to 753.2 Mgal/d in winter 2018–19 (wettest period), with a mean of 303.7 Mgal/d (table 16; fig. 23). This variability reflects the direct response of infiltration through the unsaturated zone from precipitation events-- which often occur as short, intense storm events. In contrast, base flow showed less variability than recharge, ranging from 222.7 Mgal/d during the driest period to 430.1 Mgal/d during the wettest period, with a mean discharge rate of 296.7 Mgal/d. The lower variability for the base flow budget component reflects the storage capacity of the groundwater system, where the amount of time groundwater takes to move through the aquifer along various flow paths buffers the flashy extremes of recharge inflow.

The largest changes in storage within the groundwater system ranged from a maximum loss of 246.7 Mgal/d from total aquifer storage into the flow system to a maximum gain of 335 Mgal/d from the flow system into aquifer storage (table 16). This range highlights the sensitivity to variations in recharge and discharge conditions. The mean net change in storage (“to storage” minus “from storage”) was a gain of 1.4 Mgal/d, whereas the median indicated a loss of 26.6 Mgal/d. This discrepancy reflects the impact of more variable recharge events on the mean compared to the relatively constant discharge to streams that is better captured by the median. Storage for the fractured-rock groundwater system is primarily concentrated in the overlying regolith, with minimal storage capacity in the fractured bedrock. The regolith accumulates storage during wet periods yet is prone to notable losses during dry periods. This dynamic is crucial in understanding the groundwater system response to prolonged droughts, particularly in areas dependent on groundwater withdrawals.

Groundwater withdrawal to supply wells was one of the smallest components in the county-wide simulated budget. Reported supply withdrawals ranged from 0.8 to 11.4 Mgal/d with a mean of 6.6 Mgal/d over the simulation period 2000–19 (table 16). Whereas these withdrawals appear minimal at the county scale, their localized impact at a small basin scale can be substantial where supply wells are concentrated. The magnitude of observed drawdowns in monitoring-network wells near large community-supply wells, such as USGS monitoring well WK–368 (figs. 5 and 7), highlights the limited storage capacity of the fractured-bedrock groundwater system in such areas.

Forecast Scenario Simulations

The calibrated groundwater-flow model for Wake County was used to simulate future scenarios that included 204 stress periods representing seasonal data from winter 2019–20 through fall 2070, hereinafter referred to as the “forecast period.” These scenarios incorporated projected development patterns based on the FUTURES model and climate variability under the RCP 4.5 and RCP 8.5 emissions scenarios for three GCMs that represent a range of dry to wet climatic conditions (Meentemeyer and others, 2013; Sanchez and others, 2020b). The simulations were to assess potential future impacts of climate variability and increased development on recharge and base flow, the two major components of the Wake County groundwater budget. The model forecast results were compared against simulated recharge and base flow for 2000–19 in the calibrated model (referred to as the “calibrated simulation”) to highlight potential future changes in recharge and base flow projections from 2000–19 conditions. The discussion focuses on the results of the RCP 4.5 and RCP 8.5 emissions scenarios, which reflect the aggregated results of the three GCMs used to represent the overall range in potential future climate conditions. Aggregating forecast models by RCP scenario may reduce apparent spread of variability within individual model runs, yet this presentation emphasizes the central tendency and overlapping range across all runs for each RCP scenario. The difference in the number of simulations and timeframes should be considered when comparing variability and central tendency between the calibrated and aggregated forecast simulation results. Percentiles of the simulation results, including the medians (50th percentiles), interquartile range (IQR, difference between the 25th and 75th percentiles), and the extremes (10th and 90th percentiles), are displayed as box plots in figures 24 and 25 to examine temporal and seasonal changes in the distributions of the recharge and base flow forecast results. The percentile dataset is provided in table 1.1 in appendix 1.

The temporal results for simulated recharge and base flow rates across the calibrated and forecast periods are shown in figure 24. Data were arbitrarily grouped by 5-year intervals, with a 6-year interval used for 2065–70, to aid visualization and comparisons across the entire forecast period. Examination of recharge results indicate a slight downward shift in recharge rates between the IQRs of the calibrated simulation dataset and the forecast simulation dataset (fig. 24A). This shift likely is influenced by a combination of factors, including changes in climate conditions projected by the GCMs, such as increased evapotranspiration because of warmer temperatures, and impacts of increased urban development, such as a larger impervious footprint reducing infiltration and altering runoff dynamics. Additionally, differences in the spatial resolution of the climate datasets used in the SWB model (Antolino, 2022) that provided recharge input to the groundwater model simulations may contribute to the shift. The finer 1-km resolution of the Daymet dataset, used for the calibration period, captures localized variability that would include high recharge events that contribute a notable portion of the total recharge component. In contrast, the coarser 4-km resolution of the MACAv2 dataset at 4-km scale, used for the forecast simulations, likely smooths out these events, leading to a slightly more muted recharge estimate. A similar downward shift was noted in the base flow results between the calibrated and forecast simulation datasets (fig. 24B). Because base flow is coupled to seasonal and annual recharge dynamics, changes to values in future recharge may propagate as changes in groundwater discharge to streams.

No apparent trends were noted in the forecast simulation period of 2020–70 that would suggest potential long-term changes in future recharge or base flow (figs. 24A and 24B, respectively). However, future recharge and base flow rates could have considerable yearly variability as indicated by the temporal variability within and among individual time intervals for RCP 4.5 and RCP 8.5 forecast simulations. End members in forecast recharge rates across all forecast time intervals, based on the 10th and 90th percentiles, ranged from 20.8 Mgal/d (2045–49) to 658.1 Mgal/d (2050–54) for RCP 4.5 and 23.0 Mgal/d (2020–24) to 619.5 Mgal/d (2035–39) for RCP 8.5. These ranges were similar to the low of 21.5 Mgal/d (2000–04) and high of 643.0 Mgal/d (2015–19) for the calibrated simulation. For all time intervals, median recharge values for RCP 4.5 varied by a factor of about 1.7, ranging from 161.4 to 270.1 Mgal/d. For RCP 8.5, median recharge values varied by a factor of about 1.3, ranging from 215.3 to 278.3 Mgal/d. Median recharge values for the calibrated simulation varied by a factor of about 1.6, ranging from 240.3 to 374.6 Mgal/d.

Similarly, end members in base flow rates across all forecast time intervals ranged from 183.1 Mgal/d (2045–49) to 374.3 Mgal/d (2030–34) for RCP 4.5 and 193.0 Mgal/d (2040–44) to 356.3 Mgal/d (2055–59) for RCP 8.5 (fig. 24B). The calibrated simulation exhibited a wider range in base flow end members, from a low of 226.7 Mgal/d (2010–14) to a high of 408.8 Mgal/d (2015–19). Across all time intervals, median base flow values for RCP 4.5 varied by a factor of about 1.2, ranging from 246.2 to 286.4 Mgal/d. For RCP 8.5, median base flow values varied by a factor of about 1.1, with a range of 244.0 to 274.7 Mgal/d. Median base flow values for the calibrated simulation varied by a factor of about 1.2, ranging from 269.2 to 318.5 Mgal/d. Overall variability between RCPs 4.5 and 8.5 was similar across the forecast period for both recharge and base flow. Simulated recharge rates showed more variability than did base flow rates, likely reflecting the episodic nature of precipitation that directly influences groundwater recharge. Meanwhile, groundwater discharge was moderated by aquifer storage, which buffered short-term variability and provided more consistent base flow to streams. Although the model effectively captured the central tendency of base flow conditions, it underrepresented observed low base flows (fig. 21A–D)—a limitation to be considered when interpreting forecast results of base flow, especially for periods of extended dry conditions with low recharge when extreme low flows are likely.

Seasonal analysis of recharge and base flow indicated a similar pattern of higher rates occurring in winter and fall and lower rates in spring and summer for both the calibrated and forecast simulation periods (fig. 25). Recharge rates (fig. 25A) and base flow rates (fig. 25B) were most variable during the winter season, reflecting the interaction of increased precipitation and reduced evapotranspiration. A downward shift in recharge and base flow for RCPs 4.5 and 8.5 is noted for all seasons except for recharge in winter and spring. As stated above, a downward shift between the calibrated and forecast datasets likely reflects a combination of dataset resolution differences and potential climate-driven changes. Median recharge values in winter for RCPs 4.5 and 8.5 (427.5 and 410.7 Mgal/d, respectively) were slightly higher than the winter median of 391.0 Mgal/d for the calibrated simulation that may suggest a potential for increased winter recharge under future climate scenarios. Across seasons, median recharge values varied by a factor of about 3.1 for RCP 4.5, a factor of about 3.7 for RCP 8.5, and a factor of about 3.3 for the calibrated simulation. Median base flow values for RCP 4.5, RCP 8.5, and the calibrated simulation all varied across seasons by a factor of 1.2. These seasonal patterns suggest that, while overall recharge rates do not deviate outside historical observational ranges, potential increases in winter recharge may slightly offset potential reductions during other seasons.

Forecast scenario simulations highlight the sensitivity of the regional groundwater system to potential future climate variability, with some impacts on rates of recharge and base flow. These impacts are driven by shifts in precipitation patterns, evapotranspiration rates, and overall seasonal dynamics. Winter recharge may increase under future conditions, partially offsetting reduced recharge during spring and summer. During periods of extended dry conditions when recharge becomes more intermittent and driven by specific episodic events, aquifer storage in forecast simulations gradually depletes because of limited recharge. This depletion lowers groundwater levels and consequently reduces groundwater discharge to streams. Whereas the buffering capacity of aquifer storage helps stabilize base flow rates during short-term dry conditions, extended dry periods can diminish this capacity over the duration of drought. Monitoring groundwater levels, changes in aquifer storage, recharge rates, and streamflow (especially during low flow conditions) provides valuable data to aid with tracking changes, validating and updating projections, and informing adaptive management strategies for sustaining groundwater resources under variable climate conditions.