Intermittent streamflow generation in a merokarst headwater catchment

Abstract

Intermittent headwater streams are highly vulnerable to environmental disturbances, but effective management of these water resources requires first understanding the mechanisms that generate streamflow. This study examined mechanisms governing streamflow generation in merokarst terrains, a type of carbonate terrain that covers much of the central United States yet has received relatively little attention in hydrological studies. We used high-frequency sampling of precipitation, stream water, and groundwater during summer 2021 to quantify the contributions to streamflow from different water sources and characterize their short-term dynamics in a 1.2 km² merokarst catchment at the Konza Prairie Biological Station (Kansas, USA). Mixing calculations using stable water isotopes and dissolved ions indicate that streamflow is overwhelmingly contributed by groundwater discharge from thin (1–2 m) limestone aquifers, even during wet periods, when soil water and surface runoff are generally expected to be more important. Relationships between hydraulic heads in the aquifers and their contributions to streamflow differed early in the study period compared to later, after a major storm occurred, suggesting there is a critical threshold of groundwater storage that the bedrock needs to attain before fully connecting to the stream. Furthermore, contributions from each limestone unit varied during the study period in response to differences in their hydrogeological properties and/or their stratigraphic position, which in turn impacted both the length of streamflow and its composition. Taken together, we interpret that the subsurface storage threshold and variation in aquifer properties are major controllers of flow intermittency in merokarst headwater catchments.

---

Environmental significance

Intermittent streamflow accounts for up to 50% of global discharge, but we currently lack a full understanding of the mechanisms which generate this streamflow in many physiographic settings. Characterizing streamflow generation for these settings will allow us to better manage and protect their water resources. In this study, we identify mechanisms governing streamflow generation in an intermittent headwater catchment in a merokarst terrain, an understudied subset of carbonate terrains that underlies much of the west- and south-central USA. Our results demonstrate how thin carbonate units within merokarst bedrock can serve as important hydrologic buffers and that the properties of individual units affect the timing and significance of their contribution to streamflow, and therefore the nature of intermittency in the stream.

1. Introduction

Headwater catchments with intermittent streamflow are significant sources of water, nutrients, and energy for both human populations and natural ecosystems worldwide. According to previous studies, headwater streams are estimated to comprise 66–80% of total global stream length^1,2 and intermittent flow is estimated to account for 30–50% of total river discharge.^3,4 However, the small size and high hillslope connectivity of headwater catchments⁵ make these critical waterways highly sensitive to environmental disturbances and thus they are particularly vulnerable to the effects of climate change and human developments.⁶ Successful management strategies aimed at mitigating these effects cannot be developed and applied unless we first understand the processes that govern streamflow generation,^3,7,8 but research on the mechanisms behind headwater streamflow intermittence is lacking in many climatic and physiographic settings.^7,9,10 In this study we consider the sources and timings of streamflow generation in an intermittent headwater stream within a merokarst terrain in order to further our understanding of the mechanisms behind flow intermittency.

Merokarst terrains are a subset of carbonate terrains in which thin, relatively high permeability carbonate units are interbedded with low permeability mudstones or other siliciclastic units. A high proportion of carbonate terrains are merokarst, including significant parts of central North America.^11,12 However, they remain largely understudied compared to other carbonate terrains and may be particularly sensitive to environmental change.^13,14 The rapid transport of water generally seen in carbonate terrains amplifies the sensitivity of the environment to disturbances,¹⁵ exacerbating the responsiveness of headwater catchments therein.

Our research focuses on a headwater catchment within the Konza Prairie Biological Station (hereafter referred to as Konza Prairie) in northeastern Kansas, USA (Fig. 1). Konza Prairie is a part of the broader Flint Hills physiographic province of eastern Kansas and north-central Oklahoma and is classified as a merokarst terrain.¹⁴ Previous studies conducted at Konza Prairie reveal key pieces of information about merokarst streamflow generation and flow intermittency. Costigan et al.⁹ used long-term precipitation and discharge data (1987–2011) to characterize the overall surface water hydrology of Konza Prairie and observed that headwater discharge patterns were disconnected from precipitation patterns. They concluded that threshold-driven processes were causing this disconnect and ultimately leading to intermittent flow behavior. They offer the gradual closing of soil macropores and the filling of groundwater storage zones as two possible mechanisms that might be giving rise to the threshold dynamic, but their results did not allow them to test these hypotheses.

Fig. 1 Map of Konza Prairie Biological Station and watershed N04D. Locations of wells included in this study are indicated with open red circles. Each of these wells also contains a pressure transducer used to monitor groundwater levels. In the well identifiers, EIS 1 is used for wells screened in the Lower Eiss Limestone and EIS 2 is used for wells screened in the Upper Eiss Limestone. Wells 4-6 MOR, 4-6 EIS 1, and 4-6 EIS 2 are marked by a single circle because they are located within 1 m of each other. Approximate locations where select bedrock units outcrop in the stream channel are identified based on previous research.²⁰ Limestone and shale are abbreviated Ls and Sh, respectively, in the legend.

More recently, Sullivan et al.¹⁶ employed a principal-component analysis to identify potential end members that contribute solutes to streams on Konza Prairie. They concluded that the stream water is a mixture of four different sources: surface runoff, soil water, groundwater flowing through the carbonate aquifers, and groundwater flowing along the carbonate-mudstone contacts. Sullivan et al.¹⁶ did not quantify the relative contributions to streamflow from these endmembers in their study, but data collected by Keen et al.¹⁷ suggest that the groundwater inputs are dominant. They analyzed stable isotope compositions of water from soil cores as well as archived spring and surface water samples collected during 2010–2017. Water isotope compositions of stream samples closely resembled those of groundwater discharging at the spring.

These analyses provide valuable insight into the potential water sources that contribute to streamflow generation, but they do not reveal the short-term dynamics of these water sources that might help give rise to flow intermittency. Additionally, they do not capture separate contributions to streamflow from the different stacked carbonate units within Konza Prairie bedrock stratigraphy. The different carbonate units at Konza Prairie vary chemically, meaning that shifting contributions from these aquifers to the stream can impact stream chemistry.^11,18 Furthermore, each aquifer can have significantly different hydrogeologic properties¹⁹ and outcrop patterns, suggesting that they might impart different flow pathways on water moving through them.

In this study, we seek to build on these previous analyses by investigating the relative importance of different water sources, their variation on short time scales, and how this variation might influence stream flow intermittency. To achieve this goal, we employed high frequency sampling of surface and groundwater at Konza Prairie during the 2021 growing season (April–July) to quantify streamflow source contributions and characterize their short-term dynamics. We collected stream samples on a weekly to hourly basis depending on flow conditions and applied mixing analyses to capture small timescale changes to source contributions. We established the endmembers for such mixing analyses using groundwater sampled from numerous monitoring wells along with long-term precipitation, stream, and groundwater chemistry records that are available in the extensive Konza Prairie data catalog. Lastly, we utilized sensors measuring high frequency discharge, precipitation, and groundwater elevation data to provide context for any changes to source contributions and aid in interpretations of threshold-driven behavior that might be controlling flow intermittency.

2. Study area

This study was conducted in watershed N04D within Konza Prairie, which is located in northeastern Kansas (Fig. 1). N04D is a 1.2 km² headwater catchment drained by an intermittent fourth order stream, with land cover largely dominated by native prairie grasses. Woody vegetation including burr oak (Quercus macrocarpa), hackberry (Cercis canadensis), and American elm (Ulmus americana) within the riparian zone of the principal stream while other woody species like rough-leaf dogwood (Cornus drummondii) and smooth sumac (Rhus glabra) have significantly increased in abundance along catchment hillslopes in the decades since Konza Prairie was established.^21,22 N04D is burned every four years as part of ongoing research on the connection between grassland burn frequency and increased woody plant abundance (a process referred to as encroachment), and a burn occurred three days before the first sample collection of this study on April 12, 2021.

Konza Prairie receives an average of 835 mm of precipitation annually.²³ Most of the precipitation falls during the growing season which runs from April–September²⁴ (ESI Fig. SI1A†). Discharge is measured at a gauging station at the outlet of N04D and generally peaks during May–June, coinciding with the peak of precipitation. However, streamflow typically dries out by July, despite the persistence of frequent precipitation over the remainder of the growing season²⁵ (ESI Fig. SI1B†). This drying is likely due to the watershed achieving maximum rates of evapotranspiration, which is estimated to account for 75% of the total annual water budget in Konza Prairie.²⁶ Despite this tendency for the stream to dry out by early summer, flow can return at any later month given enough rainfall.

The geology of N04D is characterized by an alternating sequence of nearly horizontal 1–2 meter-thick limestone units (abbreviated ‘Ls’) and 2–4 meter-thick mudstones/shales¹¹ (Fig. 2). The bedrock is Permian in age and the lowest (oldest) unit is the Cottonwood Ls, which crops out immediately below the gauging station (Fig. 1). The next two limestone units above the Cottonwood Ls are the Morrill Ls and Eiss Ls, with the Eiss Ls being further divided into the hydrogeologically distinct Upper Eiss Ls and Lower Eiss Ls members. The Morrill Ls and Eiss Ls crop out in the lower fourth of the catchment and have 31 groundwater monitoring wells completed within them (Fig. 1). The wells are made of 5 cm diameter PVC and have 61 cm long screens with gravel packs and bentonite seals.¹¹ Pressure transducers are installed in three of the wells that are completed in the Morrill Ls, one of the wells completed in the Upper Eiss Ls, and one of the wells completed in the Lower Eiss Ls. These transducers record the pressure and temperature of water and atmosphere within the wells every five minutes, which we use to evaluate hydraulic heads.

Fig. 2 Konza Prairie bedrock units that crop out approximately within the lower two-thirds of watershed N04D. Wells included in this study are screened within the Upper Eiss Limestone Member and Lower Eiss Limestone Member of the Bader Formation, and the Morrill Limestone Member of the Beattie Formation, Council Grove Group, Permian System. Limestone and shale are abbreviated Ls and Sh, respectively, in the figure. Unit thickness ranges are shown in gray font beneath the name of each unit. This stratigraphic column is based on Fig. 2a of Macpherson.¹¹

Groundwater flow within these limestone aquifers is mostly accommodated by fractures and solution-enlarged pores, meaning that hydraulic conductivity within each unit likely has extreme spatial heterogeneity. Indeed, well productivity varies considerably from site to site within the same unit¹¹ and recent geophysical work has shown a large degree of heterogeneity in the distribution of water within each unit.¹⁴ Previous studies have measured the hydraulic conductivities of the Morrill and Eiss Limestones via rising-head slug tests in the wells and reported values of 10⁻⁵–10⁻³ m d⁻¹ for the Morrill Ls, 10⁻⁸–10⁻⁵ m d⁻¹ for the Lower Eiss Ls, and 10⁻⁵–10⁻⁴ m d⁻¹ for the Upper Eiss Ls.^19,20 These values fall in line with those typical of karstic and non-karstic limestones and dolomites,²⁷ but given the aforementioned heterogeneity they should not be taken as measurements that can be accurately generalized across the entire spatial extent of each limestone.

Streamflow in the merokarst of Konza Prairie is predominantly generated by groundwater discharge from the limestone units via springs along the hillslopes and within the streambed. Because the limestones are interbedded with low permeability mudstones/shales, vertical movement of water between the limestones is likely restricted and thus each limestone unit acts as a semi-isolated aquifer predominately hosting horizontal groundwater flow as shown in the Graphical Abstract. Lithological differences between these aquifers can cause groundwater contained within each unit to have unique chemical compositions and flow properties.^11,18 Surface runoff can also contribute to streamflow during storms, but it is likely limited in importance due to the presence of soil macropores and bedrock fractures which route precipitation into the subsurface instead.^28,29

While the above hydrogeologic properties form the basic understanding of streamflow generation at N04D, water flow paths through the catchment can in fact be highly variable both in space and time. For instance, whereas each limestone unit is generally thought of as an isolated aquifer with limited flow across its neighboring mudstones, recent studies have shown that significant vertical transport can occur in localized areas. Barry²⁰ first found evidence of this cross-unit transport when a tracer dye injected in the Morrill Ls appeared in significant concentrations within groundwater discharging from the Cottonwood Ls below it. They posited that this flow was accommodated by a buried collapse feature. Sullivan et al.¹⁴ followed up on this finding with a geophysical survey that revealed many such areas of enhanced vertical flow, suggesting that this phenomenon might be occurring extensively throughout the catchment. As further evidence of the complicated water flow paths at N04D, the stream gauge at the outlet of the watershed generally records zero-flow between July and February of the following year, but a portion of the stream above the gauge experiences flow nearly perennially.¹¹ This perennial portion begins where the outcropping unit in the stream transitions from Upper Eiss Ls to Lower Eiss Ls and continues until the outcropping unit transitions from the Stearns Shale to the Morrill Ls. These observations suggests that a temporary flow path occurs during this portion of the year in which water discharges from the Upper Eiss Ls, flows some distance in the channel, then reinfiltrates into the subsurface and recharges the Morrill Ls.

3. Methods

3.1 Stream and groundwater sampling

In order to characterize the short-term dynamics of streamflow generation in N04D, we used a high frequency sampling strategy for surface water, groundwater, and precipitation. The sampling period ran from April 15, 2021 to July 16, 2021 to coincide with the wet season when the stream had consistent flow.

Surface water was collected from the stream immediately downstream of the flume at the catchment outlet. The frequency at which we collected stream samples varied based on the conditions in the watershed: during low flow conditions towards the end of the sampling period this frequency was roughly weekly, during the peak of the wet season this frequency was roughly daily, and during selected storms this frequency was roughly hourly. We sampled hourly during five storms, which fell on May 8–9, May 26–27, June 11, June 24–25, and July 15. In addition to the stream, for each of these events we also collected precipitation with a ball-in-funnel type collector.³⁰ These events were chosen to represent a range of storm intensities and antecedent watershed conditions. Sampling for these events began at the onset of precipitation and continued until approximately two hours after precipitation ceased. In total, we collected 44 stream samples over the 93 days sampling period.

For each stream sample, we measured water temperature, pH, specific conductivity, and dissolved oxygen using field probes placed directly into the stream. Probes were recalibrated for all parameters before each sampling trip. Immediately upstream of the probes, we collected and filtered water through 0.45 μm syringe filters into two 30 mL high density polyethylene bottles for major cation and anion analysis, a 10 mL glass vial for stable water isotope analysis, and a 60 mL amber glass bottle for analysis of non-purgeable organic carbon (NPOC) and total dissolved nitrogen (TN). The filtration equipment and polyethylene bottles that we used were either brand new or were acid washed prior to sampling. Amber glass bottles were baked overnight at 120 °C prior to sampling. Every piece of sampling equipment was rinsed with sample water three times before being filled up with sample. Within a day of sampling, we preserved cation and NPOC/TN samples with three drops of concentrated trace element grade HNO₃ and analytical grade HCl, respectively. All samples were stored on ice during extended field work outings and at 2 °C in a refrigerator upon return to the laboratory.

For each of the hourly-sampled storms, we left out a rain collector over the entire duration of the storm and, when volumes were sufficient, divided the resulting precipitation into two 10 mL aliquots for replicate analysis of stable water isotopes.

We collected groundwater samples on six different days throughout the sampling period. For each groundwater sampling event, we collected water from the three wells in the 4-6 well nest (see Fig. 1 for well locations). Each of these wells is screened in one of the three limestone aquifers which crop out in the lower reach of the stream: the Morrill Ls (well name 4-6 MOR), Lower Eiss Ls (4-6 EIS 1), and Upper Eiss Ls (4-6 EIS 2). On three of the sampling trips, we also collected groundwater from the 3-5 and 3-5-1 wells which draw from the Morrill Ls (3-5 and 3-5-1 MOR). Before sampling each well, we purged a volume of water equal to two times the standing water volume via bailing. We then bailed water into two glass beakers, rinsing each beaker three times with sample water before filling them up for collection and analysis at each well. We used one beaker for field probe measurements and the other for collection into sample bottles, as described for stream samples. At the 3-5 and 3-5-1 MOR wells, we only collected water for analysis of major ions and not stable water isotopes, NPOC, or TN. The 4-6 well completed in the Lower Eiss Ls often did not yield enough water to fill all four collection bottles, and at these times, we also only collected water for major ion analysis.

3.2 Geochemical analyses

We analyzed all samples (stream, groundwater, and precipitation) within a month of collection. We measured alkalinity in surface and groundwater samples within a week of collection using the Gran alkalinity titration method with 0.02 N H₂SO₄ as the titrant. We measured major ion concentrations (Na⁺, NH₄⁺, K⁺, Mg²⁺, Ca²⁺, Sr²⁺, F⁻, Cl⁻, NO₂⁻, Br⁻, NO₃⁻, and SO₄²⁻) using Thermo Scientific ICS-1100 Ion Chromatographs and NPOC and TN using a Shimadzu TOC-L. For both instruments, quality control procedures included requiring a correlation coefficient of at least 0.995 for successful calibration and periodic analysis of calibration verification samples during each run. Precision and detection limits for these analyses are available in the ESI (Table SI1†).

The δ¹⁸O and δD of water isotope samples were measured at 0.1 and 0.5‰ precision using a Picarro L-i2130 water analyzer in the Stable Isotope Mass Spectrometry Laboratory at the Kansas State University Division of Biology. Results are presented as per mil (‰) relative to Vienna Standard Mean Ocean Water (VSMOW).

3.3 Existing geochemical data

To establish endmembers for the subsequent mixing analysis, we supplemented our collected data with groundwater and soil water chemistry from previous studies at Konza Prairie. We extended our groundwater dataset with data from Andrews,³¹ who collected groundwater samples from the same wells used in this study between June–December 2020, and with unpublished data from samples that we collected during routine groundwater sampling at the Konza Prairie between January–March 2021. In total, we added nine samples from 4-6 MOR, eight samples from both 4-6 EIS 1 and 4-6 EIS 2, and ten samples from both 3-5 MOR and 3-5-1 MOR. These data were collected using the same methods described above. This extension allowed us to better account for some natural variability in groundwater chemistry as we defined endmember ion concentrations.

To quantify contributions to streamflow from soil water, we relied on soil water chemistry data from Tsypin and Macpherson,²⁹ who collected eight samples from the A horizon (∼17 cm deep) and nine samples from the B horizon (∼150 cm deep) in N04D during the summer of 2010. We recognize that there is some uncertainty introduced by establishing an endmember from one season of field data, but we are not aware of other available soil water chemistry data sets for our study site. Nevertheless, the fact that both Tsypin and Macpherson²⁹ and the present study were conducted during the summer suggests that their soil data may be generally applicable for this analysis.

3.4 Hydrological and hydrogeological sensor data

Stream discharge at the outlet of N04D is measured at five-minute intervals at a triangular throated flume, and daily precipitation is measured at a rain gauge located immediately adjacent to the flume. These data are available on the online Konza Prairie data catalog (https://lter.konza.ksu.edu/data). Groundwater heads were evaluated using data from pressure transducers located in each of the five wells that we sampled for this study. Pressure data was recorded in five-minute intervals. We downloaded the data from the pressure transducers onto a datalogger during each sampling and calibrated the resulting groundwater head values to manual tape-down measurements of the water level performed in each well prior to every sampling. The calibrated data was manually inspected for continuity in groundwater head between the end of one downloaded data period and the beginning of the next. If there was a discontinuity, the entire data set in the later period was shifted so that its first data point was equal to the last data point in the earlier period.

3.5 Endmember mixing analysis

We employed endmember mixing analysis (EMMA) to assess the source contributions to streamflow and their short-term dynamics which might inform on the processes governing streamflow generation. EMMA is a core component of many studies on streamflow generation^32–34 because it allows for a quantitative assessment of the contributions to streamflow from any number of potential water sources using chemical tracers. Such results are achieved by solving the following generalized equations:

X_M = X_Af_A + X_Bf_B +…+ X_nf_n (1)

f_A + f_B +…+ f_n = 1 (2)

where X represents a chemical tracer such as a solute concentration, f is the fractional contribution of endmember n to the mixed water, subscript M denotes the mixed water, and other lettered subscripts (A, B, …n) denote the various endmembers.

Selecting which tracers to employ and which water sources to use as endmembers is the key consideration for mixing analysis. While endmembers can be statistically determined from the stream chemistry data alone,³⁵ they can also be chosen ahead of time based on previous understanding of the system in question.³⁶ Based on the work of Sullivan et al.,¹⁶ we assumed that N04D stream water is composed of precipitation (incorporating precipitation falling directly into the stream channel and surface runoff), soil water, and multiple groundwater endmembers, each sourced from different limestones.

Having more than three endmembers gives rise to significant mathematical complexity.³⁵ To avoid this complexity, we separated our mixing analysis into two stages. In the first stage, we used stable water isotope ratios (δ¹⁸O and δD) as tracers in a two-endmember mixing calculation that determined the fractional contributions of precipitation and subsurface water in streamflow, an analysis known as isotope hydrograph separation.³⁷ The δ¹⁸O and δD values used as the endmember compositions varied throughout the sampling period. For calculations at each stream sampling time point, we used the isotope ratios from the most recent precipitation sample as the precipitation endmember and the mean isotope ratios from the most recent groundwater samples as the subsurface endmember. We conducted the mixing calculations separately for the two isotopic tracers and then averaged the results. As a last step in this stage of the analysis, we calculated discharge-weighted averages of the fractional contribution of precipitation to streamflow across each individual storm and the entire study period to understand the overall importance of runoff to streamflow.

Some error is likely introduced to these calculations by the fact that the precipitation endmember is established from a single composite sample of precipitation collected over the course of a storm. Isotopic ratios in precipitation are known to vary throughout storms,³⁸ so our use of a single, volume-averaged endmember means our calculation of the fractional contribution of precipitation in streamflow is not entirely accurate on the hourly time scale of some of our stream samples. While this creates uncertainty in the distribution of precipitation inputs over time, we maintain that the proportions of precipitation in streamflow averaged across the duration of each storm are accurate because our precipitation end-member values were measured on composite samples. Given that the error introduced by within-storm variability in precipitation stable water isotopic ratios is typically only 0.5–2‰ for δ¹⁸O values,^39–41 and that previous studies have reported a general insignificance of surface runoff at N04D,^42,43 we reason that the uncertainty is low in the distribution of precipitation inputs over time and thus does not impact the conclusions we make about overall behaviors in the catchment.

After accounting for the contribution of precipitation to each stream sample using stable water isotopes, the second stage of our analysis was to then use major ion concentrations as tracers to calculate the contributions from the remaining three endmembers: the soil water and two groundwater sources. We first adjusted the ion concentrations in each stream sample to remove the influence of precipitation via a rearrangement of eqn (1):

where subscript s refers to the subsurface component and subscript p refers to the precipitation component. This step ensured that the proceeding mixing calculations were only separating streamflow based on inputs from the subsurface components. We then chose tracers for the calculations that differed significantly between endmembers but varied relatively little within each endmember. We verified our choices by plotting the mean concentration of tracer ions for each endmember on a bivariate plot and ensuring that they completely bounded all stream samples. After verification, we applied the endmember mixing equations (eqn (1) and (2)) to each stream sample to calculate the fractional contribution of each water source to streamflow throughout the study period.

4. Results and discussion

The overarching goal of this study is to advance understanding of the mechanisms driving streamflow generation in the merokarst environment of N04D. To do this, we sought to quantify the relative contributions to streamflow from N04D's multiple water sources and characterize their dynamics on short time scales. Below, we present the results of our sampling campaign and subsequent mixing analysis, then discuss what these results tell us about shifting flow pathways and the influence of groundwater storage thresholds on stream intermittency.

4.1 Sampling period precipitation and discharge patterns

N04D received a total of 620.5 mm of precipitation during 2021, which is lower than the annual mean of 843 mm (period of record: 1982–2021).²⁴ Of this total amount, 57% or 353.9 mm fell during the sampling period of this study (April 15–July 16). The precipitation regime during the study period was dominated by short, intense storms rather than drawn out periods of continuous rain (Fig. 3). 49.9% of the precipitation that fell during the sampling period is accounted for by just two storms, which occurred on May 16 and July 15.

Fig. 3 Daily total precipitation (black bars) and average discharge (blue line) in N04D for 2021. We added 0.01 to each discharge value so that no-flow days could be displayed on a log scale. Sampling period is outlined by dashed lines and hourly-sampled storms are marked with red arrows.

The flume at the catchment outlet recorded flow for most of the sampling period, with the exception of a dry period between July 5 and July 14 (Fig. 3). Flow was modest at the beginning of the sampling period, averaging 0.004 m³ s⁻¹ between April 15 and May 15. A large storm on May 16 (96.5 mm as measured at the catchment outlet) caused discharge to spike to its maximum for the year of 4.69 m³ s⁻¹, which then dropped to an average of 0.077 m³ s⁻¹ on the following day. An equivalent water depth of 39.6 mm exited the catchment via the stream in the 24 hours following this large storm, and 74.2 mm exited by the beginning of the next major storm on May 26. After the May 16 storm, discharge generally decreased over the remainder of the sampling period and dried out completely on July 5. Occasional storms, such as the one on May 26, brought small spikes in discharge during this drying period. After the stream completely dried out, the large storm on July 15 caused a week-long return of modest flow. Across the whole sampling period, 1.49 × 10⁵ m³ of water exited N04D via the stream, which equates to 124.2 mm of water depth over the whole catchment.

The storms that we selected for hourly sampling are marked in Fig. 3 and covered a range of hydrological conditions. Total precipitation that fell for these events ranged from 2.80 mm to 80.0 mm and their spacing throughout the sampling period meant that stream flow levels were different for each. One event coincided with the low-flow period before the large May 16 storm, three events occurred over the progressive drying period after the May 16 storm, and one event covered the large July 15 storm that produced flow after the prolonged dry period.

4.2 Sampling period groundwater head patterns

Fig. 4 shows how groundwater head in each of the sampled wells responded to hydrologic events during the study period. The 4-6 EIS 2 well (Upper Eiss Ls) responded rapidly to all significant storms (>10 mm precipitation) with near-instantaneous spikes in head that are indicative of a pressure response to the newly added precipitation. Modest spikes in head were produced by three moderate storms in March and April, but water levels in the well recessed back to their pre-storm levels prior to the onset of each subsequent storm. The large storm on May 16 caused a 1.33 meter increase in head which then recessed over the next several weeks. Storms occurring after May 16 produced smaller spikes in water level that interrupted the recession from the May 16 event, but the well returned to its pre-storm head on June 14.

Fig. 4 Groundwater head in meters above sea level (m.a.s.l.) within each of the sampled wells from February 23, 2021 to August 7, 2021. Daily total precipitation is also shown. Well IDs refer to locations in Fig. 1, with EIS 1 corresponding to the Lower Eiss Ls and EIS 2 corresponding to the Upper Eis LS. The vertical dashed lines outline the sampling period of this study. Sharp downward spikes in the groundwater head of 4-6 EIS 1 occurred in response to groundwater removal from the well during sampling. These spikes are largely not visible in the head profiles of the other wells, reflecting differences in hydraulic conductivities between limestone units.

The wells in the Morrill Ls (3-5 MOR, 3-5-1 MOR, and 4-6 MOR) all responded similarly to one another. As with 4-6 EIS 2, heads in each of these wells increased as a result of the three moderate storms in March and April. Unlike 4-6 EIS 2, however, these increases were relatively drawn out rather than instantaneous spikes, and the increases associated with each storm did not fully recess to pre-storm levels prior to the onset of the next storm. Heads in these wells were still in the process of recessing from these March and April storms when the large May 16 storm occurred. As a result of the May 16 storm, water levels instantaneously spiked upwards similarly to 4-6 EIS 2, but the magnitudes of their increases were not as large (0.938 m in 3-5 MOR, 0.563 m in 3-5-1 MOR, and 0.339 m in 4-6 MOR). Following this storm, heads resumed their drawn out recessions which continued for well over a month and did not land on low, stable values until late July or early August. As with head in 4-6 EIS 2, moderate storms which occurred during this recession phase produced small spikes in head, but in the Morrill wells these interruptions were lower in magnitude and were again followed by slower recessions.

These differences in groundwater head behavior indicate a difference in how water moves through the Upper Eiss Ls and Morrill Ls, despite previous slug tests resulting in similar hydraulic conductivities for them.^19,20 One possibility is that the slug tests do not capture the spatial heterogeneity in fractures and solution conduits that causes the Upper Eiss Ls to have an overall higher hydraulic conductivity than the Morrill Ls. This would explain why groundwater head fluctuations in March and April were flashier for the Upper Eiss Ls than the Morrill Ls, and why the recessions following these storms took longer in the Morrill Ls. This hypothesis is complicated, however, by the fact that wells in both the Upper Eiss Ls and the Morrill Ls seem to have similar increasing pressure responses to storms occurring after April. Although the spikes produced in the Morrill Ls wells were lower than those in 4-6 EIS 2, they were instantaneous spikes nonetheless and not the relatively drawn-out peaks like those in March and April. This similarity in response times to the onset of storms in May and beyond suggests that overall hydraulic conductivities in the two units might not be significantly different.

An alternative explanation is that the differences in behavior are caused by the simple stratigraphic order of the two units. The Upper Eiss Ls is positioned above the Morrill Ls, so it would naturally be more impacted by the addition of recharge from the surface. This relationship could create a stronger pressure response in the groundwater head in 4-6 EIS 2 than in the underlying Morrill Ls wells, consistent with observed head responses to storms throughout the study period. Similarly, the slower recessions for the Morrill wells can also be explained by stratigraphic order. While horizontal groundwater flow is favored in these limestone aquifers, the presence of spatially discrete vertical connections between different units causes higher units to ‘leak’ into lower ones.^16,20 Thus, the Morrill Ls might be receiving continuous input of groundwater from the overlying Upper Eiss Ls as it recesses following a storm, effectively extending the amount of time it takes to return to a low, stable head.

Compared to the wells discussed thus far, changes to groundwater head in the 4-6 EIS 1 well (Lower Eiss Ls) occurred very slowly. Whereas the wells in the Upper Eiss Ls and Morrill Ls had near instantaneous pressure responses following the May 16 storm, 4-6 EIS 1 did not achieve its maximum head until May 22, six days after the others. Further, the well took several days to recover the water that was removed during each sampling (Fig. 4). This is in accordance with previous studies concluding that the Lower Eiss Ls has a significantly lower conductivity than the adjacent limestone units.^11,19

4.3 Water chemistry

Stream water and groundwater samples had similar composition, with concentrations of every major ion in the stream samples falling within the range measured in the groundwater samples (Fig. 5, ESI Tables SI3 and SI4†). The most abundant solutes in both surface and groundwater samples were alkalinity, Ca²⁺, Mg²⁺, and SO₄²⁻, all of which except alkalinity show considerable variability between some of the wells. For Mg²⁺ and SO₄²⁻, the 3-5-1 MOR well had the highest concentrations of any well, and 4-6 EIS 2 had the lowest concentrations. Major ion concentrations in the 4-6 MOR and 3-5 MOR wells were very similar to each other, suggesting a common source for the water that makes its way into them. Previous studies have indicated this recharge source is the stream itself.^11,20 Differences in major ion concentrations between the groundwater from these wells and the stream water may reflect mixing relationships within the aquifer, mineral precipitation and dissolution, and differences in residence times between the groundwater and stream water.

Fig. 5 Box plots of select major ion concentrations in all surface and groundwater samples. Hinges mark the interquartile range of each group of samples, and stems extend to 1.5 times the interquartile range. Individual points mark outlier samples that had concentrations greater or less than 1.5 times the interquartile range. Extreme outlier points are not shown for K⁺ and Cl⁻ to show finer variation in concentrations. Complete data are available in the ESI (Tables SI3 and SI4†).

Some samples of the stream water contained Cl⁻, Na⁺, and K⁺ concentrations that are 3–50 times higher than the median (Fig. 5). These outliers all came from samples collected during the May 8–May 9 storm, which was the first major rain event to occur after N04D was burned on April 12. Previous studies in various grassland and non-grassland environments have shown that fires can mobilize large amounts of nutrients from soil and vegetation that can later be leached into surface waters.^44–46 Runoff during the May 8–May 9 storm could have carried these recently mobilized nutrients into the stream, offering one potential explanation for the spikes in Cl⁻, Na⁺, and K⁺ levels. Why other major ions such as Ca²⁺ and Mg²⁺ did not spike in concentration as well is unclear.

Values of pH were generally higher in stream water samples (avg. 7.69) than in the groundwater samples (avg. 6.99), and vice versa for alkalinity (ESI Tables SI3 and SI4†). This difference can be attributed to CO₂ outgassing from the water after it is discharged from the subsurface.⁴⁷ Stream water NPOC values averaged 3.00 mg L⁻¹ and did not vary significantly with discharge (ESI Table SI3†). NPOC was not measured in groundwater due to use of dye tracers for previous research.²⁰ For both surface water and groundwater samples, the concentrations of NO₂⁻, NO₃⁻, Br⁻, NH₄⁺, and TN were frequently below the detection limits of our methods and thus excluded from the remainder of the analysis.

4.4 Stable water isotopes

A well-fitting linear regression (r² = 0.97) defines the δ¹⁸O and δD values that we measured in our stream, groundwater, and precipitation samples (Fig. 6). This line has a slightly lower slope and intercept than the Global Meteoric Water Line (GMWL)⁴⁸ which is typical for the temperate, mid-continental climatic setting of the Konza Prairie⁴⁹ and is similar to a Local Meteoric Water Line (LMWL) calculated for the site (personal communication, Dr Jesse Nippert). The δ¹⁸O values of our groundwater samples ranged from −6 to −5.3‰ which is a similar range to those found in other studies of the groundwater in and around Konza.^17,50 The δ¹⁸O values of our precipitation samples were relatively enriched, ranging from −4.4 to 0‰. Seasonal variations in the isotopic compositions of Konza precipitation result in summer precipitation that is isotopically enriched, as our results show, and winter precipitation that is isotopically depleted, reaching δ¹⁸O values as low as −14‰.⁵¹ The fact that our groundwater samples are somewhere in the middle of this range suggests that the water stored in the limestone aquifers is an integration of summer and winter precipitation over multiple years.

Fig. 6 Stable water isotopes in stream water, groundwater, and precipitation samples collected at N04D. Solid line denotes the global meteoric water line (GMWL) (δD = 8 × δ¹⁸O + 10), and the dashed line denotes the linear regression of all plotted samples. Isotope data are available in the ESI (Tables SI3–SI5†).

The δ¹⁸O and δD values of the stream water samples generally fell between the relatively depleted values of the groundwater and the relatively enriched values of the precipitation, but positioned closer to the former than the latter (mean δ¹⁸O of the stream samples = −5.4‰). The isotopic ratios of the stream water only deflected closer to the precipitation values at the peak of major storms. The maximum δ¹⁸O value of our stream water samples was −4.4‰, which occurs at the peak of the May 8–May 9 storm. These results suggest that, during our sampling events, the stream water was comprised mostly of groundwater with somewhat minor amounts of more recent precipitation (i.e., surface runoff) mixed in, consistent with findings of Keen et al.¹⁷

4.5 Endmember mixing analysis

4.5.1 Endmember and tracer selection. The first stage of our mixing analysis was a two-endmember calculation to evaluate contributions from precipitation and subsurface water using stable water isotope ratios as tracers. Endmember isotope compositions for each of the five hourly-sampled storm events are provided in Table 1. The event endmember compositions used for each storm are the values measured from the corresponding composite precipitation sample. The pre-event endmember compositions for each storm are an average of values measured from the most recent groundwater samples from the 4-6 EIS 2 and 4-6 MOR wells, which were within 0.2‰ δ¹⁸O of each other for each sampling date. Discharge-weighted averages of the fractional contributions to streamflow from precipitation are included in Table 1 for each storm. Subsequent mixing analyses use the fractional contributions calculated for each hourly time point rather than this overall averaged value (ESI† Table SI6†). This approach carries uncertainties which are described in Section 3.5. Given the mostly low contribution of precipitation to streamflow, however, we reason that this uncertainty has a minimal impact on our subsequent analyses.

Table 1 Stable water isotope ratios used for endmember mixing between event and pre-event water. The date ranges correspond to the time periods around the storms identified in Fig. 3. All isotope ratios are expressed per mil (‰) VSMOW

Date range		Event endmember		Pre-event endmember		Fraction event
Start	End	δ ^18O	δD	δ ^18O	δD	—
5/6	5/14	−1.8	−6.5	−5.7	−36	0.181
5/22	6/1	−1.25	−10	−5.7	−35	0.018
6/8	6/14	0	12	−5.6	−35	0.002
6/22	6/27	−1.3	−7	−5.6	−35	0.021
7/15	7/16	−4.4	−22	−5.6	−36	0.500

---

The second stage of the mixing analysis was a three-endmember calculation between contributions from soil water, groundwater from the Morrill Ls as represented in the 3-5-1 MOR well, and groundwater from the Upper Eiss Ls as represented in the 4-6 EIS 2 well (Fig. 7). We selected 3-5-1 MOR to represent the Morrill Ls because it had a stable chemistry over the entire sampling period, displaying little variation in all chemical parameters (Fig. 5 and ESI Table SI4†). The 3-5 MOR and 4-6 MOR wells receive significant amounts of recharge directly from the stream,^11,20 making them unsuitable for stream water mixing calculations and leaving 3-5-1 MOR as the better representation of ‘pure’ groundwater held in the Morrill Ls. We did not include the Lower Eiss Ls as an endmember because it has a much lower hydraulic conductivity than the other aquifers and has been shown to have little interaction with the stream.^11,19

Fig. 7 Bivariate plot of Ca²⁺/Mg²⁺vs. SO₄²⁻ for the three endmembers selected for mixing analysis, with stream water sample chemistries adjusted for the influence of precipitation as described in the text. The average values of Ca²⁺/Mg²⁺ and SO₄²⁻ that were used as the endmember compositions are marked with red diamonds, and error bars show one standard deviation from the mean. Dashed lines denote the mixing envelope. Endmember compositions (Table SI2) and an identical plot with raw stream chemistry data (Fig. SI2) are available in the ESI.†

To calculate the fractional contributions from the endmembers, we used SO₄²⁻ concentration and the mass ratio of Ca²⁺/Mg²⁺ as tracers. We selected these tracers based on the criteria laid out in the methods: SO₄²⁻ and Mg²⁺ concentrations both displayed high variation between the groundwater wells selected as endmembers and low variation within each well (Fig. 5). Using the concentration ratio Ca²⁺/Mg²⁺ instead of the concentration of either ion made the tracer more robust to the effects of evaporation, sorption, and mineral precipitation in the stream channel. Indeed, long-term data from N04D has shown that there is little variation in Ca²⁺/Mg²⁺ values between stream samples collected on the same day from different points along the stream reach.⁵² Because we already accounted for the contributions to streamflow from precipitation in the first stage of the analysis, we adjusted the SO₄²⁻ and Ca²⁺/Mg²⁺ values of all stream samples using eqn (3) to remove the influence of precipitation on these concentrations. All adjusted stream values were successfully bounded by a mixing envelope created by our selected endmembers in a SO₄²⁻vs. Ca²⁺/Mg²⁺ plot (Fig. 7). Endmember compositions (Table SI2†) and an identical plot without adjusted stream compositions (Fig. SI2) are available in the ESI.†

A major assumption that we made in the mixing analysis of this study is that groundwater contributions to the stream were only coming from the Upper Eiss Ls and the Morrill Ls. While we provided justification above for excluding the Lower Eiss Ls, N04D contains a total of six different limestone formations that outcrop above the flume,¹⁴ each of which could potentially contribute water to streamflow. However, the Eiss and Morrill members are the stratigraphically lowest of these limestones and have much more contact area with the stream channel than the higher units in the steeper upper reaches of the catchment. Several dozen springs drain the limestone units of N04D²⁰ primarily in the upper reaches of the watershed, providing a potential path for water to reach the stream from these upper limestones without needing direct contact with the channel. However, these springs only produce water ephemerally and, based on long-term monitoring of stream and groundwater chemistry at the site,⁵² their contribution to streamflow is likely very small. Specifically, the stream water chemistry adjacent to the Upper Eiss Ls is generally very similar to the groundwater from the Upper Eiss Ls in terms of major ions and temperature. Thus, we reason that much of the water from the upper limestone units likely makes its way into the Eiss or Morrill members eventually.

4.5.2 Source contributions to streamflow. The immediate aim of this study was to quantify source contributions to streamflow at N04D and characterize their short-term dynamics. The results of our mixing calculations show that source contributions varied in relative strength over the sampling period but were overwhelmingly dominated by groundwater sources (Fig. 8). Stream discharge summed across all of our sampling time points partitions into 96.3% groundwater (73.5% from the Upper Eiss Ls and 22.9% from the Morrill Ls), 3.83% surface runoff, and 0.13% soil water. Our discrete sampling points did not span every storm or moment during the study period, so the true values of these percentages would likely differ somewhat from our results. However, these results are representative of overall relationships in the stream because we distributed our sampling across the entire wet season and included a wide range of hydrological conditions.

Fig. 8 Fractional contributions to streamflow from precipitation, the two groundwater endmembers, and soil water through the sampling period. Stream discharge is marked with the light blue line for reference to flow conditions. Dotted lines indicate extended periods of time between samplings when short-term dynamics are not captured. Mixing calculation data are available in the ESI (Table SI6†).

While groundwater dominated streamflow for most of the sampling period, our mixing calculations indicate that contributions from each aquifer varied in relative importance through the study period (Fig. 8). The Upper Eiss Ls contributed the most water to streamflow among all sources for most of the sampling period, especially during high flow conditions (daily average stream discharge > 0.02 m³ s⁻¹) when its fractional contribution could exceed 0.75. Given that the Upper Eiss Ls is stratigraphically higher and outcrops over a larger surface area than the Morrill Ls, it follows that it generally contributed more to streamflow across the sampling period. However, as the stream dried following the May 16 storm, contributions from the Morrill Ls slowly overtook those from the Upper Eiss Ls and exceeded them by late June when daily average stream discharge had fallen below 0.01 m³ s⁻¹.

This pattern of dominance by the Upper Eiss Ls during high flow conditions and by the Morrill Ls during prolonged dry periods might be related to the differences between the units in pressure response and recession time that we discussed in Section 4.2 (Fig. 4). The Upper Eiss Ls generally has higher magnitude spikes in groundwater head and more rapid rates of recovery towards pre-storm levels than the Morrill Ls, which we hypothesized might be due to enhanced fracturing and dissolution in the Upper Eiss Ls, giving it a greater overall hydraulic conductivity. According to this hypothesis, the Upper Eiss Ls dominates streamflow under wet conditions by quickly routing recharge into the channel, but in the absence of further additions of precipitation this output quickly dwindles. In comparison, the slightly slower responding Morrill Ls releases its water in more of a steady trickle, which has a lower volumetric flow rate than the discharge from the Upper Eiss Ls but can be sustained for longer and thus becomes proportionally more important as the Upper Eiss Ls depletes.

Even if the overall hydraulic conductivities of the two limestone units are not significantly different, this pattern in their fractional contributions can also be explained by their stratigraphic order. As the stratigraphically lower unit, much of the storm response in the Morrill Ls relies on vertical transport of groundwater from the overlying Upper Eiss Ls and through the Stearns Shale that separates them. Though some spatially discrete structures exist throughout the catchment and may be able to accomodate rapid vertical transport,^14,20 vertical flow between limestone units will likely still be overall slower than horizontal flow within each individual limestone unit. Thus, when a storm creates a pressure response in the groundwater within the Upper Eiss Ls, much of the resulting flow will be directed horizontally towards the stream and relatively less will move vertically into the Morrill Ls. As a result, discharge from the Upper Eiss Ls will be more important to streamflow than discharge from the Morrill Ls in the immediate aftermath of storms. During dry periods, however, the slower vertical movement of water into the Morrill Ls sustains an elevated head for a longer period of time, allowing it to continue discharging into the stream after the Upper Eiss Ls has drained. The net results of either hydraulic conductivity or stratigraphic order controlling the observed behavior in streamflow generation are the same, but further studies would be helpful in evaluating the extent to which either hypothesis is supported.

Surface runoff only comprised a significant proportion of the total stream discharge during the peak of storms occurring when the stream was dry (July 15 storm) or under low-flow conditions (May 8–May 9 storm) (Fig. 8). The May 8–May 9 and May 26–May 27 storms had very similar amounts of precipitation (22.1 mm and 22.5 mm, respectively), but were comprised of an order of magnitude difference in proportions of surface runoff (0.181 and 0.018). The difference between these storms was the flow conditions in the stream prior to precipitation. In the May 8–May 9 storm, the stream was nearly dry (discharge = 0.0016 m³ s⁻¹) and in the May 26–May 27 storm flow was over an order of magnitude larger (discharge = 0.0272 m³ s⁻¹). Thus, the importance of surface runoff in the stream compared to inputs from subsurface sources seems to diminish drastically during high flow conditions.

We hypothesize that this behavior reflects a ‘fill and spill’ dynamic to the hydrology of N04D: when the watershed is dry, precipitation infiltrating into the subsurface merely ‘fills’ the subsurface storage units but does not yet push groundwater into the stream, causing the precipitation that does run over the surface to comprise a relatively large proportion of streamflow. Once these subsurface units surpass some critical thresholds of storage, however, they ‘spill’ into the stream and dwarf the surface runoff in terms of fractional contribution.

Fill and spill models were originally conceived to explain subsurface stormflow on the hillslope scale,⁵³ but subsequent studies have shown that fill and spill dynamics can operate across entire watersheds and are heavily controlled by sediment and rock properties.^54,55 Gutierrez-Jurado et al.⁸ modeled streamflow generation in a non-karstic catchment and showed that the gradual filling of groundwater storage zones can be a major controller on flow intermittency. The hydrogeologic regime of merokarst terrains is very different than that of the sandy aquifer that Gutiérrez-Jurado et al. modeled. However, Costigan et al.⁹ hypothesized that such a fill and spill dynamic was operating in N04D and concluded that streamflow and precipitation only became synchronized once the watershed attained a threshold amount of subsurface water storage. Our results on the dynamics of surface runoff agree that fill and spill behavior characterizes streamflow generation at N04D, and further analyses on groundwater inputs to the stream as discussed below seem to add support to this notion.

4.5.3 Relationships with stream discharge, groundwater head, and storage thresholds. Well-defined relationships between fractional contributions and stream discharge (F–Q) build confidence in the selection of endmembers for the EMMA and can also be used to further understand the hydrological connections between the endmembers and the stream.⁵⁶ F–Q relationships for both the Upper Eiss Ls and the Morrill Ls are defined by a well-fitting (r² > 0.7) power law function, but only for samples that were collected after the major storm on May 16 (Fig. 9A and B). The high fitness of the F–Q relationships for the samples collected after May 16 support our conclusions about the effect of differences in hydraulic conductivity and/or stratigraphic position on streamflow generation from groundwater sources, as contributions from the Upper Eiss Ls are more important during wetter conditions and contributions from the Morrill Ls are more important during drier conditions. However, the F–Q relationships fail to capture the samples collected prior to May 16, suggesting that our selection of endmembers may not be appropriate for this period and potentially pointing to threshold fill and spill behavior at N04D that controls the connection between the groundwater and the stream.

Fig. 9 Relationship between fractional contribution to streamflow and stream discharge for the (A) Upper Eiss Ls and (B) Morrill Ls. The dotted lines denote the power law relationship fit to only the samples collected after May 16. Grey triangles denote samples which were collected prior to May 14.

The SO₄²⁻ and Mg²⁺ concentrations in groundwater samples from the 4-6 EIS 2 well in early May were over one standard deviation higher than their long-term averages, potentially reflecting an increased residence time of groundwater in the subsurface. SO₄²⁻ concentration may have increased in response to dissolution of gypsum or anhydrite whereas Mg²⁺ levels may have been high as a result of exchange with Ca²⁺ in clay minerals or potentially chlorite dissolution.^11,57 The lack of increased Ca²⁺ concentrations in these same samples is likely due to the fact that the groundwater generally equilibrated with calcite, so an increase in contact time will not result in further dissolution without some accompanying change in temperature or pH.

Regardless of the mineral sources of these ions, the increased contact time with bedrock for these early May samples implies that the Upper Eiss groundwater was not being extensively flushed into the stream at this time. This implied lack of a hydrologic connection suggests that the groundwater units had not yet reached the critical storage threshold required to ‘spill’ water into the stream. As such, this might explain why our calculated fractional contributions for this early period are different from what is predicted in the F–Q relationships and it also means that the results of our mixing calculations (Fig. 8) might not be accurate for this same period. The May 16 storm input a massive amount of water into the system, however, making it likely that the subsurface units then reached their storage thresholds, fully connected to the stream, and were able to spill water into it according to the established F–Q relationships. Based on these results, we interpret that fractional contributions to streamflow from the groundwater units at N04D can be reasonably estimated from the F–Q relationships so long as the aquifers have surpassed some threshold water storage amount.

Establishing relationships between groundwater contributions and hydraulic head in each limestone unit allows us to better understand the threshold that governs subsurface connection to the stream. In the Upper Eiss Ls, volumetric discharge of water to the stream (fractional contribution multiplied by stream discharge) varies linearly with groundwater hydraulic head in the 4-6 EIS 2 well (Fig. 10A). In the Morrill Ls, a linear relationship between discharge into the stream and hydraulic head in 3-5-1 MOR only holds for samples collected after the large May 16 storm (Fig. 10B). The samples collected prior to May 16 all plot well below the linear trendline of the post-May 16 samples, indicating that discharge into the stream is much lower than predicted by the relatively high groundwater heads in 3-5-1 MOR at this time.

Fig. 10 Relationship between volumetric input to the stream of water from the (A) Upper Eiss Ls and (B) Morrill Ls, and head in their respective wells. Dotted regression lines are fit only to the samples collected after May 16. Grey triangles denote samples which were collected prior to May 14.

One interpretation of this discrepancy is that it is further evidence of a disconnect between the groundwater and the stream during these earlier sampling points. After the May 16 storm, the behavior switches to the predictable linear relationship, possibly indicating that a storage threshold was achieved. The clear difference in fractional contribution values for the samples collected from the two periods might also be related to the fact that our mixing model may be inaccurate for the pre-May 16 samples, but in this case, we would also expect there to be a similar difference in the samples for the Upper Eiss Ls. The lack of a similar behavior switch in the Upper Eiss Ls might instead suggest that the threshold dynamics are different for these two aquifers, and that the Upper Eiss Ls requires much less addition of water to connect to the stream than the Morrill Ls does. A higher overall hydraulic conductivity in the Upper Eiss Ls or its higher stratigraphic position with respect to the Morrill Ls might cause this difference in threshold behavior, much as it causes the differences in hydraulic head responses to storms as discussed previously.

Another interpretation of the discrepancy between discharge into the stream and well hydraulic head for the Morril Ls is that head, or at least head in a single well such as 3-5-1 MOR, is not a suitable measure for determining if the storage threshold has been attained in the aquifers. The wells only measure hydraulic head at a single point for each limestone unit which can crop out along several dozens of meters of stream reach. The extreme spatial heterogeneity in fractures and solution conduits in these limestones means that head measured at one point may be significantly different from head measured at another point which exists along a more conductive flow path. Studies employing greater well control or geophysical techniques to constrain groundwater head across a greater spatial extent could be used to quantify the storage thresholds that control groundwater connection to the stream.

4.6 Implications

Large portions of Kansas, Oklahoma, Missouri, and Texas in the USA are covered by merokarst terrain,¹¹ yet streamflow generation processes which control the availability of surface water in these areas remain poorly understood. Our results allow us to hypothesize some of the main hydrological processes that determine amounts and timings of streamflow in N04D which can be applied to other merokarst headwater catchments to predict future streamflow behavior.

We interpret that groundwater discharge dominates headwater streamflow generation in the merokarst of N04D, during both high and low flow conditions (Fig. 8). Contributions from surface runoff are only important during some storms, and contributions from soil water are volumetrically insignificant across the entire wet season. These behaviors are in contrast to the widespread shallow-and-deep hypothesis of streamflow generation, which states that contributions to streamflow from shallow units (i.e. surface runoff and soil water) become increasingly more important relative to deep units (i.e. groundwater) under wet conditions.^58,59 While studies conducted at non-karstic catchments of similar size to N04D generally support the shallow-and-deep hypothesis,^60,61 we hold that the karstic properties of N04D prevent the development of such a regime. Thin soils place an upper limit on the significance of shallow flow paths, while soil macropores and bedrock fractures rapidly route water to limestone aquifers ensuring that deep flow paths remain dominant even during very wet conditions. Given that groundwater flow paths are critical in sustaining streamflow, our finding that there are groundwater storage thresholds which govern subsurface connections to the stream then implies that these thresholds are major controllers on stream intermittency. Significant, sustained streamflow can only occur when enough precipitation has fallen in a short enough amount of time to surpass the groundwater storage thresholds.

The thinly bedded nature of the limestones at N04D means that its hydrological behavior also deviates from what is typically seen in massive karst systems. Confining mudstone layers largely restrict vertical movement of water between limestones and the unique depositional history of each individual limestone bed means that each one might possess unique hydrogeological properties. Whether caused by differences in hydraulic conductivity or simple stratigraphic position, our results show that each of our studied limestone aquifers, which are vertically separated by only a couple of meters from one another, respond to recharge and interact with the stream in different ways (Fig. 4). This heterogeneity expands the capacity of the subsurface to act as a hydrological buffer, in that each limestone unit releases recharge into the stream at a different rate to effectively extend the amount of time that the stream can flow. Consistent with this observation, streams across the Flint Hills region of Kansas have greater proportions of baseflow in total streamflow than streams in adjacent, non-merokarstic regions,⁶² suggesting that the merokarst bedrock does improve hydrological buffering. Because the sensitivity of a catchment to different lengths of drought is strongly influenced by the hydraulic conductivity of the underlying rocks/sediment,⁶³ the existence of a spectrum of conductivities might make merokarst terrains such as N04D more resilient to drought overall. This effect should not be overstated, as flow at N04D is still highly intermittent, but it might act to mitigate some of the risk that climate change poses to intermittent streams worldwide.^6,64

6. Conclusions

The processes governing headwater streamflow generation are poorly understood despite their critical importance in planning effective ecosystem management strategies.^3,7,8 In this study, we used mixing calculations and high-frequency sampling of stream and groundwater at the Konza Prairie Biological Station, USA to quantify source contributions to streamflow at an intermittent merokarst catchment and identify controls on stream intermittency. Our results indicate: (1) streamflow is overwhelmingly dominated by inputs from limestone aquifers, although surface runoff can be locally important during large storms that occur after prolonged dry periods, (2) the stream channel does not become fully connected to its groundwater sources until some threshold of groundwater storage in the subsurface is surpassed, and (3) contributions to streamflow from each aquifer vary through time based on the unique hydrogeological properties and/or stratigraphic position of each limestone. Because our study indicates both that groundwater is the dominant source of streamflow and that its connection to the stream is controlled by a water storage threshold, we conclude that this threshold dynamic is a major controller of stream intermittency. Additionally, we hypothesize that the merokarst terrain itself expands hydrological buffering in the system. The presence of many thinly bedded and heterogenous strata gives the catchment a higher resilience to drought than if it had a more homogeneous lithology. These findings are important to consider for future water management in merokarst regions such as those in the central United States, where changing climate threatens the ecosystem health of many intermittent headwater streams.

Author contributions

CH and MK designed the study and CH prepared the initial manuscript draft. All authors contributed to data collection, analysis, and manuscript revision.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We thank the Konza Prairie Biological Station and its staff for site access and management. We are also grateful for constructive conversations with Gwendolyn Macpherson, Pamela Sullivan, Rachel Keen, Walter Dodds, Benjamin Sikes, Erin Seybold, Samuel Zipper, and Logan Swenson, and for helpful input from two anonymous reviewers. This project was supported with funding from the National Science Foundation (DEB-2025849; EPSCoR OIA-1656006) and the State of Kansas.

References

1. L. B. Leopold, M. G. Wolman and J. P. Miller, Fluvial processes in geomorphology, Dover Publications, Inc., New York, NY, 1964.
2. R. J. Naiman, H. Décamps and M. E. McClain, Riparia: Ecology, Conservation, and Management of Streamside Communities, Elsevier, Burlington, MA, 1st edn, 2005.
3. T. Datry, S. T. Larned and K. Tockner, Intermittent Rivers: A Challenge for Freshwater Ecology, BioScience, 2014, 64, 229–235.
4. S. Tooth, Process, form and change in dryland rivers: a review of recent research, Earth-Sci. Rev., 2000, 51, 67–107.
5. L. Benda, M. A. Hassan, M. Church and C. L. May, Geomorphology of the steepland headwaters: the transition from hillslopes to channels, J. Am. Water Resour. Assoc., 2005, 41, 835–851.
6. K. L. Jaeger, J. D. Olden and N. A. Pelland, Climate change poised to threaten hydrologic connectivity and endemic fishes in dryland streams, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 13894–13899.
7. M. Shanafield, S. A. Bourke, M. A. Zimmer and K. H. Costigan, An overview of the hydrology of non-perennial rivers and streams, Wiley Interdiscip. Rev.: Water, 2021, 8, e1504.
8. K. Y. Gutiérrez-Jurado, D. Partington, O. Batelaan, P. Cook and M. Shanafield, What Triggers Streamflow for Intermittent Rivers and Ephemeral Streams in Low-Gradient Catchments in Mediterranean Climates, Water Resour. Res., 2019, 55, 9926–9946.
9. K. H. Costigan, M. D. Daniels and W. K. Dodds, Fundamental spatial and temporal disconnections in the hydrology of an intermittent prairie headwater network, J. Hydrol., 2015, 522, 305–316.
10. K. Shook and J. Pomeroy, Changes in the hydrological character of rainfall on the Canadian prairies, Hydrol. Processes, 2012, 26, 1752–1766.
11. G. L. Macpherson, Hydrogeology of thin limestones: The Konza Prairie Long-Term Ecological Research Site, Northeastern Kansas, J. Hydrol., 1996, 186, 191–228.
12. D. J. Weary and D. H. Doctor, Karst in the United States: A Digital Map Compilation and Database, USGS, Reston, VA, 2014.
13. S. E. Vero, G. L. Macpherson, P. L. Sullivan, A. E. Brookfield, J. B. Nippert, M. F. Kirk, S. Datta and P. Kempton, Developing a Conceptual Framework of Landscape and Hydrology on Tallgrass Prairie: A Critical Zone Approach, Vadose Zone J., 2018, 17, 11.
14. P. L. Sullivan, C. Zhang, M. Behm, F. Zhang and G. L. Macpherson, Toward a new conceptual model for groundwater flow in merokarst systems: Insights from multiple geophysical approaches, Hydrol. Processes, 2020, 34, 4697–4711.
15. D. C. Ford and P. Williams, Karst Hydrogeology and Geomorphology, John Wiley & Sons, Ltd, Hoboken, NJ, 1st edn, 2007.
16. P. L. Sullivan, M. W. Stops, G. L. Macpherson, L. Li, D. R. Hirmas and W. K. Dodds, How landscape heterogeneity governs stream water concentration-discharge behavior in carbonate terrains (Konza Prairie, USA), Chem. Geol., 2019, 527, 118989.
17. R. M. Keen, J. B. Nippert, P. L. Sullivan, Z. Ratajczak, B. Ritchey, K. O’Keefe and W. K. Dodds, Impacts of Riparian and Non-riparian Woody Encroachment on Tallgrass Prairie Ecohydrology, Ecosystems, 2022 DOI:10.1007/s10021-022-00756-7.
18. H. K. Wood and G. L. Macpherson, Sources of Sr and implications for weathering of limestone under tallgrass prairie, northeastern Kansas, Appl. Geochem., 2005, 20, 2325–2342.
19. M. L. Pomes, PhD thesis, University of Kansas, 1995.
20. E. R. Barry, MS thesis, University of Kansas, 2018.
21. J. M. Briggs, A. K. Knapp, J. M. Blair, J. L. Heisler, G. A. Hoch, M. S. Lett and J. K. McCarron, An ecosystem in transition. Causes and consequences of the conversion of mesic grassland to shrubland, Bioscience, 2005, 55, 243–254.
22. A. M. Veach, W. K. Dodds and A. Skibbe, Correction: Fire and Grazing Influences on Rates of Riparian Woody Plant Expansion along Grassland Streams, PLoS One, 2014, 10, e0129409.
23. B. P. Hayden, in Grassland Dynamics, ed. A. K. Knapp, J. M. Briggs, D. C. Hartnett and S. L. Collins, Oxford University Pres, New York, 1998, pp. 19–34.
24. J. Nippert, APT01 Daily precipitation amounts measured at multiple sites across konza prairie,  DOI:10.6073/pasta/dc20ec67bc74157cd37ae7989f44cc4c.
25. W. Dodds, ASD02 Stream discharge measured at the flumes on watershed N04D at Konza Prairie,  DOI:10.6073/pasta/14bad446298b9892f27ab9fc1b1dfd55.
26. K. O'Keefe, D. M. Bell, K. A. McCulloh and J. B. Nippert, Bridging the Flux Gap: Sap Flow Measurements Reveal Species-Specific Patterns of Water Use in a Tallgrass Prairie, J. Geophys. Res.: Biogeosci., 2020, 125, e2019JG005446.
27. M. A. Lewis, C. S. Cheney and B. É. Ó Dochartaigh, Guide to Permeability Indices, British Geological Survey, 2006, p. 29, (CR/06/160N), unpublished, https://nora.nerc.ac.uk/id/eprint/7457/.
28. G. L. Macpherson, J. A. Roberts, J. M. Blair, M. A. Townsend, D. A. Fowle and K. R. Beisner, Increasing shallow groundwater CO2 and limestone weathering, Konza Prairie, USA, Geochim. Cosmochim. Acta, 2008, 72, 5581–5599.
29. M. Tsypin and G. L. Macpherson, The effect of precipitation events on inorganic carbon in soil and shallow groundwater, Konza Prairie LTER Site, NE Kansas, USA, Appl. Geochem., 2012, 27, 2356–2369.
30. U. Prechsl, A. Gilgen, A. Kahmen and N. Buchmann, Reliability and quality of water isotope data collected with a low-budget rain collector, Rapid Commun. Mass Spectrom., 2014, 28, 879–885.
31. K. Andrews, MS thesis, Kansas State University, 2021.
32. R. M. Cowie, J. F. Knowles, K. R. Dailey, M. W. Williams, T. J. Mills and N. P. Molotch, Sources of streamflow along a headwater catchment elevational gradient, J. Hydrol., 2017, 549, 163–178.
33. R. Dwivedi, T. Meixner, J. C. McIntosh, P. A. T. Ferré, C. J. Eastoe, G. Niu, R. L. Minor, G. A. Barron-Gafford and J. Chorover, Hydrologic functioning of the deep critical zone and contributions to streamflow in a high-elevation catchment: Testing of multiple conceptual models, Hydrol. Processes, 2019, 33, 476–494.
34. F. Liu, M. W. Williams and N. Caine, Source waters and flow paths in an alpine catchment, Colorado Front Range, United States, Water Resour. Res., 2004, 40, W09401.
35. N. Christophersen and R. P. Hooper, Multivariate analysis of stream water chemical data: The use of principal components analysis for the end-member mixing problem, Water Resour. Res., 1992, 28, 99–107.
36. S. Inamdar, in Ecological Studies, Dordrecht, 2011, pp. 163–183.
37. J. Klaus and J. J. McDonnell, Hydrograph separation using stable isotopes: Review and evaluation, J. Hydrol., 2013, 505, 47–64.
38. J. J. McDonnell, M. Bonell, M. K. Stewart and A. J. Pearce, Deuterium variations in storm rainfall: Implications for stream hydrograph separation, Water Resour. Res., 1990, 26, 455–458.
39. S. Uhlenbrook and S. Hoeg, Quantifying uncertainties in tracer-based hydrograph separations: a case study for two-, three- and five-component hydrograph separations in a mountainous catchment, Hydrol. Processes, 2003, 17, 431–453.
40. M. R. McHale, J. J. McDonnell, M. J. Mitchell and C. P. Cirmo, A field-based study of soil water and groundwater nitrate release in an Adirondack forested watershed, Water Resour. Res., 2002, 38, 2-1.
41. S. K. Carey and W. L. Quinton, Evaluating runoff generation during summer using hydrometric, stable isotope and hydrochemical methods in a discontinuous permafrost alpine catchment, Hydrol. Processes, 2005, 19, 95–114.
42. L. J. Gray, G. L. Macpherson, J. K. Koelliker, W. K. Dodds, A. K. Knapp, J. M. Briggs, D. C. Hartnett and S. L. Collins, Hydrology and aquatic chemistry, Grassl. Dyn. Long-Term Ecol. Res. Tallgrass Prairie, Oxf. Univ. Press N. Y., 1998, pp. 159–176.
43. D. R. Steward, X. Yang, S. Y. Lauwo, S. A. Staggenborg, G. L. Macpherson and S. M. Welch, From precipitation to groundwater baseflow in a native prairie ecosystem: a regional study of the Konza LTER in the Flint Hills of Kansas, USA, Hydrol. Earth Syst. Sci., 2011, 15, 3181–3194.
44. C. Swindle, P. Shankin-Clarke, M. Meyerhof, J. Carlson and J. Melack, Effects of Wildfires and Ash Leaching on Stream Chemistry in the Santa Ynez Mountains of Southern California, Water, 2021, 13, 2402.
45. D. M. Larson, W. K. Dodds, M. R. Whiles, J. N. Fulgoni, T. R. Thompson and Y. Cao, A before-and-after assessment of patch-burn grazing and riparian fencing along headwater streams, J. Appl. Ecol., 2016, 53, 1543–1553.
46. P. Pereira, A. Cerda, D. Martin, X. Úbeda, D. Depellegrin, A. Novara, J. F. Martínez-Murillo, E. C. Brevik, O. Menshov, J. R. Comino and J. Miesel, Short-term low-severity spring grassland fire impacts on soil extractable elements and soil ratios in Lithuania, Sci. Total Environ., 2017, 578, 469–475.
47. K. S. Ohmes, G. Macpherson and B. L. Huff, Low stream-flow measurement and stream CO₂ emission calculation, Konza Prairie LTER site, Northeastern Kansas, USA, Geol. Soc. Am. Abstr. Prog., 2009, 41, 664.
48. H. Craig, Isotopic variations in meteoric waters, Science, 1961, 133, 1702–1703.
49. A. L. Putman, R. P. Fiorella, G. J. Bowen and Z. Cai, A Global Perspective on Local Meteoric Water Lines: Meta-analytic Insight Into Fundamental Controls and Practical Constraints, Water Resour. Res., 2019, 55, 6896–6910.
50. N. Clauer, I. Techer and S. Chaudhuri, Geochemical Tracing of Potential Hydraulic Connections between Groundwater and Run-Off Water in Northeastern Kansas, USA, Hydrology, 2017, 4, 56.
51. J. B. Nippert and A. K. Knapp, Linking water uptake with rooting patterns in grassland species, Oecologia, 2007, 153, 261–272.
52. G. L. Macpherson, AGW01 Long-term measurement of groundwater physical and chemical properties from wells on watershed N04D at Konza, Prairie, 2021.
53. H. J. Tromp-van Meerveld and J. J. McDonnell, Threshold relations in subsurface stormflow: 2. The fill and spill hypothesis, Water Resour. Res., 2006, 42, W02411.
54. L. Guo, H. Lin, B. Fan, J. Nyquist, L. Toran and G. J. Mount, Preferential flow through shallow fractured bedrock and a 3D fill-and-spill model of hillslope subsurface hydrology, J. Hydrol., 2019, 576, 430–442.
55. T. Sayama, J. J. McDonnell, A. Dhakal and K. Sullivan, How much water can a watershed store?, Hydrol. Processes, 2011, 25, 3899–3908.
56. F. Liu, M. H. Conklin and G. D. Shaw, Insights into hydrologic and hydrochemical processes based on concentration-discharge and end-member mixing analyses in the mid-Merced River Basin, Sierra Nevada, California, Water Resour. Res., 2017, 53, 832–850.
57. G. L. Macpherson and P. L. Sullivan, Watershed-scale chemical weathering in a merokarst terrain, northeastern Kansas, USA, Chem. Geol., 2019, 527, 118988.
58. W. Zhi and L. Li, The Shallow and Deep Hypothesis: Subsurface Vertical Chemical Contrasts Shape Nitrate Export Patterns from Different Land Uses, Environ. Sci. Technol., 2020, 54, 11915–11928.
59. W. Zhi, L. Li, W. Dong, W. Brown, J. Kaye, C. Steefel and K. H. Williams, Distinct Source Water Chemistry Shapes Contrasting Concentration-Discharge Patterns, Water Resour. Res., 2019, 55, 4233–4251.
60. P. J. Mulholland, Hydrometric and stream chemistry evidence of three storm flowpaths in Walker Branch Watershed, J. Hydrol., 1993, 151, 291–316.
61. B. Stewart, J. B. Shanley, J. W. Kirchner, D. Norris, T. Adler, C. Bristol, A. A. Harpold, J. N. Perdrial, D. M. Rizzo, G. Sterle, K. L. Underwood, H. Wen and L. Li, Streams as mirrors: reading subsurface water chemistry from stream chemistry, Water Resour. Res., 2022, 58, e2021WR029931.
62. B. R. Waterman, G. Alcantar, S. G. Thomas and M. F. Kirk, Spatiotemporal variation in runoff and baseflow in watersheds located across a regional precipitation gradient, J. Hydrol. Reg. Stud., 2022, 41, 101071.
63. M. Stoelzle, K. Stahl, A. Morhard and M. Weiler, Streamflow sensitivity to drought scenarios in catchments with different geology, Geophys. Res. Lett., 2014, 41, 6174–6183.
64. V. Acuña, T. Datry, J. Marshall, D. Barceló, C. N. Dahm, A. Ginebreda, G. McGregor, S. Sabater, K. Tockner and M. A. Palmer, Why Should We Care About Temporary Waterways?, Proc. Am. Assoc. Adv. Sci., 2014, 343, 1080–1081.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2va00191h

---