Continental-scale analysis of shallow and deep groundwater contributions to streams

Groundwater discharge zones establish hydrologic connectivity between aquifers and streams, controlling water quantity and quality, especially during low-flow periods. Such discharge shapes stream thermal regimes and supports groundwater-dependent ecosystems that are vulnerable to climate change and contamination. Groundwater can buffer stream temperatures, creating thermal refugia, but recent observations show widespread, spatially heterogeneous stream warming, partly due to variable groundwater contributions to flow. The depth of contributing groundwater is critical: shallow groundwater (within ~6 m of land surface) exhibits strong seasonal thermal variability and is more sensitive to land use and seasonal drawdown, whereas deeper groundwater shows minimal annual thermal variation, distinct chemistry, and more seasonally stable discharge. Efficient, scalable methods to infer groundwater source-depth are limited. Many existing approaches quantify baseflow or temperature sensitivity but do not distinguish depth, while chemical/isotopic methods are resource-intensive. In the absence of groundwater influence, annual stream temperatures closely track air temperatures; departures in amplitude and timing (phase) indicate influences from shallow or deep groundwater or dam operations. Shallow groundwater contributions impart phase lags and variable damping; deep groundwater contributions dampen the annual signal with minimal phase shift. This study applies a refined, process-based signal analysis to classify 1729 U.S. stream sites by groundwater signature (shallow, deep, atmospheric, or major dam) using multi-year paired air–stream temperature records. Objectives are to: (1) compare temperature-signal-based categorization to baseflow indices, (2) map spatial patterns and landscape drivers of groundwater discharge characteristics, and (3) assess long-term stream temperature trends (14–30 years) among signature categories, providing continental-scale inference to inform predictions of stream responses to climate and land-use change.

Prior work shows that groundwater substantially influences stream thermal regimes and habitats, providing thermal buffering and refugia for cold-water species. Observed stream warming is widespread but spatially heterogeneous, influenced by groundwater contributions. Groundwater depth governs thermal variability: shallow aquifers transmit pronounced annual temperature signals and are sensitive to land use, contamination, evapotranspiration, and drought, whereas deeper groundwater is thermally stable on annual scales yet responds to long-term climate trends and often has distinct chemistry that can affect surface water quality. Existing broad-scale methods (e.g., hydrograph separation, linear air–water temperature sensitivity) quantify baseflow or temperature sensitivity but typically cannot resolve source-depth. Chemical and isotopic end-member approaches can infer mixing but are resource intensive and may not isolate shallow flow paths without additional hydrologic characterization. Signal processing of paired air–water annual temperature has emerged as a scalable approach to infer hydrogeologic controls on stream thermal regimes and to differentiate shallow versus deep groundwater influences based on amplitude damping and phase lags.

Study design: Classified 1,729 stream/river sites across the conterminous U.S. into atmospheric, shallow groundwater, deep groundwater, or major dam signatures using multi-year paired air and stream temperature annual signals. Sites span 1st–9th order streams across 21 physiographic provinces. Major dam-influenced sites were identified and excluded from groundwater signature analyses. Data sources and screening: From ~4,000 candidate stream temperature stations, 1,811 met criteria: within 25 km of a NOAA air temperature station and ≥2 consecutive post-2010 years of data with no ≥30-day gaps; after QA/QC, 1,729 sites remained. Stream temperature data came from USGS NWIS, NorWeST, and SHEDS repositories. Daily air temperatures were obtained from NOAA GHCN-Daily via the rnoaa package, preferentially using the nearest station meeting 75% annual completeness and 75% overlap with stream records; a second nearest station was used when needed (n=191). Watershed attributes were linked using NHDPlus COMIDs and EPA StreamCat; GAGES II provided distance to major dams, watershed slope, and Hydrologic Disturbance Index (HDI). Signal processing: For each site, average daily air and stream temperatures were fit with a static sinusoid (a·sin(t)+β·cos(t)+C) using SciPy curve fitting to minimize RMSE. From fitted coefficients, annual amplitude (A) and phase (φ) were computed. Derived metrics: amplitude ratio Ar = A_stream / A_air and phase lag Δφ = φ_stream − φ_air, converted to days (positive values indicate stream lagging air). Daily stream temperatures <1°C were excluded to avoid freeze–thaw decoupling; stream values >60°C were removed. Negative phase lags between 0 and −10 days were set to 0; lags < −10 days (n=25) were removed due to likely management artifacts or poor fits. Signature thresholds: Conservative, theory- and field-informed thresholds were used: deep groundwater signature if Ar < 0.65 with minimal phase lag (no explicit lag threshold); shallow groundwater signature if phase lag Δφ ≥ 10 days (indicative of shallow groundwater mixing) with variable Ar; atmospheric signature if 0.65 ≤ Ar ≤ 1.1 and Δφ ≈ 0; records with Ar > 1.1 were removed as probable pairing or measurement errors. Sites with Δφ > 40 days or within 25 km downstream of major dams were categorized as major dam signatures and excluded from groundwater signature analyses. Very low Ar values (<0.4) were manually checked for proximity to major dams upstream (within ~30 km). Baseflow comparison: For the subset of 554 NWIS sites with concurrent discharge data, baseflow index (BFI; percent of daily baseflow averaged over the temperature record) was computed using the USGS-R DVstats bfi function to compare hydrograph-separation-derived baseflow with temperature-based signatures. Trend analysis: Long-term trends were assessed for 184 sites with ≥14 complete years within 1990–2019 (after removing anomalous/managed series). Monthly mean temperatures were analyzed for annual and summer (June–August) periods using nonparametric Theil–Sen slopes with deseasoning (R openair TheilSen), robust to outliers and accounting for seasonality. This 14–30 year window covers at least two ENSO cycles; potential PDO influences were noted but not found to differentiate western from other sites in this dataset. Human disturbance metrics: Impervious cover (StreamCat) and HDI (GAGES II; composite of dams, reservoir storage change, canals, road density, proximity to pollutant discharges, freshwater withdrawals, undeveloped land fragmentation) were compared across signature categories.

• Scope and classification: Of 1,729 sites, 305 were classified as major dam signatures and excluded from groundwater signature analyses. Among the remaining 1,424 sites, 39% (n=556) exhibited pronounced groundwater signatures: 47% deep groundwater (n=264) and 53% shallow groundwater (n=292). Atmospheric signatures comprised 868 sites.
• Signal metrics by category:
  • Deep groundwater: mean amplitude ratio Ar = 0.54 (σ = 0.10); mean phase lag = 3.8 days (σ = 3.4).
  • Shallow groundwater: mean Ar = 0.59 (σ = 0.18); mean phase lag = 16.6 days (σ = 6.6).
  • Atmospheric: mean Ar = 0.85 (σ = 0.12); mean phase lag = 2.3 days (σ = 2.7), not different from zero.
• Baseflow comparison (n=554 with discharge): Atmospheric sites had significantly lower median BFI (0.69) than shallow groundwater (0.79) and deep groundwater (0.86) sites, consistent with groundwater driving baseflow contributions.
• Spatial heterogeneity: Groundwater signatures occur across stream sizes, physiographic provinces, and within sub-watersheds. Within the North Fork Clearwater–Lake Creek watershed (ID–MT), headwaters tended toward shallow signatures while mainstem valleys showed deep signatures; the watershed outlet shifted to an atmospheric signature. All four examined sites in upper Lake Creek had shallow signatures (>15-day phase lags), aligning with observed temperature impairments potentially linked to warming shallow groundwater.
• Human disturbance patterns: Atmospheric signature sites were associated with higher impervious cover and higher HDI (median HDI: atmospheric 16; deep groundwater 9; shallow groundwater 5.5), suggesting reduced stream–groundwater connectivity in more developed, low-slope catchments.
• Long-term temperature trends (14–30 years; n=184):
  • Atmospheric (n=132): >50% warming (70 sites); warming rates 0.01–0.09 °C yr−1 (mean 0.04 °C yr−1).
  • Shallow groundwater (n=29): 45% warming; rates 0.01–0.10 °C yr−1 (mean 0.04 °C yr−1).
  • Deep groundwater (n=23): 52% stable; 26% warming (6 sites) at 0.01–0.05 °C yr−1 (mean 0.01 °C yr−1); 22% cooling.
  • Summer trends: 70% of shallow groundwater sites warmed; 61% atmospheric; 43% deep groundwater.
• Interpretation: Shallow groundwater-dominated streams exhibit more frequent and stronger warming trends and lower baseflow fractions than deep groundwater-dominated streams, indicating greater vulnerability to thermal stress and low-flow conditions. Deep groundwater contributions impart greater thermal buffering and sometimes cooling trends, emphasizing resistance to short- to medium-term climate warming.

The study demonstrates that the depth of groundwater sources exerts a first-order control on stream thermal regimes at continental scales. Deep groundwater contributions dampen annual stream temperature signals with minimal phase shift and are associated with higher baseflow indices and greater thermal stability, including a larger fraction of stable or even cooling trends over recent decades. In contrast, shallow groundwater contributions impart significant phase lags, lower baseflow fractions, and a higher prevalence of warming trends, particularly in summer when ecological thermal stress is most acute. These results empirically validate theoretical and modeling predictions regarding the differential buffering capacity of shallow versus deep groundwater and highlight the sensitivity of shallow systems to seasonal heatwaves, drought, and withdrawals. Human landscape modification correlates with reduced groundwater signature prevalence: sites draining catchments with greater impervious cover and higher HDI more often show atmospheric signatures, consistent with diminished infiltration, reduced recharge, and altered shallow groundwater dynamics due to pumping and stormwater infrastructure. Heterogeneity within basins underscores the need for fine-scale characterization; even adjacent reaches can differ in dominant source-depths because of geomorphology and hydrologic connectivity. Collectively, the findings inform process-based management, suggesting that assumptions equating high baseflow with thermal stability can be misleading unless groundwater source-depth is explicitly considered.

This work introduces a scalable, process-based framework using paired air–stream annual temperature signals to classify the dominant groundwater source-depth influencing streams across the conterminous United States. Nearly 40% of non-dam sites exhibit pronounced groundwater signatures, split between shallow and deep sources. Deep groundwater-dominated streams show stronger thermal buffering and higher baseflow fractions, while shallow groundwater-dominated streams are more prone to warming and reduced baseflow, increasing vulnerability to thermal stress and low flows under climate change and water use. The approach leverages widely available datasets, enabling managers and modelers to incorporate groundwater source-depth into predictions of stream temperature, flow, and water quality. Future directions include: refining regional/geomorphic thresholds for amplitude ratio and phase lag; integrating discharge variability and groundwater level data to disentangle climatic versus withdrawal effects; expanding high-resolution mapping within priority basins for cold-water species; and coupling with chemical/isotopic tracers to quantify mixing fractions and residence times. Incorporating source-depth explicitly into hydrologic and ecological models will improve forecasting and guide targeted mitigation, such as protecting deep-groundwater-fed refugia and managing shallow aquifer withdrawals and recharge.

• Temporal coverage: Long-term trend analysis spans 14–30 years, sufficient for multiple ENSO cycles but short relative to PDO and other decadal variability, potentially affecting trend detectability.
• Signal model simplification: A static sinusoidal fit extracts the dominant annual component but may not capture complex nonstationary or sub-annual dynamics; fitting imprecision can introduce small non-physical phase lags that were zeroed or filtered.
• Threshold generality: Conservative amplitude ratio (0.65) and phase-lag (10 days) thresholds are based on select well-studied systems; optimal thresholds may vary by region, hydrogeomorphology, and channel conditions (e.g., shading, confinement).
• Data pairing and quality: Air–water station pairing within 25 km and removal of extreme Ar (>1.1) mitigate mismatches but some pairing or measurement errors may persist; sub-daily variability is not used (daily means only). Negative phase lags beyond −10 days were excluded as likely managed or poor-fit cases.
• Dam influence handling: Sites with Δφ > 40 days or within 25 km downstream of major dams were categorized as dam-influenced and excluded from groundwater signature analyses; residual management effects upstream or beyond set distances may remain.
• Baseflow comparison subset: BFI analyses were limited to 554 sites with concurrent discharge records and use an empirical separation method, which carries its own uncertainties and may not directly quantify groundwater source-depth.
• Confounding hydrologic factors: Discharge variability, droughts, and groundwater withdrawals can modulate stream temperatures and baseflow, particularly in shallow systems, complicating attribution solely to source-depth.