Study Area

The Fountain Creek alluvial aquifer, which includes the area previously designated as the Widefield aquifer (Edelmann and Cain, 1985), is composed of unconsolidated and permeable Quaternary-age sand and gravel alluvium deposited along the main valley of Fountain Creek. The alluvial aquifer is located in an urban area of Colorado Springs and surrounding communities in a semiarid climate with a mean annual precipitation of approximately 18.48 inches per year based on gridded datasets (PRISM Climate Group, 2020). The alluvium is generally coarse grained and lacks substantial intervening clay confining layers (AECOM, 2017), with a mean hydraulic conductivity (K) of about 830 feet per day (ft/d; Edelmann and Cain, 1985). The Fountain Creek alluvial aquifer is underlain by the relatively impermeable bedrock of the Late Cretaceous Pierre Shale, which prevents the downward movement of groundwater and promotes horizontal flow in the permeable alluvial deposits (Edelmann and Cain, 1985). The top of the Pierre Shale beneath the Fountain Creek alluvial aquifer is a Quaternary erosional surface exhibiting paleochannels and ridges (Radell and others, 1994; Lewis, 1995; AECOM, 2017). Four named tributary paleochannels drain from the Peterson Space Force Base to the main stem of the Fountain Creek alluvial aquifer (fig. 1): Sand Creek, Broadway, Schlage, and Windmill Gulch (AECOM, 2017). The Sand Creek, upper Broadway, and Windmill Gulch paleochannels coincide with existing stream channels, whereas the lower Broadway and Schlage paleochannels do not (AECOM, 2017). Groundwater flow in the center of the aquifer generally follows the top of the impermeable bedrock surface, and variably saturated groundwater conditions occur in several tributary alluvium paleochannels (Lewis, 1995; AECOM, 2017). The alluvium thickness ranges from 0 to about 100 feet, and the saturated thickness ranges from 0 to about 45 feet (Radell and others, 1994). Saturated thickness within each paleochannel is generally greatest along their central axis (AECOM, 2017). In the main valley of Fountain Creek, saturated thickness in the alluvial aquifer is greatest along a northwest to southeast trending axis between the tributaries of Sand Creek and Windmill Gulch (Lewis, 1995). Alluvial, colluvial, and eolian deposits located on ridges, outside of the alluvial aquifer extent on the top of bedrock surface were often unsaturated (Lewis, 1995; AECOM, 2017).

Study area streamflow varies seasonally because of snowmelt, contributions from groundwater base flow, and discharges from wastewater treatment plants (WWTP). Peak streamflow in Fountain Creek generally occurs May through July (USGS, 2021) coincident with snowmelt runoff from the adjacent mountains. Low streamflow periods in Fountain Creek were characterized by inputs from WWTP (Mau and others, 2007) and by groundwater discharge. In addition to Fountain Creek, Canal 4 is an artificial surface water body in the study area. Canal 4 runs along the northeastern extent of the alluvial aquifer and indicated by previous studies to be a source of groundwater recharge (fig. 1; Newman and others, 2021).

Previous investigations indicated variable spatial representations of the saturated extent of the Fountain Creek alluvial aquifer. For instance, in early investigations by Edelmann and Cain (1985) and Cain and Edelmann (1986), the saturated extent was indicated to closely follow the modern stream channel, with minor saturation in paleochannels. Topper and others (2003) mapped alluvial aquifers on a regional scale and in the vicinity of the study area indicated saturated conditions could be expected farther up the paleochannels than the original extents by Edelmann and Cain (1985) and Cain and Edelmann (1986). More recent investigations specific to the area near Peterson Space Force Base (AECOM, 2017) further delineated saturated conditions in both the paleochannels and in the main stem. The estimate of saturated extent from AECOM (2017), however, does not include paleochannels on the western side of the aquifer extending toward the mountain front mapped by Topper and others (2003). In this investigation, the datasets of Topper and others (2003) and AECOM (2017) were combined to create a new estimate of the extent of the alluvial aquifer, which is illustrated in figure 1 and included in Russell and Newman (2024).

Groundwater and surface-water interaction is a principal aspect of the water budget for the Fountain Creek alluvial aquifer because of the potential hydrologic connection to Fountain Creek. Previous water budgets for the Fountain Creek alluvial aquifer indicated the primary source of inflow to the aquifer is recharge from Fountain Creek streamflow losses that occur mainly in the reach between the tributaries of Spring Creek and Sand Creek (fig. 1; Edelmann and Cain, 1985). Other sources of recharge included infiltration from precipitation, irrigation return flows, and irrigation canal seepage; underflow from tributary alluvium and upgradient parts of the Fountain Creek alluvial aquifer; and infiltration from septic systems and sewage lagoons (Edelmann and Cain, 1985). Edelmann and Cain (1985) focused on the occurrence and distribution of nitrogen in groundwater, because at that time, a large percentage of the surface-water flow in Fountain Creek was derived from the Colorado Springs WWTP, and wastewater contained elevated nitrogen concentrations (Edelmann and Cain, 1985). Groundwater and surface-water interaction, as well as the distribution of nitrate in the Fountain Creek alluvial aquifer, were the focuses of the Lewis (1995) groundwater-quality study in the Fountain Creek valley.

Groundwater discharge components from the Fountain Creek alluvial aquifer used in this study include evapotranspiration (ET), discharge to Fountain Creek, discharge from groundwater wells, and groundwater flow to downgradient parts of the aquifer. Discharge to ET is facilitated by a shallow groundwater table and dense stands of cottonwood trees in the Fountain Creek valley. Groundwater and surface-water interaction is hypothesized to be controlled by depth to the groundwater table and the top of bedrock, which results in areas of groundwater recharge and discharge along Fountain Creek. Between the tributaries of Spring and Sand Creeks, Fountain Creek was previously considered a losing stream that recharged the aquifer based on water-budget calculations (Edelmann and Cain, 1985). Downstream from Sand Creek, groundwater generally discharged to the stream, and Fountain Creek was considered a gaining stream (Edelmann and Cain, 1985). Groundwater discharge from wells occurs primarily in the main valley of the Fountain Creek alluvial aquifer. Pumping wells provided water to the communities of Security, Widefield, and Fountain, and domestic, stock, and commercial wells also pump groundwater from the Fountain Creek alluvial aquifer (U.S. Army Corps of Engineers, 2018a; U.S. Army Corps of Engineers, 2018b; U.S. Army Corps of Engineers, 2018c).

Previous Studies

Several previous studies and ongoing data collection by the USGS in the Fountain Creek valley contributed to understanding the groundwater and surface-water conditions in the study area. Groundwater hydrology and water quality were investigated by Edelmann and Cain (1985), Cain and Edelmann (1986), and Lewis (1995) to support an understanding of water budgets and groundwater-flow directions. In the 1980s, the USGS developed and has since maintained a computer program estimating streamflow gains and losses in Fountain Creek, known as the transit-loss model (Kuhn, 1988, 1991; Kuhn and others, 1998, 2007; Kuhn and Arnold, 2006; Colarullo and Miller, 2019). The transit-loss model is used by the Colorado Division of Water Resources and local water providers for daily native water accounting, transmountain diversions, transit losses, and water allocations along Fountain Creek and Monument Creek (not shown on any maps). Results from the transit-loss model provide a detailed accounting of stream gains and losses, which can be used to quantify gaining and losing reaches of Fountain Creek.

In 1998, the USGS, in cooperation with Colorado Springs City Engineering, began a study of the Fountain Creek and Monument Creek Basins to characterize the effects of wastewater-treatment effluent and storm runoff on water quality and suspended-sediment conditions (Mau and others, 2007). Water-quality samples were collected at 11 main-stem sites between 1981 and 2001 and at 14 tributary sites in 2003. Suspended-sediment samples were collected daily at 7 main-stem sites from 1998 through 2001, daily at 6 main-stem sites from 2003 through 2006, and intermittently at 13 tributary sites from 2003 through 2006 (Mau and others, 2007). Water-quality samples were analyzed for nutrients, fecal coliform, and trace elements, and the results were evaluated in terms of base flow, normal flow, and stormflow. Stormflow samples exhibited concentrations of fecal coliform that frequently exceeded water-quality standards 1998–2006 on main-stem and tributary sites, and suspended sediment, nitrogen, phosphorous, and trace-element load substantially increased during stormflow, indicating these constituents were transported to surface water primarily by storm runoff (Mau and others, 2007).

As part of the study discussed in this report, an extensive water-quality dataset was collected from groundwater wells and surface waters September 2018 through May 2019 to evaluate sources of recharge to the aquifer as well as sources of solutes. This dataset included major and trace elements, rare earth elements, pharmaceutical and wastewater indicator compounds, and groundwater age tracers. The water-quality dataset was published in the National Water Information System (NWIS) database (USGS, 2021) and interpretations were reported in Newman and others (2021). In general, Newman and others (2021) found surface waters were sources of aquifer recharge based on the presence of pharmaceutical and wastewater indicator compounds in groundwater wells in close proximity to both Fountain Creek and Canal 4 (fig. 1). Spatial distribution of groundwater age estimates also indicates Fountain Creek and Canal 4 were sources of recharge, because groundwater near these surface water bodies was generally younger than groundwater in the center of the aquifer. Groundwater ages estimated in Newman and others (2021) ranged from several months to approximately 21 years old. The findings of Newman and others (2021) assist in conceptualizing the hydrologic conditions discussed in this report.

Purpose and Scope

The purpose of this report is to characterize the groundwater flow and movement in the Fountain Creek alluvial aquifer 2018 through 2020 during decreased groundwater withdrawals from pumping wells. Additionally, this report provides hydrologic conditions for the alluvial aquifer pertinent to the potential for transport of solutes. An integrated approach was used to evaluate groundwater hydrologic processes, assess groundwater and surface-water interactions, estimate hydraulic aquifer properties using aquifer testing, and complete numerical groundwater-flow simulations. Groundwater-flow directions were determined from continuous and periodically measured groundwater-level elevations, and groundwater and surface-water interactions were evaluated using synoptic streamflow measurements. Groundwater-flow processes were evaluated quantitatively using statistical water-level modeling approaches and a numerical groundwater-flow model; the latter is used to simulate recharge and discharge from the alluvial aquifer, interactions with surface water, and groundwater-flow paths. Input and model files facilitating quantitative analyses using statistical approaches and the numerical groundwater-flow model were compiled into a USGS data release (Kisfalusi and Newman, 2022) and a model archive associated with this report (Russell and Newman, 2024). This data release and model archive may be used to reproduce this report results or be updated for future conditions. Synoptic water-quality data were collected and interpreted to indicate physical and chemical processes occurring in the aquifer, and water-quality sampling analyses and results are reported in Newman and others (2021) and available in the NWIS database (USGS, 2021).

Groundwater Hydrology

An understanding of the groundwater-flow directions in the Fountain Creek alluvial aquifer was determined from groundwater-level elevation measurements. As part of the study, groundwater levels were measured at 31 locations (fig. 2; table 1), many were measured in previous USGS studies (for example, Radell and others, 1994; Lewis, 1995). Discrete depth-to-groundwater measurements were collected at each location and then converted to groundwater-level elevations in feet (ft) above the North American Vertical Datum of 1988 (NAVD 88) using the depth below land surface of the groundwater level and the land-surface elevation of the location in NAVD 88. Groundwater levels were measured using either an electric or a steel tape according to standard methods described in Cunningham and Schalk (2011).

Three locations were additionally equipped with pressure transducers to collect continuous groundwater-level data. Pressure transducers were installed at wells MW1-1, T11-MW001, and TH-22 (fig. 2). Pressure transducers were set to record groundwater levels at one-hour intervals. During each site visit, the continuous data were downloaded and discrete depth to groundwater was measured to allow drift correction of continuous data recorded by the pressure transducer according to methods described in Cunningham and Schalk (2011). Both discrete and continuous groundwater-level elevation data were used to assess the distribution and movement of groundwater within the study area.

Statistical groundwater-level modeling of groundwater-level elevations was completed to determine the primary physical processes controlling groundwater flow. Specifically, a transfer-function-noise (TFN) approach was useful for indicating the effect of external hydrologic stressors (for example, recharge, groundwater pumping, and streamflow) on groundwater-level elevation fluctuations in individual wells (Bakker and Schaars, 2019). The TFN approach uses statistical signal analysis algorithms to decompose the observed fluctuation in groundwater-level elevations into response functions, which represent different hydraulic stressors. For example, the effect of a groundwater pumping well on groundwater-level elevation fluctuations in an observation well may be evaluated using a Hantush function (Bakker and Schaars, 2019). Using the TFN approach, the proportion of the total observed variation in groundwater-level elevation fluctuations may be attributed to different hydrologic stressors.

Groundwater-level modeling was conducted with the TFN approach for the Fountain Creek alluvial aquifer using the Pastas package (Collenteur and others, 2019) implemented in the Python programming language. The Pastas package was used because of the number of response functions included and the ability to programmatically conduct TFN modeling. Input requirements were time series of groundwater-level observations, estimated recharge, and pumping discharge. The groundwater-level model was completed using discrete data from groundwater-level elevation observation locations. Groundwater-level data that predated the study data for this report were also included in the analysis. Two groundwater-level models were completed: one in which only recharge was considered and one in which only groundwater pumping was considered. These two groundwater-level models allow for the two primary hypothesized stressors in the study area to be evaluated individually. The groundwater-level model considering only recharge was constructed by estimating a daily time series of precipitation and ET for the study area using the Climate Engine online tool (Huntington and others, 2017). A time series of point recharge at each observation well was calculated as precipitation minus ET. The groundwater-level model using only groundwater pumping was constructed by aggregating pumping rates for nearby observation wells on a daily basis. Groundwater-level model code, inputs, and outputs are provided in a USGS data release and a model archive associated with this study (Russell and Newman, 2024).

Groundwater and Surface-Water Interactions

Evaluation of the interaction between groundwater and Fountain Creek was included in this study because previous studies indicated the exchange of groundwater and surface water were major terms in the local water budget (Edelmann and Cain, 1985; Cain and Edelmann, 1986). When groundwater-level elevations near the streambed are higher than the stream-surface elevation, groundwater may flow into the stream, known as a gaining stream. Conversely, if the groundwater-level elevation is lower than the stream-surface elevation, then water may flow from the stream into the groundwater, known as a losing stream. Gaining streams are groundwater discharge locations, whereas losing streams are groundwater recharge locations (Winter and others, 1998).

To assess streamflow gain or loss in the study area, synoptic streamflow measurements were made in March 2019 at sites corresponding to the Fountain Creek transit-loss model where established USGS streamgages do not exist (fig. 3 of this report; Kuhn and Arnold, 2006). Streamflow was measured using an acoustic Doppler velocimeter according to methods described in Rehmel (2007) and Turnipseed and Sauer (2010). To supplement discrete streamflow data, the daily mean streamflow from three continuous USGS streamgages (Fountain Creek at Colorado Springs, CO, 07105500; Fountain Creek below Janitell Road below Colorado Springs, CO, 07105530; and Fountain Creek at Security, CO, 07105800) was extracted from the National Water Information System (NWIS) database (USGS, 2021) on the date of the synoptic streamflow measurements. Data supplemented from the NWIS database provided a more spatially detailed understanding of groundwater and surface-water interactions.

Using streamflow measured from the transit-loss model at upstream and downstream locations, net streamflow gain or loss from each reach was calculated according to equation 1 (modified from Simonds and Sinclair, 2002)

Values of ΔQ less than zero indicate the stream reach is losing (a source of groundwater recharge), and values of ΔQ greater than zero indicate the stream reach is gaining (a point of groundwater discharge).

The net streamflow gain or loss along each reach is subject to errors in streamflow measurement, which are accumulated by different sources of error in each measurement (Rosenberry and LaBaugh, 2008). Total error for each reach gain loss was calculated according to the error propagation formula (Harmel and others, 2006):

Error in individual measurements was quantified using the interpolated variance estimator described by Cohn and others (2013) and recorded by the acoustic Doppler velocimeter. Individual measurement and net errors were calculated in the unit of the streamflow analysis, in ft³/s, by multiplying the quantitative error (5 percent) by the measured streamflow, yielding an estimate in ft³/s of the absolute error of the measurement.

A number of potential diversions and inflows to Fountain Creek in the study area, primarily from irrigation ditches divert surface water from Fountain Creek, and WWTP discharge treated effluent into Fountain Creek (Kuhn and others, 2007; Mau and others, 2007). Streamflow, at each diversion and inflow, was not measured as part of this study, but values for each were extracted from the transit-loss model and the NWIS database for the date of the synoptic measurements (March 3, 2019). These values were used for computing streamflow gain or loss assuming an error in each diversion and inflow value of 10 percent, which is the likely amount of error, assuming the diversion and inflow values had qualitative indicators of poor given by the field data collector (Turnipseed and Sauer, 2010).

Aquifer Testing

Hydraulic properties of an aquifer govern the flow of groundwater through the aquifer from the recharge area to the discharge area. Hydraulic properties such as K, specific storage (S_s), and specific yield (S_y) may be estimated from reference materials for a given aquifer type (Freeze and Cherry, 1979), computed with inverse methods from spatially distributed groundwater-level observations (Shapoori and others, 2015), or directly estimated through aquifer testing (Kruseman and de Ridder, 1990; Butler, 1997). Aquifer testing consists of stressing the hydrologic system, typically by either adding or removing water from a well and observing the time-dependent system recovery to the pretest conditions. Data from the test may then be compared to various analytical solutions selected based on the hydrogeologic conceptualization of the aquifer and used to calculate hydraulic properties.

Single-well aquifer tests were conducted in 11 observation wells during July 2019 (fig. 4; table 2). Three wells were located in the Broadway paleochannel, two wells were located in the Windmill Gulch paleochannel, and six wells were located within the main stem of the Fountain Creek alluvial aquifer. Six wells were tested with traditional slug tests, and five wells had a combination of traditional slug tests along with single-well pumping tests (table 2). Each traditional slug test consisted of a falling-head (slug-in) and a rising-head (slug-out) test. In these tests, different initial displacements were used to evaluate the presence of complicating factors such as unidirectional flow or low-K well skins (materials within the well bore that inhibit connection between the well and the aquifer; Butler, 1997). Each solid slug was dropped below the groundwater table for the falling-head test, and when the groundwater level reached the equilibrated initial groundwater level, the slug was pulled above the groundwater table for the rising-head test (Butler, 1997; Cunningham and Schalk, 2011). Single-well pumping and recovery tests were conducted using a Grundfos Redi-Flo 2 pump. Pumping was initiated once static groundwater levels were confirmed by pressure transducer measurements. When groundwater-level displacement was deemed substantial based on the length of pumping period and observed groundwater-level change rate, pumping was discontinued, and recovery was observed until approximate static conditions (within 0.06 ft of initial groundwater level) were reached.

For each test type, a vented pressure transducer was placed near the bottom of the well (below the solid slug or the pump) to measure groundwater-level displacement (Cunningham and Schalk, 2011). Pressure transducers were placed as far as possible away from the groundwater-level perturbation to minimize signal noise during data collection. Groundwater levels were measured at a frequency of 0.5 seconds with observations made before, during, and after an aquifer test to differentiate drawdown and recovery from groundwater-level changes caused by environmental factors. Well completion diagrams for those observation wells included in aquifer testing are included as part of a USGS data release (Kisfalusi and Newman, 2022).

The aquifer testing goal was to quantify hydraulic properties in the Fountain Creek alluvial aquifer and to evaluate spatial variability in the hydraulic properties. The primary hydraulic property estimated using the aquifer-test results was K, although in some instances estimation of storage properties (for example, S_s and S_y) may also be possible. Aquifer-testing data were analyzed using the software AQTESOLV version 4.5 (Duffield, 2007). Analysis of slug-test data used the methods of Bouwer and Rice (1976), Springer and Gelhar (1991), and model of Hyder and others (1994). The Bouwer and Rice (1976) method is best suited for application to unconfined aquifers and has various assumptions relating to aquifer penetration, well screen position, and well flow behavior (Butler, 1997). The Springer and Gelhar (1991) method was used to account for inertial effects in the wells and the oscillatory responses seen in the study area. The Springer and Gelhar (1991) solution also incorporates frictional well loss in small-diameter wells (Butler, 1997). The Hyder and others (1994) model was useful when the presence of a variable behavior between rising-head and falling-head tests was observed, potentially because of interferences between the aquifer material and the well bore (Butler, 1997), commonly known as “well-skin effects,” which were evident during one test. The single-well pumping tests were analyzed using the Cooper and Jacob (1946) method to evaluate both the pumping and recovery data. Halford and others (2006) indicate the Cooper and Jacob (1946) method was the analytical solution most applicable to single-well pumping tests for unconfined conditions, regardless of partial penetration or other potentially complicating factors.

Groundwater-Flow Simulations

Given the complex nature of the Fountain Creek alluvial aquifer and multiple potential hydrologic stressors in the system (for example, groundwater pumping and groundwater and surface-water interactions) it was necessary to apply a groundwater-flow model to simulate and predict groundwater movement. In this setting, a groundwater-flow model may be used to evaluate conceptual groundwater recharge models and discharge sources, quantify changing groundwater budgets under transient conditions, and quantitatively predict groundwater-flow paths. A numerical groundwater-flow model was constructed for the Fountain Creek alluvial aquifer using the finite-difference code MODFLOW with the Newton formulation solver (MODFLOW-NWT; Niswonger and others, 2011). The model was spatially discretized into individual cells whose size and properties were specified through user input, with user-specified time steps for temporal discretization of the MODFLOW model. Discretization of the groundwater-flow model and application of boundary conditions are described in detail in the following paragraphs and illustrated in figure 5. The numerical groundwater-flow model and all associated input files are included as a USGS data release associated with this report (Russell and Newman, 2024).

The numerical model of the Fountain Creek alluvial aquifer has 1 layer of 291 rows and 254 columns of 200-ft by 200-ft size cells, with a total of 17,610 active cells (fig. 5). A one-layer model was used because the alluvium overlies Upper Cretaceous Pierre Shale with low permeability, which serves as a vertical no-flow boundary for the groundwater-flow system. The model grid was projected and aligned relative to the Colorado State Plane Coordinate System (North American Datum of 1983, units: feet, zone: Colorado Central Zone).

The active model grid extent was based on the extent of the Quaternary alluvium deposits that compose the Fountain Creek alluvial aquifer. The active extent of the model grid did not include the paleochannels present to the northeast of the primary aquifer. The paleochannels were excluded because of transient unsaturated conditions and thin saturated thickness, as indicated by groundwater-level monitoring, though these paleochannels may be dominant sources of chemical constituents to the aquifer (Newman and others, 2021). To evaluate groundwater flow from the paleochannels, the General-Head Boundary (GHB) package (Harbaugh and others, 2000) was used to represent flow originating in the Broadway and Windmill Gulch paleochannels.

The model simulates time periods as model stress periods, and during each model stress period, hydrologic stressors were specified and held constant for the stress period duration. The Fountain Creek alluvial aquifer model contains 241 stress periods. Of these 241 stress periods, the initial stress period represents a steady state period, and the other 240 stress periods were transient, simulating each month from 2000 to 2019. For the initial stress period of the numerical model, hydrologic stressors and groundwater-flow rates were assumed to be constant, and the stress period represents mean annual data throughout the transient stress periods. The outputs from the steady-state model stress period were used as inputs for the subsequent transient stress period.

The Fountain Creek alluvial aquifer system numerical groundwater-flow model is composed of multiple hydrologic boundaries representing areas of inflow or outflow. These hydrologic boundaries were simulated through MODFLOW packages (fig. 5). The hydrologic boundaries simulated in this numerical model were separated into two different types: head-dependent and specified-flux boundaries. Head-dependent boundaries allow flow into or out of the model based on differences between user-specified groundwater levels. Specified-flux boundaries allow flow into or out of the model based on user-specified flux rates. Recharge and well withdrawals were simulated in the model using specified-flux boundaries, and evapotranspiration, streams, reservoirs, and lateral groundwater flow were simulated using head-dependent boundaries.

Recharge for the numerical model was simulated using the Recharge (RCH) package (Harbaugh and others, 2000) and precipitation data from the Parameter-elevation Regressions on Independent Slopes Model (PRISM) dataset from Oregon State University (PRISM Climate Group, 2020). The PRISM data were modified using the ArcGIS Resample tool (Esri, 2021), which interpolated the precipitation data to fit the model grid cell size. Distributed recharge is caused by infiltration of precipitation (rain and snow) into the aquifer system. The amount of water infiltrating the land surface is affected by several factors, such as precipitation rates, evapotranspiration rates, unsaturated-zone permeability and moisture capacity, and land surface slope; therefore, only a fraction of the amount of precipitation from the PRISM dataset recharges the aquifer. This fraction is applied to the recharge arrays as a multiplier for each simulated stress period.

To determine the recharge amount applied to the aquifer, the RORA program (Rutledge, 2000) was applied to USGS streamgage streamflow data from the NWIS database (USGS, 2021) throughout and nearby the active model area. RORA is a program that estimates recharge for the drainage basin of a specified streamgage using a recession-curve displacement method. This method uses the upward displacement of the streamflow-recession curve during groundwater recharge periods, such as storms, to determine the rate of groundwater recharge to an idealized, homogenous aquifer (Rutledge, 2000). Estimation of recharge using the RORA program differs from the point estimates of recharge used as inputs to the TFN model because RORA estimates spatially integrated recharge rates for a drainage basin above a streamgage, whereas point estimates of recharge at wells used in the TFN model are not spatially integrated but vary through time. The RORA program was used for 11 USGS streamgages from 2000 to 2019 (fig. 5). The recharge estimates output by RORA ranged from 2.14 to 3.47 inches per year (in/yr). These outputs were then weighted based on the basin size and period of record of the streamflow data, with the mean weighted estimate of recharge being 2.53 in/yr, or about 13.72 percent of the mean annual precipitation (18.48 in/yr) from the PRISM dataset. The recharge multipliers within the recharge package were constrained so recharge for the numerical model would not deviate too far from the mean weighted estimate of recharge, with a final mean annual recharge rate of 2.68 in/yr used for the numerical model transient stress periods.

Evapotranspiration is the process by which water is directly transferred into the atmosphere through evaporation and through plant transpiration; evapotranspiration represents a water-budget outflow from the alluvial aquifer. The Evapotranspiration (EVT) package (Harbaugh and others, 2000) was used to simulate this process in the numerical model. Evapotranspiration within the model was limited to cells with cottonwood trees identified based on satellite imagery from Google Earth (Russell and Newman, 2024) and found mostly along the main reach of Fountain Creek (fig. 5). A simulated root-zone depth of 0.40 ft was used, and the simulated ET rate varied each month based on seasonal ET rates from remotely sensed imagery showing cottonwood trees in the lower Colorado River Basin by Jetton (2008).

Given the evidence of groundwater and surface-water interactions in the aquifer as indicated by previous studies (Edelmann and Cain, 1985; Cain and Edelmann, 1986; Lewis, 1995), it is hypothesized streamflow losses may be a substantial groundwater recharge source. A detailed water-quality analysis of groundwater and surface water, carried out in conjunction with this study, indicated wells near Fountain Creek displayed water-quality signatures consistent with streamflow recharge (Newman and others, 2021). To evaluate exchanges between groundwater and surface water, streams within the active model area were simulated using the Streamflow-Routing (SFR) (Niswonger and Prudic, 2005), the Drain (DRN) and the River (RIV) (Harbaugh and others, 2000) packages. Concrete-lined streams within the active model area, as shown on satellite imagery from Google Earth (Russell and Newman, 2024), were not simulated, as the concrete acts as an impermeable boundary between the stream and the underlying aquifer. An example of a concrete-lined stream not simulated as interacting with the aquifer is Windmill Gulch between Canal 4 and Fountain Creek (fig. 1). The main reach of Fountain Creek and its tributaries were simulated using the SFR package. The SFR package uses Darcy’s Law, streambed conductance, groundwater-table elevation, and stream stage to determine flow between the stream and the surrounding aquifer. Specified inflows within the SFR package allow for base flows from tributaries outside the model active area to be incorporated into the model, allowing for more accurate simulated streamflows. SFR calculates the simulated base flows along specified segments and routes the base flows downstream until they were outside of the active model area or captured by a different hydrologic boundary.

All SFR stream segments had a thickness of 1.0 ft and a Manning’s roughness coefficient of 0.025, based on common assumptions for the SFR package (Niswonger and Prudic, 2005). Streambed conductivity for each stream segment ranged from 0.01 to 50.0 ft/d. This wide range of conductivity is because of the Pierre Shale outcropping in some parts of the stream, whereas the rest of the streambed is predominantly fine to coarse sand. Simulated stream widths ranged from 20.0 to 95.0 ft and were estimated using satellite imagery from Google Earth (Russell and Newman, 2024). Streambed elevations were modified using USGS 5-meter digital elevation model values to allow for accurate simulated base flows (USGS, 2020). Monthly base flow inflows from outside the active model area included Fountain Creek, Shooks Run, Spring Creek, Sand Creek, B Ditch Drain, and Clover Ditch and were input using the SFR package. Data gaps in the monthly base flow data were filled with the mean base flow for the corresponding month and location. The main reach of Fountain Creek has several points of artificial inflows and diversions (Kuhn and Arnold, 2006). However, these points were not simulated in the numerical model, because they were determined to be minor components of the water budget.

The SFR package routes flow downstream to connected stream reaches. Because of this functionality, streams can only be simulated within the package if the upstream reaches were higher than or equal in elevation to the downstream reaches. This elevation requirement can cause an issue with artificial canals, such as Canal 4. Therefore, the canals and other streams not included in the SFR package were simulated using the DRN and RIV packages (fig. 5). The difference between the two packages is the RIV package allows for flow into and out of the groundwater systems to a simulated surface-water feature, whereas the DRN package only allows for outflow from the groundwater system to the simulated surface-water feature. Flow between the aquifer and surface-water features in the RIV package is dependent on stage of the surface-water feature, hydraulic conductivity of the feature-aquifer interconnection, and the head at the node in the cell underlying the surface-water feature. Flow in the DRN package is dependent on the feature’s conductance, elevation, and head in the cell. Flow from the aquifer into the drain feature will only be simulated if the head in the drain cell is above the drain elevation. Because the DRN package only allows for flow out of the aquifer and into surface-water features, the RIV package was predominantly used for streams and canals. However, in areas where the numerical model experienced hydraulic heads above land surface, the DRN package was implemented.

All the RIV and DRN package cells were channelized, meaning the cell elevations were set below land surface to simulate a stream channel. For the RIV package, the cell elevations were set to 4 ft below land surface, and the DRN package cell elevations were typically set to 1.0 ft above the bottom of the model (representing the interface between the alluvial aquifer and the Pierre Shale). The DRN cells have the potential to remove water from the alluvial aquifer based on the gradient between the aquifer and DRN in each cell location. However, because the aquifer is typically thin in the vicinity of most DRN cells (AECOM, 2017), the DRN package is limited in its ability to cause unrealistic water budget outflows from the model. The ability for the RIV and DRN packages to transmit water is governed by hydraulic conductivity values. The hydraulic conductivity values for the RIV package varied based on location and ranged from 107 to 17,500 feet per day (ft/d). The hydraulic conductivity values for the DRN package also varied based on location and ranged from 2,978 to 44,800 ft/d. For the RIV package, a stage elevation of 1.0 ft above the streambed was set for all simulated streams.

Lateral inflows into the model from external water sources (Big Johnson Reservoir and the Broadway and Windmill Gulch paleochannels) were simulated using the GHB package (Harbaugh and others, 2000). Flow between the GHB and aquifer is the product of hydraulic conductivity of the boundary and the difference between the head in the aquifer and head in the GHB. Big Johnson Reservoir, an approximately 300-acre artificial lake and one of the major surface-water features in the active model area (fig. 5), was simulated using a GHB head of 5,800 ft and a hydraulic conductivity of 15.0 ft/d. Water-quality data indicated leakage from Big Johnson Reservoir may recharge the aquifer (Newman and others, 2021). Although water levels in Big Johnson Reservoir declined beginning in 2016 when draining the reservoir began for repairs (Steiner, 2017), noble gas data in groundwater from wells downgradient from the reservoir indicated recent equilibrium with the atmosphere, whereas other nearby groundwater wells did not (Newman and others, 2021). Recent atmospheric equilibrium indicates the reservoir is a source of recharge to the aquifer.

Lateral inflow and outflow were other substantial components of the modeled water budget, because the paleochannels were not directly simulated in the active model area and because the alluvial aquifer continues to the south, outside of the active model domain. Groundwater flow from paleochannels to the part of the aquifer simulated in the numerical model was simulated using several GHBs (fig. 5), and a GHB was used to represent the southern end of the model domain. Based on a previous water budget for the aquifer, Edelmann and Cain (1985) found groundwater flow likely exits the study area to the south. Head elevations for the GHBs were set using data from nearby USGS groundwater-level observation wells 384653104451901 and 383952104405601. The monthly means of observed groundwater levels, recorded in feet below land surface for 384653104451901 were used for the GHB heads in the Broadway and Windmill Gulch paleochannels, and the monthly means of observed groundwater levels for 383952104405601 were used for the GHB heads at the southern boundary. Hydraulic conductivity of 15.0 ft/d was used, and model inputs were compiled into the GHB package (Russell and Newman, 2024).

A primary hydrologic stressor in the Fountain Creek alluvial aquifer is groundwater pumping for municipal use. Several water-supply districts in the area operated pumping wells, prior to the PFAS compounds discovery (U.S. Army Corps of Engineers, 2018a; U.S. Army Corps of Engineers, 2018b; U.S. Army Corps of Engineers, 2018c; McDonough and others, 2021). The distribution of 38 known pumping wells in the study area is illustrated in figure 6A, along with reported groundwater pumpage through time (fig. 6B). Groundwater pumping data are included in the data release associated with this study (Russell and Newman, 2024). Most pumping wells were located in the central part of the alluvial aquifer, near Fountain Creek. Total groundwater pumpage, for 2000 to 2018 (fig. 6B) varied, but trends in specific water-supply districts are apparent, with groundwater pumpage from the Security district decreasing beginning in 2013, coinciding with an increase in pumping from the Widefield district. The effect of groundwater pumping on specific wells was assessed via the TFN groundwater-level model (as discussed in the “Groundwater Hydrology” section of this report), and the variation in groundwater pumping was used to inform the transient calibration of the groundwater-flow model, as described in the next paragraph. Groundwater pumping wells were simulated in the model using the Well (WEL) package (Harbaugh and others, 2000). The mean value of all well pumping during the transient stress period was used in the steady-state stress period.

Following initial parameterization, the numerical model was calibrated to achieve better agreement between groundwater-level observations and base flows. This model calibration process includes manual and automated calibration steps. Manual calibration was done by changing input parameter values until a user specified and subjective reasonable fit for groundwater-level elevations was achieved, after which automated calibration using computer programs was used to complete the calibration process based on quantitative criteria. The computer programs used in the automated calibration step continuously change user-specified parameters to minimize the differences between simulated and observed groundwater-level elevations and base flows. These user-specified parameters were set within a predetermined range that conceptually match the properties of the Fountain Creek alluvial aquifer. To complete this parameter estimation, the program PEST++ iterative ensemble smoother (PEST++IES) was used, which applies an iterative ensemble smoother algorithm to minimize a user-defined objective function that describes discrepancies between the model simulations and observations (White and others, 2020).

The automated calibration process for the alluvial aquifer used groups of input parameters and groups of observations. Input parameters were split into 8 groups of 477 total parameters, and observations were split into 5 groups of 2,100 total observations. Of the 477 parameters, there were 241 parameters for recharge, 93 parameters for horizontal hydraulic conductivity (K_h), 93 parameters for the vertical to horizontal hydraulic conductivity ratio (K_v/h), 21 parameters for SFR streambed conductance, 14 parameters for ET, 4 parameters for SFR stream widths, 4 parameters for RIV riverbed conductance, 3 parameters for GHB conductance, and 1 parameter each for S_s, S_y, SFR stream thickness, and DRN conductance. The recharge parameters were multipliers applied to each recharge array for every stress period within the numerical model, and the range of the recharge multipliers were limited to only allow realistic rates of recharge to the simulated groundwater system. The 93 parameters for K_h were evenly spaced points throughout the active model domain using a pilot points approach (Doherty, 2003), and the range of K_h was based on aquifer testing data collected as part of this study. The 93 parameters for K_v/h were fixed, meaning they were not changed during the automated calibration process. Most ET parameters were multipliers (13 of 14) used similarly to recharge array multipliers. Twelve of the multipliers represented monthly ET rates, and the other multiplier represented the ET rate for the steady-state stress period of the numerical model. The last evapotranspiration parameter adjusted the root-zone depth in the EVT package. All boundary conductance parameters were limited to realistic boundary conductance values calculated using K values and area of model cells.

Of the 2,100 total observations used in the automated calibration process, 1,136 groundwater-level elevation were observations from 42 wells and 964 were base flow observations from 4 USGS streamgages (Fountain Creek below Janitell Road below Colorado Springs, CO, 07105530; Sand Creek above Mouth at Colorado Springs, CO, 07105600; Fountain Creek at Security, CO, 07105800; and Clover Ditch Drain near Widefield, CO, 07105820; fig. 5). All groundwater-level elevation data were retrieved from the USGS NWIS database (USGS, 2021). The 964 base flow observations were obtained from the USGS NWIS database and processed using the USGS Groundwater Toolbox (Barlow and others, 2017). Base flow was estimated for the two streamgages located along Fountain Creek (07105530 and 07105800; fig. 5) throughout the entire transient stress period of the numerical model, 2000–19. Base flow was estimated for the streamgage along Sand Creek (07105600) seasonally for the months April–September between 2003 and 2014 and transitioned to estimates for each month throughout the year between 2014 and 2019. Base flow was estimated monthly for the streamgage along Clover Ditch Drain (07105820) beginning in November 2017 and extending through 2019. Any periods of missing base flow estimates were filled with the mean monthly base flow for the site using the available data (USGS, 2021). Most of the streamflow measurements made at these streamgages, during the transient period of the numerical model, were rated as fair by the hydrographer conducting the measurement, meaning streamflow measurements were within 8 percent of true streamflow (Turnipseed and Sauer, 2010).

These observations were categorized into five groups, one group for the groundwater-level elevation observations and a separate group for each of the four streamgages. The contribution of the observation groups to the objective function were adjusted to maintain a balance between the observation groups during the calibration process. The base flow observations accounted for 81 percent of the observation group contribution to the objective function, meaning the model was most sensitive to changes in base flow.

The ability of the model to emulate observed groundwater-level elevations was assessed quantitatively using the root-mean-squared-error (RMSE) and the normalized RMSE (nRMSE) calculated according to

The RMSE and nRMSE were useful model calibration indicators, because they quantitatively indicate the ability of the model to simulate the hydrologic system. The calibration goal is to reduce both the RMSE and nRMSE. An nRMSE less than or equal to 10 percent is generally considered to represent a well-calibrated model (Anderson and others, 2015). All results of model calibration and simulations are described in the “Groundwater-Flow Simulations” section of this report.

The approximate groundwater-flow paths of water sourced in pumping wells was assessed using the calibrated model and the particle tracking software MODPATH version 7 (Pollock, 2016). MODPATH is a software package that computes the approximate trajectory of hypothetical particles within the aquifer based on advective transport. Particle tracking may be useful for understanding source zones of groundwater discharged from pumping wells (Moeck and others, 2020). In the context of the Fountain Creek alluvial aquifer, the groundwater to pumping wells source is pertinent because of the occurrence of PFAS compounds in drinking water in the area (Hu and others, 2016; McDonough and others, 2021). MODPATH simulations were created using the steady-state numerical model, as is common for particle-tracking applications in complex aquifers (Moeck and others, 2020). Particle tracking was conducted in reverse mode, meaning a hypothetical particle was placed in each model cell containing a pumping well. The model then was run in reverse to calculate the trajectory of each particle from the pumping wells to the origin of the particle. This analysis essentially outlines the areas of the aquifer contributing groundwater to pumping wells. The Python scripts used for MODPATH model creation and processing are included in the USGS data release associated with this report (Russell and Newman, 2024).

Model Limitations

Limitations to the numerical groundwater-flow model for the Fountain Creek alluvial aquifer include the nature of the spatial and temporal variability of the model discretization, the necessary simplifications used in representing hydrologic boundaries, and the necessary simplifications of the approach used to evaluate groundwater capture by pumping wells. Each of these limitations is discussed in detail in the following paragraphs.

The process of numerical groundwater-flow model creation results in the discretization both spatially and temporally of the physical system into a framework that allows for model calculations to converge. This means hydraulic parameters of the aquifer must be smoothed in a manner fitting the spatial discretization of the system; the result is generally a simplified distribution of K with respect to reality. Temporally variable input parameters, such as precipitation recharge and groundwater pumping, were also smoothed and averaged in a manner allowing these parameters to be used by the model, but this is a simplification of real-world behavior. Spatial and temporal discretization of the model was done to minimize these effects on the ultimate model outputs, but these simplifications mean the model is not a perfect analog for the hydrologic system.

Hydrologic boundaries, such as surface water boundaries, may be implemented in the model in different ways, and the way they were incorporated controls the ability of boundaries to interact with the remainder of the model domain. Several different MODFLOW packages were used to represent surface water: GHB, SFR, RIV, and DRN. The GHB boundary represented the Big Johnson Reservoir, which is reasonable for representation of lakes and reservoirs that were sources of groundwater recharge (Anderson and others, 2015). The SFR package was used for streams expected to have the most direct connection to groundwater, such as along the main stem of Fountain Creek, because this package allows for two-directional flow between groundwater and surface water. The SFR package, however, can create convergence issues in some portions of the model. The RIV and DRN packages were used to represent areas where groundwater and surface-water interactions were hypothesized to compose a smaller water-budget component of the model when compared to areas simulated with the SFR package, which can represent bidirectional flow into and out of the aquifer. This is a necessary simplification to achieve model convergence, and means groundwater and surface-water exchange was unidirectional in areas of the model simulated with the DRN package. This simplification is not expected to cause inaccurate representations of boundaries in the model, based on spatial analysis of water-quality data presented in Newman and others (2021).

The way advective groundwater-flow paths were simulated also included necessary simplifications. One simplification was the steady-state simulation was used to simulate flow paths, which is a common approach in evaluating sources of water to pumping wells (Moeck and others, 2020). Using transient models in particle tracking may result in complications because of starting time and location of particles, the rate of change of heads in the transient model, and the complexity of the model framework (Rayne and others, 2014). For the Fountain Creek alluvial aquifer, results of the steady-state simulated groundwater-level elevations closely resemble the long-term median groundwater-level elevations (fig. 8), indicating the steady-state model is a reasonable representation of the long-term groundwater system and groundwater-flow paths from the steady-state model were likely reasonable representations of actual flow paths under most conditions. Some seasonal variations in flow paths are likely because of seasonal changes in groundwater-level elevations, but these are expected to cause only minor variation in the groundwater-flow paths as illustrated in figure 21. Another simplification of the advective particle tracking model is only the calibrated model was used in simulating flow paths, as opposed to integrating a full uncertainty analysis of flow path density from multiple uncalibrated models (for example, Moeck and others, 2020). In their analysis of particle tracking, Moeck and others (2020) illustrated uncertainty in hydraulic properties of the aquifer results in uncertainty as to flow paths and the delineation of groundwater source areas to wells. The Fountain Creek alluvial aquifer model incorporated an inverse model analysis using PEST++IES to understand uncertainty in the K of the aquifer, but not every uncalibrated model realization from the uncertainty analysis was carried forward for particle tracking. Qualitatively, sensitivity analysis (fig. 17) indicates base flows and groundwater-level elevation simulations were most sensitive near the Sand Creek paleochannel, but sensitivity is low elsewhere in the model domain. Groundwater-flow paths were therefore unlikely to change substantially in the majority of the model domain as a result of uncertainty in the hydraulic property values.

Because the MODPATH simulations use the steady-state model, they cannot be used to estimate time of travel within the transient groundwater system. Use of the steady-state model was needed because of numerical constraints and issues arising from transient simulations, as described in the previous paragraph. Although time of travel between recharge sources and wells cannot be computed with the current model, these results can be used in conjunction with groundwater age dating described in Newman and others (2021). The combination of these methods generally indicates solutes occurring in wells have likely been recharged within the last 10 years, and these solutes could be derived from the locations of the particle sources on figure 21.